diff --git a/DESCRIPTION b/DESCRIPTION
index 3f400b3..1f18a94 100644
--- a/DESCRIPTION
+++ b/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
diff --git a/NAMESPACE b/NAMESPACE
index e711666..2031efb 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/R/RcppExports.R b/R/RcppExports.R
index cf9c280..313348a 100644
--- a/R/RcppExports.R
+++ b/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)
}
diff --git a/R/irf.R b/R/irf.R
new file mode 100644
index 0000000..291944f
--- /dev/null
+++ b/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
+}
diff --git a/R/plots.R b/R/plots.R
index 102a445..862ef66 100644
--- a/R/plots.R
+++ b/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)
+}
diff --git a/R/utility_functions.R b/R/utility_functions.R
index f370a3d..83e03ab 100644
--- a/R/utility_functions.R
+++ b/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)
}
diff --git a/man/irf.Rd b/man/irf.Rd
new file mode 100644
index 0000000..dada9fc
--- /dev/null
+++ b/man/irf.Rd
@@ -0,0 +1,33 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/irf.R
+\name{irf}
+\alias{irf}
+\title{Impulse response functions}
+\usage{
+irf(x, shock, ahead = 8)
+}
+\arguments{
+\item{x}{An object of type \code{bayesianVARs_bvar}.}
+
+\item{shock}{A vector of shocks. If a factor model was used, then \code{shock} is
+applied to the factors, and its dimension must be the number of factors.
+
+\if{html}{\out{}}\preformatted{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.
+}\if{html}{\out{
}}}
+}
+\value{
+Returns a \code{bayesianVARs_irf} object
+}
+\description{
+Effect of a shock on the factors or variables over time
+}
+\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)
+
+}
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 98fffac..68cc4c5 100644
--- a/src/RcppExports.cpp
+++ b/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},
diff --git a/src/irf.cpp b/src/irf.cpp
new file mode 100644
index 0000000..ff1a3b6
--- /dev/null
+++ b/src/irf.cpp
@@ -0,0 +1,59 @@
+#include
+
+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;
+}