Impulse Reponse Functions

DESCRIPTION

@@ -49,4 +49,4 @@ Config/testthat/edition: 3
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2

NAMESPACE

@@ -7,6 +7,7 @@ S3method(fitted,bayesianVARs_bvar)
 S3method(pairs,bayesianVARs_predict)
 S3method(plot,bayesianVARs_bvar)
 S3method(plot,bayesianVARs_fitted)
+S3method(plot,bayesianVARs_irf)
 S3method(plot,bayesianVARs_predict)
 S3method(predict,bayesianVARs_bvar)
 S3method(print,bayesianVARs_bvar)
@@ -18,6 +19,7 @@ S3method(summary,bayesianVARs_draws)
 S3method(summary,bayesianVARs_predict)
 S3method(vcov,bayesianVARs_bvar)
 export(bvar)
+export(irf)
 export(my_gig)
 export(posterior_heatmap)
 export(specify_prior_phi)

R/RcppExports.R

@@ -23,6 +23,10 @@ my_gig <- function(n, lambda, chi, psi) {
     .Call(`_bayesianVARs_my_gig`, n, lambda, chi, psi)
 }
 
+irf_cpp <- function(coefficients, factor_loadings, U_vecs, shock, ahead) {
+    .Call(`_bayesianVARs_irf_cpp`, coefficients, factor_loadings, U_vecs, shock, ahead)
+}
+
 sample_PHI_cholesky <- function(PHI, PHI_prior, Y, X, U, d_sqrt, V_prior) {
     .Call(`_bayesianVARs_sample_PHI_cholesky`, PHI, PHI_prior, Y, X, U, d_sqrt, V_prior)
 }

R/irf.R

@@ -0,0 +1,50 @@
+#' Impulse response functions
+#'
+#' Effect of a shock on the factors or variables over time
+#'
+#' @param x An object of type `bayesianVARs_bvar`.
+#' @param shock A vector of shocks. If a factor model was used, then `shock` is
+#'		applied to the factors, and its dimension must be the number of factors.
+#'
+#'		If the Cholesky model was used, then the dimension of `shock` is expected
+#'		to be equals the number of equations. Note that the contemporaneous effects
+#'		of a shock in the Cholesky model depend on the ordering of equations.
+#'
+#' @return Returns a `bayesianVARs_irf` object
+#'
+#' @examples
+#' train_data <- 100 * usmacro_growth[,c("GDPC1", "PCECC96", "GPDIC1", "AWHMAN", "GDPCTPI", "CES2000000008x", "FEDFUNDS", "GS10", "EXUSUKx", "S&P 500")]
+#' prior_sigma <- specify_prior_sigma(data=train_data, type="cholesky")
+#' mod <- bvar(train_data, lags=2L, draws=2000, prior_sigma=prior_sigma)
+#' ir <- irf(mod, shock=c(0,0,0,0,0,0,1,0,0,0))
+#' plot(ir, n_col=2)
+#'
+#' @export
+irf <- function(x, shock, ahead=8) {
+	if (x$sigma_type == "factor") {
+		number_of_factors <- dim(x$facload)[2]
+		if (length(shock) != number_of_factors) {
+			stop("the shock vector must have dimension equals to the number of factors")
+		}
+	}
+	else if (x$sigma_type == "cholesky") {
+		if (length(shock) != ncol(mod$PHI)) {
+			stop("the shock vector must have dimension equals to the number of variables")
+		}
+	}
+	else {
+		stop("unknown model")
+	}
+
+	ret <- irf_cpp(
+		x$PHI,
+		x$facload,
+		x$U,
+		shock,
+		ahead
+	)
+
+	colnames(ret) <- colnames(x$Y)
+	class(ret) <- "bayesianVARs_irf"
+	ret
+}

R/plots.R

@@ -623,3 +623,78 @@ plot.bayesianVARs_predict <- function(x, dates = NULL, vars = "all", ahead = NUL
 
   invisible(x)
 }
