Add functions and tests for forecast and forecast errors

DESCRIPTION

@@ -7,11 +7,12 @@ Description: Identify VAR models by finding one shock that explains most of the
 License: MIT + file LICENSE
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 Imports: 
     expm,
     ggplot2,
     pracma,
+    purrr,
     stats,
     svars,
     vars

NAMESPACE

@@ -28,9 +28,11 @@ S3method(plot,fevdirf)
 S3method(plot,fevfd)
 S3method(plot,irffd)
 export(bootstrap)
+export(fe)
 export(fev)
 export(fevdfd)
 export(fevfd)
+export(forecast)
 export(hd)
 export(hs)
 export(id_fevdfd)

R/fe.R

@@ -0,0 +1,51 @@
+#' Calculate the historical forecast error for a VAR out to a specified horizon
+#'
+#' @param var vars::VAR or svars object
+#' @param horizon number of steps out to calculate fe
+#'
+#' @return forecast error
+#'
+#' @export
+#'
+#' @importFrom  expm %^%
+#'
+#' @examples
+#' x <- svars::USA
+#' v <- vars::VAR(x, p = 2)
+#' fe(v, 10)
+#'
+fe <- function(var, horizon = 1) {
+
+  ## Declare this to avoid "no visible binding" warning
+  `%^%` <- expm::`%^%`
+
+  if (inherits(var, "varest")) {
+    svar <- svars::id.chol(var)
+    res <- stats::residuals(var)
+  } else if (inherits(var, "svars")) {
+    svar <- var
+    res <- stats::residuals(svar$VAR)
+  }
+
+  ssv <- as_statespace_var(svar$A_hat, svar$B)
+  res <- cbind(res, matrix(0, nrow(res), var$K * (var$p - 1)))
+
+  fe <- matrix(NA, nrow(res), var$K)
+  names(fe) <- names(res)
+
+  for (t in seq_len(nrow(res) - (horizon - 1))) {
+    fe_t <- 0
+
+    for (h in seq_len(horizon)) {
+      irf_h <- ssv$mx %^% (h - 1)
+      res_h <- res[t + (horizon - 1) - (h - 1), ]
+
+      e <- irf_h %*% res_h
+
+      fe_t <- fe_t + e
+    }
+    fe[t + (horizon - 1), ] <- ssv$my %*% fe_t
+  }
+
+  return(fe)
+}

R/fev.R

