R/bootstrap.R

In pli2016/forecast: Forecasting Functions for Time Series and Linear Models

#Bootstrap functions

# Trend estimation like STL without seasonality.
# Non-robust version
tl <- function(x, ...)
{
  x <- as.ts(x)
  tspx <- tsp(x)
  n <- length(x)
  tt <- 1:n
  fit <- supsmu(tt, x)
  out <- ts(cbind(trend=fit$y, remainder=x-fit$y))
  tsp(out) <- tsp(x)

  out <- structure(list(time.series=out),class="stl")
  return(out)
}

# Function to return some bootstrap samples of x
# based on LPB
lpb <- function(x, nsim=100)
{
  n <- length(x)
  meanx <- mean(x)
  y <- x - meanx
  gamma <- wacf(y, lag.max=n)$acf[,,1]
  s <- length(gamma)
  Gamma <- matrix(1, s, s)
  d <- row(Gamma) - col(Gamma)
  for(i in 1:(s-1))
    Gamma[d==i | d==(-i)] <- gamma[i+1]
  L <- t(chol(Gamma))
  W <- solve(L) %*% matrix(y,ncol=1)
  out <- ts(L %*% matrix(sample(W, n*nsim, replace=TRUE), nrow=n, ncol=nsim) + meanx)
  tsp(out) <- tsp(x)
  return(out)
}



# Bootstrapping time series (based on Bergmeir et al., 2016, IJF paper)
# Author: Fotios Petropoulos

MBB <- function(x, window_size) {
  
  bx = array(0, (floor(length(x)/window_size)+2)*window_size)
  for (i in 1:(floor(length(x)/window_size)+2)){
    c <- sample(1:(length(x)-window_size+1),1)
    bx[((i-1)*window_size+1):(i*window_size)] <- x[c:(c+window_size-1)]
  }
  start_from <- sample(0:(window_size-1),1) + 1
  bx[start_from:(start_from+length(x)-1)]
}


bld.mbb.bootstrap <- function(x, num, block_size = if(frequency(x)>1) 2*frequency(x) else 8 ) {

  freq <- frequency(x)
  
  xs <- list()
  xs[[1]] <- x # the first series is the original one
  
  if (num>1){
    # Box-Cox transformation
    lambda <- BoxCox.lambda(x, lower=0, upper=1)
    x.bc <- BoxCox(x, lambda)
    
    if (freq>1){
      # STL decomposition
      x.stl <- stl(ts(x.bc, frequency=freq), "per")$time.series
      seasonal <- x.stl[,1]
      trend <- x.stl[,2]
      remainder <- x.stl[,3]
    } else {
      # Loess
      trend <- 1:length(x) 
      suppressWarnings(x.loess <- loess(x.bc ~ trend, span=6/length(x), degree=1))
      seasonal <- rep(0, length(x))
      trend <- x.loess$fitted
      remainder <- x.loess$residuals
    }
    
    # Bootstrap some series, using MBB
    for (i in 2:num){      
      xs[[i]] <- InvBoxCox(trend + seasonal + MBB(remainder, block_size), lambda)
    }
  }
  
  xs
}