R/rwrf.R

In twolodzko/boots: Nonparametric Bootstrap/Resampling-Based Time-Series Forecasting

#' Random walk (with drift) resampling-based forecast
#' 
#' @param x          a numeric vector or time series
#' @param h          number of periods for forecasting
#' @param R          number of bootstrap replications
#' @param level      confidence levels for prediction intervals
#' @param lambda     Box-Cox transformation parameter. Ignored
#'                   if \code{NULL}. Otherwise, forecasts
#'                   back-transformed via an inverse Box-Cox
#'                   transformation
#' 
#' @importFrom stats median sd
#' @importFrom stats ts time frequency tsp start end
#' @importFrom forecast BoxCox InvBoxCox
#' 
#' @export

rwrf <- function(x, h = 10L, R = 500L, level = c(80, 95)) {
  
  if (!is.simple.vector(x))
    stop("x is not a vector")

  seed <- .Random.seed
  xname <- deparse(substitute(x))
  method <- "Random walk (with drift) resampling-based forecast"
  
  x <- as.ts(x)
  origx <- x
  n <- length(x)
  f <- frequency(x)
  xtsp <- tsp(x)

  # sample
  
  dx <- diff(x)
  
  smpl <- repeatFun(R, function(i) {
    cumsum(c(x[n], sample(dx, h, replace = TRUE)))[-1L]
  }, workers = 1L)
  
  # prepare results

  mean <- colMeans(smpl)
  sd <- apply(smpl, 2L, sd)
  
  quant <- t(apply(smpl, 2L, quantile, c(0.5, 0.5-level/200, 0.5+level/200)))
  median <- ts(quant[, 1L], start = xtsp[2L] + 1/f, frequency = f)
  q <- length(level)   
  lower <- ts(quant[, 2:(q+1), drop = FALSE], start = xtsp[2L] + 1/f, frequency = f)
  upper <- ts(quant[, (q+2):(2*q+1), drop = FALSE], start = xtsp[2L] + 1/f, frequency = f)
  
  fcast <- ts(mean, start = xtsp[2L] + 1/f, frequency = f)
  out <- list(method = method, level = level, x = origx, xname = xname,
              mean = fcast, sd = sd, median = median, lower = lower, upper = upper,
              R = R, seed = seed)
  structure(out, class = c("forecast", "boots"))
}