Faster forecasts

R/RcppExports.R

@@ -17,6 +17,10 @@
     .Call(`_bgvars_bgvectvpalg`, object)
 }
 
+.draw_forecast <- function(i, k, a0, a, b_, c_, sigma, pred) {
+    .Call(`_bgvars_draw_forecast`, i, k, a0, a, b_, c_, sigma, pred)
+}
+
 .gir <- function(A, h, impulse, response) {
     .Call(`_bgvars_gir`, A, h, impulse, response)
 }

R/predict.bgvar.R

@@ -43,7 +43,7 @@
 #'                  country_data = country_data,
 #'                  global_data = global_data,
 #'                  countries = c("US", "JP", "CA", "NO", "GB", "EA"), 
-#'                  domestic = list(variables = c("y", "Dp", "r"), lags = 1),
+#'                  domestic = list(variables = c("y", "Dp", "r"), lags = 2),
 #'                  foreign = list(variables = c("y", "Dp", "r"), lags = 1),
 #'                  global = list(variables = c("poil"), lags = 1),
 #'                  deterministic = list(const = TRUE, trend = FALSE, seasonal = FALSE),
@@ -87,9 +87,6 @@
 #' @export
 predict.bgvar <- function(object, ..., n.ahead = 10, new_x = NULL, new_d = NULL, ci = .95) {
 
-  # Dev specs
-  # n.ahead = 10; new_x = NULL; new_d = rep(1, 10); ci = .95
-
   k <- length(object[["model"]][["endogen"]][["variables"]]) # Endogenous variables
   p <- object[["model"]][["endogen"]][["lags"]]  # Determine number of lags
   tot <- k * p # Total number of endogenous regressors
@@ -107,13 +104,9 @@ predict.bgvar <- function(object, ..., n.ahead = 10, new_x = NULL, new_d = NULL,
                                    exogen = object[["data"]][["global"]], s = s)
   tt <- nrow(model_data[["data"]][["Y"]]) # Number of observations
 
-  # Endogenous variables
-  a <- object[["a"]]
-
   # Global variables
   m <- 0
   if (global) {
-    a <- cbind(a, object[["b"]])
     n_glo <- k_glo * (s + 1)
     tot <- tot + n_glo
     if (is.null(new_x)) {
@@ -128,7 +121,6 @@ predict.bgvar <- function(object, ..., n.ahead = 10, new_x = NULL, new_d = NULL,
 
   # Determinisitc terms
   if ("c" %in% names(object)) {
-    a <- cbind(a, object[["c"]])
     n <- ncol(object[["c"]]) / k
     tot <- tot + n
     if (is.null(new_d)) {
@@ -170,27 +162,21 @@ predict.bgvar <- function(object, ..., n.ahead = 10, new_x = NULL, new_d = NULL,
     pred[pos_d, ] <- cbind(object[["data"]][["deterministic"]][tt, ], t(new_d))
   }
 
+  # Determine number of draws
+  draws <- NA
+  for (i in c("a0", "a", "b", "c", "sigma")) {
+    if (is.na(draws)) {
+      if (!is.null(object[[i]])) {
+        draws <- nrow(object[[i]])
+      }
+    }
+  }
+
   # Calculate forecasts
-  draws <- nrow(a)
   result <- array(NA, dim = c(k, n.ahead, draws))
   pb <- utils::txtProgressBar(style = 3)
   for (draw in 1:draws) {
-    for (i in 1:n.ahead) {
-
-      a0_i <- solve(matrix(object[["a0"]][draw, ], k))
-
-      # Generate random error
-      temp <- eigen(matrix(object[["sigma"]][draw, ], k))
-      u <- temp$vectors %*% diag(sqrt(temp$values), k) %*% t(temp$vectors) %*% stats::rnorm(k)
-
-      pred[1:k, i + 1] <- a0_i %*% matrix(a[draw, ], k) %*% pred[, i] + a0_i %*% u
-      if (p > 1) {
-        for (j in 1:(p - 1)) {
-          pred[j * k + 1:k, i + 1] <- pred[(j - 1) * k + 1:k, i]
-        } 
-      }
-    }
-    result[,, draw] <- pred[1:k, -1]
+    result[,, draw] <- .draw_forecast(draw, k, object[["a0"]], object[["a"]], object[["b"]], object[["c"]], object[["sigma"]], pred)[1:k, -1]
     utils::setTxtProgressBar(pb, draw / draws)
   }
   close(pb)

src/RcppExports.cpp

@@ -55,6 +55,24 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// draw_forecast
+arma::mat draw_forecast(int& i, int& k, arma::mat& a0, arma::mat& a, Rcpp::Nullable<Rcpp::NumericMatrix>& b_, Rcpp::Nullable<Rcpp::NumericMatrix>& c_, arma::mat& sigma, arma::mat pred);
+RcppExport SEXP _bgvars_draw_forecast(SEXP iSEXP, SEXP kSEXP, SEXP a0SEXP, SEXP aSEXP, SEXP b_SEXP, SEXP c_SEXP, SEXP sigmaSEXP, SEXP predSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< int& >::type i(iSEXP);
+    Rcpp::traits::input_parameter< int& >::type k(kSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type a0(a0SEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type a(aSEXP);
+    Rcpp::traits::input_parameter< Rcpp::Nullable<Rcpp::NumericMatrix>& >::type b_(b_SEXP);
+    Rcpp::traits::input_parameter< Rcpp::Nullable<Rcpp::NumericMatrix>& >::type c_(c_SEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type sigma(sigmaSEXP);
+    Rcpp::traits::input_parameter< arma::mat >::type pred(predSEXP);
+    rcpp_result_gen = Rcpp::wrap(draw_forecast(i, k, a0, a, b_, c_, sigma, pred));
+    return rcpp_result_gen;
+END_RCPP
+}
 // gir
 arma::mat gir(Rcpp::List A, int h, int impulse, int response);
 RcppExport SEXP _bgvars_gir(SEXP ASEXP, SEXP hSEXP, SEXP impulseSEXP, SEXP responseSEXP) {
@@ -103,6 +121,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_bgvars_bgvartvpalg", (DL_FUNC) &_bgvars_bgvartvpalg, 1},
     {"_bgvars_bgvecalg", (DL_FUNC) &_bgvars_bgvecalg, 1},
     {"_bgvars_bgvectvpalg", (DL_FUNC) &_bgvars_bgvectvpalg, 1},
+    {"_bgvars_draw_forecast", (DL_FUNC) &_bgvars_draw_forecast, 8},
     {"_bgvars_gir", (DL_FUNC) &_bgvars_gir, 4},
     {"_bgvars_ir", (DL_FUNC) &_bgvars_ir, 5},
     {"_bgvars_vardecomp", (DL_FUNC) &_bgvars_vardecomp, 3},

src/draw_forecast.cpp

@@ -0,0 +1,60 @@
+// [[Rcpp::depends(RcppArmadillo)]]
+
+#include <RcppArmadillo.h>
+
+// [[Rcpp::export(.draw_forecast)]]
+arma::mat draw_forecast(int &i, // index of draw
+                        int &k, // number of endogenous vars
+                        arma::mat &a0, // A0
+                        arma::mat &a, // A
+                        Rcpp::Nullable<Rcpp::NumericMatrix> &b_, // B
+                        Rcpp::Nullable<Rcpp::NumericMatrix> &c_, // C
+                        arma::mat &sigma, // Sigma
+                        arma::mat pred) { // Data matrix for prediction
+
+  const int n_a = a.n_cols / k;
+  const int p = n_a / k;
+  int n_b = 0, n_c = 0;
+  const int n_tot = pred.n_rows;
+  const int n_ahead = pred.n_cols - 1;
+  arma::mat a0_i = arma::zeros<arma::mat>(k, k);
+  arma::vec eigval, u;
+  arma::mat eigvec;
+
+  // Collect posterior draws in one coefficient matrix  
+  arma::mat coef = arma::zeros<arma::mat>(k, n_tot);
+  coef.submat(0, 0, k - 1, n_a - 1) = arma::reshape(a.row(i - 1), k, n_a);
+  if (b_.isNotNull()) {
+    arma::mat b = Rcpp::as<arma::mat>(b_).row(i - 1);
+    n_b = b.n_cols / k;
+    coef.submat(0, n_a, k - 1, n_a + n_b - 1) = arma::reshape(b, k, n_b);
+  }
+  if (c_.isNotNull()) {
+    arma::mat c = Rcpp::as<arma::mat>(c_).row(i - 1);
+    n_c = c.n_cols / k;
+    coef.submat(0, n_a + n_b, k - 1, n_a + n_b + n_c - 1) = arma::reshape(c, k, n_c);
+  }
+
+  // Invert A0
+  a0_i = arma::solve(arma::reshape(a0.row(i - 1), k, k), arma::eye<arma::mat>(k, k));
+
+  // Forecast iterations
+  for (int j = 0; j < n_ahead; j++) {
+
+    // Generate random error
+    arma::eig_sym(eigval, eigvec, arma::reshape(sigma.row(i - 1), k, k));
+    u = eigvec * arma::diagmat(arma::sqrt(eigval)) * eigvec.t() * arma::randn(k);
+
+    // Forecast for next period
+    pred.submat(0, j + 1, k - 1, j + 1) = a0_i * coef * pred.col(j) + a0_i * u;
+
+    // Update lags
+    if (p > 1) {
+      for (int l = 0; l < (p - 1); l++) {
+        pred.submat((l + 1) * k, j + 1, (l + 2) * k - 1, j + 1) = pred.submat(l * k, j, (l + 1) * k - 1, j); 
+      }
+    }
+  }
+
+  return pred;
+}