R/predict.R

In tsDyn: Nonlinear Time Series Models with Regime Switching

#' Predict method for objects of class \sQuote{\code{nlar}}.
#' 
#' Forecasting a non-linear model object of general class \sQuote{\code{nlar}},
#' including \sQuote{\code{setar}} and \sQuote{\code{star}}.
#' 
#' @aliases predict predict.nlar
#' @param object An object of class \sQuote{\code{nlar}}; generated by
#' \command{setar()} or \command{lstar()}. 
#' @param newdata Optional. A new data frame to predict from.
#' @param n.ahead An integer specifying the number of forecast steps.
#' @param type Type of forecasting method used. See details.
#' @param nboot The number of replications for type MC or bootstrap.
#' @param ci The forecast confidence interval (available only with types MC and
#' bootstrap).
#' @param block.size The block size when the block-bootstrap is used.
#' @param boot1Zero Whether the first innovation for MC/bootstrap should be set
#' to zero.
#' @param seed optional, the seed for the random generation.
#' @param \dots Further arguments passed to the internal \sQuote{\code{oneStep}} function. 
#' Mainly argument \sQuote{\code{thVar}} if an external threshold variable was provided
#' 
#' @details
#' The forecasts are obtained recursively from the estimated model.  Given that
#' the models are non-linear, ignoring the residuals in the 2- and more steps
#' ahead forecasts leads to biased forecasts (so-called naive).  Different
#' resampling methods, averaging \code{n.boot} times over future residuals, are
#' available:
#' 
#' \describe{ \item{naive}{No residuals} \item{MC}{Monte-Carlo method, where
#' residuals are taken from a normal distribution, with a standard deviation equal to the
#' residuals sd. } \item{bootstrap}{Residuals are resampled from the empirical
#' residuals from the model.} \item{block-bootstrap}{Same as bootstrap, but
#' residuals are resampled in block, with size \code{block.size}} }
#' 
#' The \code{MC} and \code{bootstrap} methods correspond to equations 3.90 and
#' 3.91 of Franses and van Dijk (2000, p. 121).  The bootstrap/MC is initiated
#' either from the first forecast, \code{n.ahead}=1 (set with \code{boot1zero}
#' to TRUE), or from the second only.
#' 
#' When the forecast method is based on resampling, forecast intervals are
#' available. These are obtained simply as empirical \code{ci} quantiles of the
#' resampled forecasts (cf Method 2 in Franses and van Dijk, 2000, p. 122).
#' 
#' @return A \sQuote{\code{ts}} object, or, in the case of MC/bootstrap, a list
#' containing the prediction (pred) and the forecast standard errors
#' (\code{se}).
#' @author Matthieu Stigler
#' @seealso The model fitting functions \code{\link{setar}},
#' \code{\link{lstar}}. 
#' 
#' A more sophisticated predict function, allowing to do sub-sample rolling
#' predictions: \code{\link{predict_rolling}}.
#' @references Non-linear time series models in empirical finance, Philip Hans
#' Franses and Dick van Dijk, Cambridge: Cambridge University Press (2000).
#' @keywords regression
#' @examples
#' 
#' 
#' x.train <- window(log10(lynx), end = 1924)
#' x.test <- window(log10(lynx), start = 1925)
#' 
#' 
#' ### Use different forecasting methods:
#' mod.set <- setar(x.train, m=2, thDelay=0)
#' pred_setar_naive <- predict(mod.set, n.ahead=10)
#' pred_setar_boot <- predict(mod.set, n.ahead=10, type="bootstrap", n.boot=200)
#' pred_setar_Bboot <- predict(mod.set, n.ahead=10, type="block-bootstrap", n.boot=200)
#' pred_setar_MC <- predict(mod.set, n.ahead=10, type="bootstrap", n.boot=200)
#' 
#' ## Plot to compare results:
#' pred_range <- range(pred_setar_naive, pred_setar_boot$pred, pred_setar_MC$pred, na.rm=TRUE)
#' plot(x.test, ylim=pred_range, main="Comparison of forecasts methods from same SETAR")
#' lines(pred_setar_naive, lty=2, col=2)
#' lines(pred_setar_boot$pred, lty=3, col=3)
#' lines(pred_setar_Bboot$pred, lty=4, col=4)
#' lines(pred_setar_MC$pred, lty=5, col=5)
#' 
#' legLabels <- c("Observed", "Naive F", "Bootstrap F","Block-Bootstrap F", "MC F")
#' legend("bottomleft", leg=legLabels, lty=1:5, col=1:5)
#' 

