Use multiple ARIMA models to bootstrap a confidence interval for the total

bootstrapForecast.R

bootstrapForecast = function(..., newxreglist=NULL, inv.transform=NULL, runs=1000, n.ahead=90, CI.level = .95) {
  models = list(...)
  nummodels = length(models)
  sims = vapply(1:runs, function(i) {
    sim = vapply(1:nummodels, function(j) {
      if (is.data.frame(newxreglist)) newxreg = newxreglist
      else newxreg = newxreglist[[j]]
      if (!is.null(inv.transform)) {
        if(is.list(inv.transform)) inv = inv.transform[[j]]
        else inv = inv.transform
        inv(simulate(models[[j]], nsim=n.ahead, bootstrap=T, xreg=newxreg))
      }
      else simulate(models[[j]], nsim=n.ahead, bootstrap=T, xreg=newxreg)
    }, rep(0, n.ahead))
    rowSums(sim)
  }, rep(0, n.ahead))
  list(bootstrap_mean=apply(sims,1,mean),
       bootstrap_pointwise_lower = apply(sims,1,function(i) quantile(i, (1-CI.level)/2)),
       bootstrap_pointwise_upper = apply(sims,1,function(i) quantile(i, 1 - (1-CI.level)/2)),
       total_lower = quantile(colSums(sims), (1-CI.level)/2),
       total_upper = quantile(colSums(sims), 1 - (1-CI.level)/2))
}