Added generate.AR() and forecast.AR(bootstrap = TRUE)

cc @robjhyndman

NAMESPACE

@@ -21,6 +21,7 @@ S3method(forecast,VAR)
 S3method(forecast,croston)
 S3method(forecast,fable_theta)
 S3method(forecast,model_mean)
+S3method(generate,AR)
 S3method(generate,ARIMA)
 S3method(generate,ETS)
 S3method(generate,NNETAR)

NEWS.md

@@ -5,8 +5,10 @@
 * Added the `selection_metric` argument to `ARIMA()`, which allows more control
   over the measure used to select the best model. By default this function will
   extract the information criteria specified by the `ic` argument.
-* Added `refit()` method for NNETAR (#287, @Tim-TU)
-* Added `refit()` method for MEAN (#289, @Tim-TU)
+* Added `refit()` method for NNETAR (#287, @Tim-TU).
+* Added `refit()` method for MEAN (#289, @Tim-TU).
+* Added `generate()` method for AR, the `forecast()` method now supports
+  bootstrap forecasting via this new method.
 
 ## Bug fixes
 

R/ar.R

@@ -186,13 +186,23 @@ estimate_ar <- function(x, p, xreg, constant, fixed) {
   E <- E * x_sd
   x <- x * x_sd
 
+  nx <- length(coef) - p
+  if (!is.null(xreg)) {
+    xcoef <- coef[seq_len(nx)]
+    xm <- drop(X[,seq_len(nx),drop=FALSE] %*% xcoef)
+  } else {
+    xm <- rep(0, nrow(new_data))
+  }
+  xm <- c(rep.int(NA_real_, p), xm)
+
   # Output model
   structure(
     list(
       coef = coef,
       coef.se = coef_se,
       fits = c(rep.int(NA_real_, p), YH),
       resid = c(rep.int(NA_real_, p), E),
+      reg_resid = x - xm*x_sd,
       last = x[(length(E)+1):length(x)],
       sigma2 = drop(varE),
       aic = aic,
@@ -212,9 +222,12 @@ estimate_ar <- function(x, p, xreg, constant, fixed) {
 #'   forecast()
 #' @export
 forecast.AR <- function(object, new_data = NULL, specials = NULL,
-                         bootstrap = FALSE, times = 5000, ...) {
+                        bootstrap = FALSE, times = 5000, ...) {
   if (bootstrap) {
-    abort("Bootstrapped forecasts for ARs are not yet implemented.")
+    sim <- map(seq_len(times), function(x) generate(object, new_data, specials, bootstrap = TRUE)[[".sim"]]) %>%
+      transpose() %>%
+      map(as.numeric)
+    return(distributional::dist_sample(sim))
   }
 
   h <- NROW(new_data)
@@ -248,6 +261,56 @@ forecast.AR <- function(object, new_data = NULL, specials = NULL,
   distributional::dist_normal(fc, se)
 }
 
+#' @inherit generate.ARIMA
+#' @examples
+#' as_tsibble(lh) %>%
+#'   model(AR(value ~ order(3))) %>%
+#'   generate()
+#' @export
+generate.AR <- function(object, new_data = NULL, specials = NULL, 
+                        bootstrap = FALSE, ...) {
+  n <- length(object$last)
+  p <- object$p
+  coef <- object$coef
+
+  # Get xreg
+  h <- max(map_int(key_data(new_data)[[".rows"]], length))
+  xreg <- specials$xreg[[1]]$xreg
+  if(object$constant){
+    xreg <- cbind(constant = rep(1, h), xreg)
+  }
+
+  # Predict xreg
+  nx <- length(coef) - p
+  if (!is.null(xreg)) {
+    xcoef <- coef[seq_len(nx)]
+    xm <- drop(xreg %*% xcoef)
+  } else {
+    xm <- rep(0, nrow(new_data))
+  }
+
+  ar <- coef[nx + seq_len(p)]
+  # Generate future innovations if missing
+  if(!(".innov" %in% names(new_data))){
+    if(bootstrap){
+      res <- residuals(x)
+      new_data$.innov <- sample(na.omit(res) - mean(res, na.rm = TRUE),
+                                nrow(new_data), replace = TRUE)
+    }
+    else{
+      new_data$.innov <- stats::rnorm(nrow(new_data), sd = sqrt(object$sigma2))
+    }
+  }
+
+  new_data %>% 
+    group_by_key() %>% 
+    transmute(".sim" := stats::filter(
+      !!sym(".innov"), ar, method = "recursive", 
+      init = rev(object$reg_resid)[seq_along(ar)])) %>% 
+    dplyr::ungroup() %>% 
+    mutate(".sim" := as.numeric(!!sym(".sim") + !!xm))
+}
+
 #' Refit an AR model
 #'
 #' Applies a fitted AR model to a new dataset.