#' @export
#### Predict nlar
predict.nlar <- function(object, newdata, n.ahead=1, type=c("naive", "MC", "bootstrap", "block-bootstrap"),
                         nboot=100, ci=0.95, block.size=3, boot1Zero=TRUE, seed = NULL, ...)
{

  type <- match.arg(type)

  if(missing(newdata)) 
    newdata <- object$str$x

  res <- newdata
  n.used <- length(res)
  m <- object$str$m
  d <- object$str$d
  sd.res <- sqrt( mse(object) )
  steps <- object$str$steps
  resid <- as.numeric(residuals(object))
  resid <- resid[!is.na(resid)]

  tsp(res) <- NULL
  class(res) <- NULL

  res <- c(res, rep(0, n.ahead))
  xrange <- (m-1)*d + steps - ((m-1):0)*d

  if(type=="naive") nboot <-1

  ## Loop for new predictions using specific oneStep()
  pred_fun <- function(res){ 
    innov <- switch(type, 
                    "naive"= rep(0, n.ahead), 
                    "MC"= rnorm(n.ahead, mean=0, sd=sd.res), 
                    "bootstrap" = sample(resid, size=n.ahead, replace=TRUE),
                    "block-bootstrap" = sample.block(resid, block.size= block.size))
    if(boot1Zero) innov[1] <- 0
    for(i in n.used + (1:n.ahead)) {
      res[i] <- oneStep(object, newdata = t(as.matrix(res[i - xrange])),
			itime=(i - n.used), ...)
      res[i] <- res[i] + innov[i-n.used]
    }
    return(res)
  }

  if(!is.null(seed)) set.seed(seed)
  res <- replicate(nboot, pred_fun(res=res))
  res_means <- rowMeans(res, na.rm=TRUE)
  pred <- res_means[n.used + 1:n.ahead]

## Compute SE if MC/boot
  if(type!="naive"){
    ci_half <- (1-ci)/2
    SE <- t(apply(res[n.used + 1:n.ahead, ,drop=FALSE], 1 , sd, na.rm=TRUE))
    CI <- t(apply(res[n.used + 1:n.ahead, ,drop=FALSE], 1 , quantile, prob=sort(c(ci_half, 1-ci_half)), na.rm=TRUE))
    if(is.ts(newdata)){
      SE <- ts(SE, start = tsp(newdata)[2] + deltat(newdata),
               frequency=frequency(newdata))
      CI <- ts(CI, start = tsp(newdata)[2] + deltat(newdata),
               frequency=frequency(newdata))
    }
  }

## Return result
  if(is.ts(newdata)) pred <- ts(pred, start = tsp(newdata)[2] + deltat(newdata),
				frequency=frequency(newdata))
  if(type=="naive"){
    result <- pred
  } else {
    result <- list(pred=pred, se=SE, ci=CI)
  }

  return(result)
  
}


#### Predict TVECM/TVAR
#' @export
predict.TVECM <- function(object, ...) stop("predict() not available for TVECM")




############# TESTS
if(FALSE){
  library(tsDyn)
  
  set <- setar(lynx, m=2)
  tsDyn:::predict.nlar(set, n.ahead=1)
  
  environment(predict.nlar) <- environment(star)
  environment(oneStep) <- environment(star)
  
## n-ahead: 1
  tsDyn:::predict.nlar(set, n.ahead=1)
  tsDyn:::predict.nlar(set, n.ahead=1, type="naive")
  tsDyn:::predict.nlar(set, n.ahead=1, type="MC")
  tsDyn:::predict.nlar(set, n.ahead=1, type="bootstrap")
  
## n-ahead: 2
  tsDyn:::predict.nlar(set, n.ahead=2)
  tsDyn:::predict.nlar(set, n.ahead=2, type="naive")
  tsDyn:::predict.nlar(set, n.ahead=2, type="MC", nboot=100)
  tsDyn:::predict.nlar(set, n.ahead=2, type="bootstrap", nboot=100)
  
## n-ahead: 5
  tsDyn:::predict.nlar(set, n.ahead=5)
  tsDyn:::predict.nlar(set, n.ahead=5, type="naive")
  tsDyn:::predict.nlar(set, n.ahead=5, type="MC", nboot=5)
  
###### LSTAR
  mod.lst <- lstar(lynx, m=2, th=1300)
  
  predict(mod.lst, n.ahead=5)
  predict(mod.lst, n.ahead=5, type="MC", nboot=100)
  predict(mod.lst, n.ahead=5, type="bootstrap", nboot=100)
  
  ###### AAR
  mod.aar <- aar(lynx, m=2)
  predict(object=mod.aar , n.ahead=5)
  
  ## NNet
  mod.nnet <- nnetTs(log10(lynx), m=1, size=3, control=list(maxit=1000))
  
  predict(mod.nnet, n.ahead=10)
  predict(mod.nnet, n.ahead=10, type="MC", nboot=2)
  predict(mod.nnet, n.ahead=10, type="bootstrap", nboot=2)
  
  
}