+
+#' @export
+plot.bayesianVARs_irf <- function(x, vars = "all",
+                                      quantiles = c(0.05,0.25,0.5,0.75,0.95),
+                                      n_col = 1L,
+                                      ...){
+
+  n_ahead <- nrow(x)
+
+  if(length(vars)==1L & any(vars == "all")){
+    vars <- 1:ncol(x)
+  }else if(is.character(vars)){
+    if(any(base::isFALSE(vars %in% colnames(x)))){
+      stop("Elements of 'vars' must coincide with 'colnames(x)'!")
+    }else{
+      vars <- which(colnames(x) %in% vars)
+    }
+  }else if(is.numeric(vars)){
+    vars <- as.integer(vars)
+    if(any(vars > ncol(x))){
+      stop("'max(vars)' must be at most 'ncol(x)'!")
+    }
+  }
+  var_names <- colnames(x)
+
+  t <- seq_len(n_ahead)
+  dates <- t-1
+
+  quantiles <- sort(quantiles)
+  nr_quantiles <- length(quantiles)
+  nr_intervals <- floor(nr_quantiles/2)
+  even <- nr_quantiles%%2 == 0
+
+  pred_quants <- apply(x, 1:2, quantile, quantiles)
+
+  oldpar <- par(no.readonly = TRUE)
+  on.exit(par(oldpar), add = TRUE)
+  M <- length(vars)
+  par(mfrow=c(min(5,ceiling(M/n_col)),n_col), mar = c(2,2,2,1), mgp = c(2,.5,0))
+  for(j in seq_along(vars)){
+    plot(t, rep(0, n_ahead), type = "n", xlab="", ylab = "",
+         xaxt="n", ylim = range(pred_quants[,,vars[j]]))
+
+    axis(side=1, at = t, labels = dates[t])
+    mtext(var_names[j], side = 3)
+
+    if(nr_intervals>0){
+      for(r in seq.int(nr_intervals)){
+        alpha_upper <- 0.4
+        alpha_lower <- 0.2
+        alphas <- seq(alpha_lower, alpha_upper, length.out = nr_intervals)
+        if(nr_intervals==1){
+          alphas <- alpha_upper
+        }
+        polygon(x = c(t, rev(t)),
+                y = c(pred_quants[r,,vars[j]], rev(pred_quants[r+1,,vars[j]])),
+                col = scales::alpha("red", alphas[r]),
+                border = NA)
+        if(length(quantiles)>2){
+          polygon(x = c(t, rev(t)),
+                  y = c(pred_quants[nrow(pred_quants)+1-r,,vars[j]],
+                        rev(pred_quants[nrow(pred_quants)-r,,vars[j]])),
+                  col = scales::alpha("red", alphas[r]),
+                  border = NA)
+        }
+      }
+      if(!even){
+        lines(t, pred_quants[ceiling(length(quantiles)/2),,vars[j]],
+              col = "red", lwd = 2)
+      }
+    }
+  }
+
+  invisible(x)
+}

R/utility_functions.R