@@ -26,10 +26,14 @@ fev <- function(var, n_ahead = 20) {
   ## Calculate IRFs out to horizon (then adj to 3-dim matrix from DF)
   if (inherits(var, "fevdvar")) {
     irf <- vars::irf(var, n.ahead = n_ahead, as_matrix = TRUE)
-  } else {
+  } else if (inherits(var, "svars")) {
     irf <- vars::irf(var, n.ahead = n_ahead)[[1]][, -1] |>
       apply(1, matrix, simplify = FALSE, nrow = k, ncol = k, byrow = TRUE) |>
       simplify2array()
+  } else if (inherits(var, "varest")) {
+    stop("Please pass an svar object instead of a varest.")
+  } else {
+    stop("This class of VAR not supported")
   }
 
   # Forecast Error Variance calculation

R/forecast.R

@@ -0,0 +1,58 @@
+#' Calculate the forecasts for a VAR out to a specified horizon
+#' Conditions on information at time "t"
+#'
+#' @param var vars::VAR or svars object
+#' @param horizon number of steps out to calculate forecast
+#'
+#' @return forecast
+#'
+#' @export
+#'
+#' @importFrom  expm %^%
+#'
+#' @examples
+#' x <- svars::USA
+#' v <- vars::VAR(x, p = 2)
+#' forecast(v, 10)
+#'
+forecast <- function(var, horizon = 1) {
+
+  ## Declare this to avoid "no visible binding" warning
+  `%^%` <- expm::`%^%`
+
+  if (inherits(var, "varest")) {
+    svar <- svars::id.chol(var)
+  } else if (inherits(var, "svars")) {
+    svar <- var
+  }
+
+  ssv <- as_statespace_var(svar$A_hat, svar$B)
+
+  x <- var$datamat[, 1:var$K]
+
+  x_lag1 <- var$datamat[1, (var$K + 1):(var$K * (var$p + 1))]
+  x_lag2 <- var$datamat[, 1:(var$K * var$p)]
+  names(x_lag1) <- names(x_lag2)
+
+  x_lag <- rbind(x_lag1, x_lag2)
+
+  if (var$type == "const") {
+    const_a <-
+      purrr::map(0:(horizon - 1), ~ ssv$mx %^% .x) |>
+      simplify2array() |>
+      rowSums(dims = 2)
+
+    const <- const_a %*% c(svar$A_hat[, "const"], rep(0, var$K * (var$p - 1)))
+  } else {
+    const <- matrix(0, var$K * var$p, 1)
+  }
+
+  non_const <- (ssv$mx %^% (horizon)) %*% t(as.matrix(x_lag))
+
+  x_hat_h <-  t(c(ssv$my %*% const) + ssv$my %*% non_const)
+
+  x_hat_h <- rbind(matrix(NA, var$p + (horizon - 1), var$K), x_hat_h)
+  colnames(x_hat_h) <- colnames(x)
+
+  return(x_hat_h)
+}

R/irf.R

@@ -151,9 +151,16 @@ irf.bvartools <- function(
 #' @param ssv Statespace VAR class.
 #' @param n_ahead integer. How far out to calculate the impulse response.
 #'
+#' @return array of irfs, dimensions = c(response, shock, horizon)
+#'
 irf_ssv <- function(
-    ssv,
-    n_ahead = 10) {
+  ssv,
+  n_ahead = 10
+) {
+
+  ## Declare this to avoid "no visible binding" warning
+  `%^%` <- expm::`%^%`
+
   k <- nrow(ssv$my)
 
   irf <- array(0, dim = c(k, k, n_ahead))

tests/testthat/test_fe.R

@@ -0,0 +1,27 @@
+test_fe_forecast <- function(x, var, horizon) {
+  forecast <- forecast(var, horizon = horizon)
+  fe_forecast <- rbind(x, matrix(NA, horizon, var$K)) - forecast
+
+  fe <- fe(var, horizon)
+
+  testthat::expect_equal(na.omit(fe_forecast), na.omit(fe), ignore_attr = TRUE, tolerance = 1e-5)
+}
+
+test_that("forecast matches predict function", {
+
+  x <- svars::USA
+  x2 <- svars::LN
+
+  var <- vars::VAR(x, p = 2)
+  var_const <- vars::VAR(x, p = 2, type = "const")
+  var_more_lags <- vars::VAR(x, p = 4, type = "const")
+  var_ln <- vars::VAR(x2, p = 4, type = "const")
+
+  for (h in c(1, 2, 3, 5, 10, 20, 50)) {
+    test_fe_forecast(x, var, h)
+    test_fe_forecast(x, var_const, h)
+    test_fe_forecast(x, var_more_lags, h)
+    test_fe_forecast(x2, var_ln, h)
+  }
+
+})

tests/testthat/test_fev.R

@@ -0,0 +1,37 @@
+test_fe_fev <- function(svar, horizon) {
+
+  fe <- fe(svar, horizon)
+  fe_fev_h <- diag(cov(na.omit(fe)))
+
+  fev_fev <- fev(svar, n_ahead = horizon + 1)$fev
+  fev_fev_h <- fev_fev[fev_fev$h == horizon, "fev"] |>
+    matrix(svar$K, svar$K) |>
+    colSums()
+
+
+  testthat::expect_equal(fe_fev_h, fev_fev_h, tolerance = 0.2)
+}
+
+test_that("forecast matches predict function", {
+
+  x <- svars::USA
+  x2 <- svars::LN
+
+  var <- vars::VAR(x, p = 2)
+  var_const <- vars::VAR(x, p = 2, type = "const")
+  var_more_lags <- vars::VAR(x, p = 4, type = "const")
+  var_ln <- vars::VAR(x2, p = 4, type = "const")
+
+  svar <- svars::id.chol(var)
+  svar_const <- svars::id.chol(var_const)
+  svar_more_lags <- svars::id.chol(var_more_lags)
+  svar_ln <- svars::id.chol(var_ln)
+
+  for (h in c(1, 2, 3, 5, 10, 20, 50)) {
+    test_fe_fev(svar, h)
+    test_fe_fev(svar_const, h)
+    test_fe_fev(svar_more_lags, h)
+    test_fe_fev(svar_ln, h)
+  }
+
+})

tests/testthat/test_forecast.R

@@ -0,0 +1,28 @@
+test_forecast_predict <- function(var, horizon) {
+  forecast <- forecast(var, horizon = horizon)
+  forecast <- forecast[nrow(forecast), ]
+  predict <- predict(var, n.ahead = horizon)$fcst |>
+    purrr::map(~.x[nrow(.x), "fcst"]) |>
+    simplify2array()
+
+  testthat::expect_equal(forecast, predict, ignore_attr = TRUE)
+}
+
+test_that("forecast matches predict function", {
+
+  x <- svars::USA
+  x2 <- svars::LN
+
+  var <- vars::VAR(x, p = 2)
+  var_const <- vars::VAR(x, p = 2, type = "const")
+  var_more_lags <- vars::VAR(x, p = 4, type = "const")
+  var_ln <- vars::VAR(x2, p = 4, type = "const")
+
+  for (h in c(1, 2, 3, 5, 10, 20, 50)) {
+    test_forecast_predict(var, h)
+    test_forecast_predict(var_const, h)
+    test_forecast_predict(var_more_lags, h)
+    test_forecast_predict(var_ln, h)
+  }
+
+})