@@ -1733,7 +1733,6 @@ predict.bayesianVARs_bvar <- function(object, ahead = 1L, each = 1L, stable = TR
 
   # relevant mod settings
   sv_indicator <- which(object$heteroscedastic==TRUE)
-  intercept <- object$intercept
 
   # data preparation
   variables <- colnames(object$Y)
@@ -2038,10 +2037,10 @@ predict_old <- function(object, n.ahead, stable = TRUE, LPL = FALSE, Y_obs = NA,
 #' predictions <- predict(mod, ahead = 1L)
 #' print(predictions)
 print.bayesianVARs_predict <- function(x, ...){
-  cat(paste("\nGeneric functions for bayesianVARs_predict objects:\n",
-            " - summary.bayesianVARs_predict(),\n",
-            " - pairs.bayesianVARs_predict(),\n",
-            " - plot.bayesianVARs_predict() (alias for pairs.bayesianVARs_predict()).\n"))
+  cat("\nGeneric functions for bayesianVARs_predict objects:\n",
+      " - summary.bayesianVARs_predict(),\n",
+      " - pairs.bayesianVARs_predict(),\n",
+      " - plot.bayesianVARs_predict() (alias for pairs.bayesianVARs_predict()).\n")
   invisible(x)
 }
 

src/RcppExports.cpp

@@ -56,6 +56,21 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// irf_cpp
+arma::cube irf_cpp(const arma::cube& coefficients, const arma::cube& factor_loadings, const arma::mat& U_vecs, const arma::colvec& shock, const arma::uword ahead);
+RcppExport SEXP _bayesianVARs_irf_cpp(SEXP coefficientsSEXP, SEXP factor_loadingsSEXP, SEXP U_vecsSEXP, SEXP shockSEXP, SEXP aheadSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const arma::cube& >::type coefficients(coefficientsSEXP);
+    Rcpp::traits::input_parameter< const arma::cube& >::type factor_loadings(factor_loadingsSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type U_vecs(U_vecsSEXP);
+    Rcpp::traits::input_parameter< const arma::colvec& >::type shock(shockSEXP);
+    Rcpp::traits::input_parameter< const arma::uword >::type ahead(aheadSEXP);
+    rcpp_result_gen = Rcpp::wrap(irf_cpp(coefficients, factor_loadings, U_vecs, shock, ahead));
+    return rcpp_result_gen;
+END_RCPP
+}
 // sample_PHI_cholesky
 arma::mat sample_PHI_cholesky(const arma::mat PHI, const arma::mat& PHI_prior, const arma::mat& Y, const arma::mat& X, const arma::mat& U, const arma::mat& d_sqrt, const arma::mat& V_prior);
 RcppExport SEXP _bayesianVARs_sample_PHI_cholesky(SEXP PHISEXP, SEXP PHI_priorSEXP, SEXP YSEXP, SEXP XSEXP, SEXP USEXP, SEXP d_sqrtSEXP, SEXP V_priorSEXP) {
@@ -167,6 +182,7 @@ END_RCPP
 static const R_CallMethodDef CallEntries[] = {
     {"_bayesianVARs_bvar_cpp", (DL_FUNC) &_bayesianVARs_bvar_cpp, 21},
     {"_bayesianVARs_my_gig", (DL_FUNC) &_bayesianVARs_my_gig, 4},
+    {"_bayesianVARs_irf_cpp", (DL_FUNC) &_bayesianVARs_irf_cpp, 5},
     {"_bayesianVARs_sample_PHI_cholesky", (DL_FUNC) &_bayesianVARs_sample_PHI_cholesky, 7},
     {"_bayesianVARs_dmvnrm_arma_fast", (DL_FUNC) &_bayesianVARs_dmvnrm_arma_fast, 4},
     {"_bayesianVARs_predh", (DL_FUNC) &_bayesianVARs_predh, 6},

src/irf.cpp

@@ -0,0 +1,59 @@
+#include <RcppArmadillo.h>
+
+using namespace Rcpp;
+using namespace arma;
+
+inline void shift_and_insert(
+	arma::mat& X, //the columns of X should be y1,y2,y3, y1.l1,y2.l1,y3.l1,...,1
+	const arma::mat& new_y //what to insert in y1,y2,y3
+) {
+	for (uword i = X.n_cols-2; new_y.n_cols <= i; i--) {
+		X.col(i) = X.col(i-new_y.n_cols);
+	}
+	X.head_cols(new_y.n_cols) = new_y;
+}
+
+// [[Rcpp::export]]
+arma::cube irf_cpp(
+	const arma::cube& coefficients, //rows: lagged variables + intercept, columns: variables, slices: draws
+	const arma::cube& factor_loadings, //rows: variables, columns: factors, slices: draws,
+	const arma::mat& U_vecs, //rows: entries of a (variables x variables)-upper triagonal matrix with ones on the diagonals, cols: draws
+	const arma::colvec& shock, //columns: how much each of the factors is shocked
+	const arma::uword ahead //how far to predict ahead
+) {
+	const uword n_variables = coefficients.n_cols;
+	const uword n_posterior_draws = coefficients.n_slices;
+	const bool is_factor_model = factor_loadings.n_cols > 0;
+
+	arma::cube irf(ahead+1, n_variables, n_posterior_draws, arma::fill::none);
+
+	// trace out n_posterior_draws paths
+	arma::mat current_predictors(n_posterior_draws, coefficients.n_rows, arma::fill::zeros);
+	// all paths start with the some shock at t=0 to the variables
+	const arma::uvec upper_indices = arma::trimatu_ind(arma::size(n_variables, n_variables), 1);
+	for (uword r = 0; r < n_posterior_draws; r++) {
+		rowvec y_shock;
+		if (is_factor_model) {
+			y_shock = (factor_loadings.slice(r) * shock).t();
+		}
+		else {
+		    arma::mat U(n_variables, n_variables, arma::fill::eye);
+			U(upper_indices) = U_vecs.col(r);
+			U = U.t();
+			y_shock = solve(arma::trimatl(U), shock).t();
+		}
+
+		irf.slice(r).row(0) = y_shock;
+		current_predictors.head_cols(n_variables).row(r) = y_shock;
+	}
+
+	for (uword t = 1; t <= ahead; t++) {
+		for(uword r = 0; r < n_posterior_draws; r++) {
+			irf.slice(r).row(t) = current_predictors.row(r) * coefficients.slice(r);
+		}
+		// shift everything and make predictions the new predictors at lag zero
+		shift_and_insert(current_predictors, irf.row_as_mat(t));
+	}
+
+	return irf;
+}