diff --git a/DESCRIPTION b/DESCRIPTION index e364d28..d607a5a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -18,8 +18,13 @@ Imports: Rcpp (>= 1.0.12), RcppProgress, RcppTN -LinkingTo: bsvars, Rcpp, RcppArmadillo, RcppProgress, RcppTN +LinkingTo: + bsvars, + Rcpp, + RcppArmadillo, + RcppProgress, + RcppTN Suggests: tinytest LazyData: true Encoding: UTF-8 -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.2 diff --git a/NAMESPACE b/NAMESPACE index 7cc25aa..e6b5088 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,6 @@ # Generated by roxygen2: do not edit by hand +S3method(compute_variance_decompositions,PosteriorBVARPANEL) S3method(estimate,BVARPANEL) S3method(estimate,PosteriorBVARPANEL) S3method(forecast,PosteriorBVARPANEL) @@ -13,6 +14,7 @@ importFrom(R6,R6Class) importFrom(Rcpp,sourceCpp) importFrom(RcppTN,dtn) importFrom(RcppTN,rtn) +importFrom(bsvars,compute_variance_decompositions) importFrom(bsvars,estimate) importFrom(bsvars,forecast) useDynLib(bvarPANELs, .registration = TRUE) diff --git a/NEWS.md b/NEWS.md index b83ec97..1899ac9 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,4 +7,5 @@ 5. The package includes possibility of conditional forecasting given projected paths of some of the variables. [#8](https://github.com/bsvars/bvarPANELs/issues/8) 6. The package includes estimation and forecasting with exogenous variables. [#12](https://github.com/bsvars/bvarPANELs/issues/12) +7. The package includes forecast error variance decomposition analysis [#14](https://github.com/bsvars/bvarPANELs/issues/14) diff --git a/R/RcppExports.R b/R/RcppExports.R index 8e5f7af..8dbd966 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -5,6 +5,14 @@ bvarPANEL <- function(S, Y, X, prior, starting_values, thin, show_progress, rate .Call(`_bvarPANELs_bvarPANEL`, S, Y, X, prior, starting_values, thin, show_progress, rate_target_start_initial) } +Sigma2B_c <- function(posterior_Sigma_c, lower = TRUE) { + .Call(`_bvarPANELs_Sigma2B_c`, posterior_Sigma_c, lower) +} + +panel_variance_decompositions <- function(posterior_Sigma, posterior_A, global_Sigma, global_A, horizon, p, lower = TRUE) { + .Call(`_bvarPANELs_panel_variance_decompositions`, posterior_Sigma, posterior_A, global_Sigma, global_A, horizon, p, lower) +} + rmniw1 <- function(A, V, S, nu) { .Call(`_bvarPANELs_rmniw1`, A, V, S, nu) } diff --git a/R/bvarPANELs-package.R b/R/bvarPANELs-package.R index b67257d..2bc7e9f 100644 --- a/R/bvarPANELs-package.R +++ b/R/bvarPANELs-package.R @@ -30,7 +30,7 @@ #' @aliases bvarPANELs-package bvarPANELs #' @docType package #' @useDynLib bvarPANELs, .registration = TRUE -#' @importFrom bsvars estimate forecast +#' @importFrom bsvars estimate forecast compute_variance_decompositions #' @importFrom Rcpp sourceCpp #' @importFrom R6 R6Class #' @importFrom RcppTN rtn dtn diff --git a/R/compute_variance_decompositions.PosteriorBVARPANEL.R b/R/compute_variance_decompositions.PosteriorBVARPANEL.R new file mode 100644 index 0000000..61ad84d --- /dev/null +++ b/R/compute_variance_decompositions.PosteriorBVARPANEL.R @@ -0,0 +1,92 @@ + +#' @title Computes posterior draws of the forecast error variance decomposition +#' @description For each country, each of the draws from the posterior estimation +#' of the model is transformed into a draw from the posterior distribution of the forecast +#' error variance decomposition. +#' +#' @method compute_variance_decompositions PosteriorBVARPANEL +#' +#' @param posterior posterior estimation outcome - an object of class +#' \code{PosteriorBVARPANEL} obtained by running the \code{estimate} function. +#' +#' @param horizon a positive integer number denoting the forecast horizon for +#' the forecast error variance decompositions. +#' +#' @return An object of class \code{PosteriorFEVDPANEL}, that is, a list with +#' \code{C} elements containing \code{NxNx(horizon+1)xS} arrays of class +#' \code{PosteriorFEVD} with \code{S} draws of country-specific forecast error +#' variance decompositions. +#' +#' @seealso \code{\link{estimate.PosteriorBVARPANEL}} +#' +#' @author Tomasz Woźniak \email{wozniak.tom@pm.me} +#' +#' @references +#' Lütkepohl, H. (2017). Structural VAR Tools, Chapter 4, In: Structural vector autoregressive analysis. Cambridge University Press. +#' +#' @examples +#' # upload data +#' data(ilo_cubic_panel) +#' +#' # specify the model and set seed +#' set.seed(123) +#' specification = specify_bvarPANEL$new(ilo_cubic_panel, p = 1) +#' +#' # run the burn-in +#' burn_in = estimate(specification, 10) +#' +#' # estimate the model +#' posterior = estimate(burn_in, 20) +#' +#' # compute forecast error variance decomposition 4 years ahead +#' fevd = compute_variance_decompositions(posterior, horizon = 4) +#' +#' # workflow with the pipe |> +#' ############################################################ +#' set.seed(123) +#' ilo_cubic_panel |> +#' specify_bvarPANEL$new(p = 1) |> +#' estimate(S = 10) |> +#' estimate(S = 20) |> +#' compute_variance_decompositions(horizon = 4) -> fevd +#' +#' @export +compute_variance_decompositions.PosteriorBVARPANEL <- function(posterior, horizon) { + + posterior_Sigma = posterior$posterior$Sigma_c_cpp + posterior_A = posterior$posterior$A_c_cpp + posterior_Sg = posterior$posterior$Sigma + posterior_Ag = posterior$posterior$A + N = dim(posterior_Ag)[2] + C = dim(posterior_A[1][[1]])[3] + S = dim(posterior_A)[1] + p = posterior$last_draw$p + c_names = names(posterior$last_draw$data_matrices$Y) + + fff = .Call(`_bvarPANELs_panel_variance_decompositions`, + posterior_Sigma, + posterior_A, + posterior_Sg, + posterior_Ag, + horizon, + p, + TRUE + ) + + fevd = list() + for (c in 1:(C + 1)) { + fevd_c = array(NA, c(N, N, horizon + 1, S)) + for (s in 1:S) { + fevd_c[,,,s] = fff[c, s][[1]] + } + na_check = apply(fevd_c, 4, function(x) any(is.na(x))) + fevd_c = fevd_c[,,, !na_check] + class(fevd_c) = "PosteriorFEVD" + fevd[[c]] = fevd_c + } + + names(fevd) = c(c_names, "Global") + + class(fevd) <- "PosteriorFEVDPANEL" + return(fevd) +} \ No newline at end of file diff --git a/R/estimate.bvarPANEL.R b/R/estimate.bvarPANEL.R index 5ae49a0..0e737a9 100644 --- a/R/estimate.bvarPANEL.R +++ b/R/estimate.bvarPANEL.R @@ -83,7 +83,7 @@ estimate.BVARPANEL <- function( estimate.PosteriorBVARPANEL <- function( specification, S, - thin = 10, + thin = 1, show_progress = TRUE ) { diff --git a/inst/tinytest/test_compute_variance_decomposition.R b/inst/tinytest/test_compute_variance_decomposition.R new file mode 100644 index 0000000..da5f9cd --- /dev/null +++ b/inst/tinytest/test_compute_variance_decomposition.R @@ -0,0 +1,33 @@ + +data(ilo_cubic_panel) + +set.seed(1) +suppressMessages( + specification_no1 <- specify_bvarPANEL$new(ilo_cubic_panel) +) +run_no1 <- estimate(specification_no1, 3, 1, show_progress = FALSE) +fevd <- compute_variance_decompositions(run_no1, horizon = 2) + +set.seed(1) +suppressMessages( + fevd2 <- ilo_cubic_panel |> + specify_bvarPANEL$new() |> + estimate(S = 3, thin = 1, show_progress = FALSE) |> + compute_variance_decompositions(horizon = 2) +) + +expect_error( + compute_variance_decompositions(run_no1), + info = "compute_variance_decompositions: specify horizon." +) + +expect_equal( + sum(fevd$POL[1,,1,1]), 100, + info = "compute_variance_decompositions: sum to 100%." +) + +expect_identical( + fevd$POL[3,3,3,3], fevd2$POL[3,3,3,3], + info = "compute_variance_decompositions: identical for normal and pipe workflow." +) + diff --git a/man/compute_variance_decompositions.PosteriorBVARPANEL.Rd b/man/compute_variance_decompositions.PosteriorBVARPANEL.Rd new file mode 100644 index 0000000..fb50a74 --- /dev/null +++ b/man/compute_variance_decompositions.PosteriorBVARPANEL.Rd @@ -0,0 +1,63 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in +% R/compute_variance_decompositions.PosteriorBVARPANEL.R +\name{compute_variance_decompositions.PosteriorBVARPANEL} +\alias{compute_variance_decompositions.PosteriorBVARPANEL} +\title{Computes posterior draws of the forecast error variance decomposition} +\usage{ +\method{compute_variance_decompositions}{PosteriorBVARPANEL}(posterior, horizon) +} +\arguments{ +\item{posterior}{posterior estimation outcome - an object of class +\code{PosteriorBVARPANEL} obtained by running the \code{estimate} function.} + +\item{horizon}{a positive integer number denoting the forecast horizon for +the forecast error variance decompositions.} +} +\value{ +An object of class \code{PosteriorFEVDPANEL}, that is, a list with +\code{C} elements containing \code{NxNx(horizon+1)xS} arrays of class +\code{PosteriorFEVD} with \code{S} draws of country-specific forecast error +variance decompositions. +} +\description{ +For each country, each of the draws from the posterior estimation +of the model is transformed into a draw from the posterior distribution of the forecast +error variance decomposition. +} +\examples{ +# upload data +data(ilo_cubic_panel) + +# specify the model and set seed +set.seed(123) +specification = specify_bvarPANEL$new(ilo_cubic_panel, p = 1) + +# run the burn-in +burn_in = estimate(specification, 10) + +# estimate the model +posterior = estimate(burn_in, 20) + +# compute forecast error variance decomposition 4 years ahead +fevd = compute_variance_decompositions(posterior, horizon = 4) + +# workflow with the pipe |> +############################################################ +set.seed(123) +ilo_cubic_panel |> + specify_bvarPANEL$new(p = 1) |> + estimate(S = 10) |> + estimate(S = 20) |> + compute_variance_decompositions(horizon = 4) -> fevd + +} +\references{ +Lütkepohl, H. (2017). Structural VAR Tools, Chapter 4, In: Structural vector autoregressive analysis. Cambridge University Press. +} +\seealso{ +\code{\link{estimate.PosteriorBVARPANEL}} +} +\author{ +Tomasz Woźniak \email{wozniak.tom@pm.me} +} diff --git a/man/estimate.PosteriorBVARPANEL.Rd b/man/estimate.PosteriorBVARPANEL.Rd index 5e591b5..6e1e6f0 100644 --- a/man/estimate.PosteriorBVARPANEL.Rd +++ b/man/estimate.PosteriorBVARPANEL.Rd @@ -5,7 +5,7 @@ \title{Bayesian estimation of a Bayesian Hierarchical Panel Vector Autoregression using Gibbs sampler} \usage{ -\method{estimate}{PosteriorBVARPANEL}(specification, S, thin = 10, show_progress = TRUE) +\method{estimate}{PosteriorBVARPANEL}(specification, S, thin = 1, show_progress = TRUE) } \arguments{ \item{specification}{an object of class PosteriorBVARPANEL generated using diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 36f9d88..49d44a1 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -71,6 +71,35 @@ RcppExport SEXP _bvarPANELs_forecast_bvarPANEL(SEXP posterior_A_c_cppSEXP, SEXP UNPROTECT(1); return rcpp_result_gen; } +// Sigma2B_c +arma::cube Sigma2B_c(arma::cube& posterior_Sigma_c, const bool lower); +RcppExport SEXP _bvarPANELs_Sigma2B_c(SEXP posterior_Sigma_cSEXP, SEXP lowerSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::cube& >::type posterior_Sigma_c(posterior_Sigma_cSEXP); + Rcpp::traits::input_parameter< const bool >::type lower(lowerSEXP); + rcpp_result_gen = Rcpp::wrap(Sigma2B_c(posterior_Sigma_c, lower)); + return rcpp_result_gen; +END_RCPP +} +// panel_variance_decompositions +arma::field panel_variance_decompositions(arma::field& posterior_Sigma, arma::field& posterior_A, arma::cube& global_Sigma, arma::cube& global_A, const int horizon, const int p, const bool lower); +RcppExport SEXP _bvarPANELs_panel_variance_decompositions(SEXP posterior_SigmaSEXP, SEXP posterior_ASEXP, SEXP global_SigmaSEXP, SEXP global_ASEXP, SEXP horizonSEXP, SEXP pSEXP, SEXP lowerSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::field& >::type posterior_Sigma(posterior_SigmaSEXP); + Rcpp::traits::input_parameter< arma::field& >::type posterior_A(posterior_ASEXP); + Rcpp::traits::input_parameter< arma::cube& >::type global_Sigma(global_SigmaSEXP); + Rcpp::traits::input_parameter< arma::cube& >::type global_A(global_ASEXP); + Rcpp::traits::input_parameter< const int >::type horizon(horizonSEXP); + Rcpp::traits::input_parameter< const int >::type p(pSEXP); + Rcpp::traits::input_parameter< const bool >::type lower(lowerSEXP); + rcpp_result_gen = Rcpp::wrap(panel_variance_decompositions(posterior_Sigma, posterior_A, global_Sigma, global_A, horizon, p, lower)); + return rcpp_result_gen; +END_RCPP +} // rmniw1 arma::field rmniw1(const arma::mat& A, const arma::mat& V, const arma::mat& S, const double& nu); RcppExport SEXP _bvarPANELs_rmniw1(SEXP ASEXP, SEXP VSEXP, SEXP SSEXP, SEXP nuSEXP) { @@ -241,6 +270,8 @@ RcppExport SEXP _bvarPANELs_RcppExport_registerCCallable() { static const R_CallMethodDef CallEntries[] = { {"_bvarPANELs_bvarPANEL", (DL_FUNC) &_bvarPANELs_bvarPANEL, 8}, {"_bvarPANELs_forecast_bvarPANEL", (DL_FUNC) &_bvarPANELs_forecast_bvarPANEL, 6}, + {"_bvarPANELs_Sigma2B_c", (DL_FUNC) &_bvarPANELs_Sigma2B_c, 2}, + {"_bvarPANELs_panel_variance_decompositions", (DL_FUNC) &_bvarPANELs_panel_variance_decompositions, 7}, {"_bvarPANELs_rmniw1", (DL_FUNC) &_bvarPANELs_rmniw1, 4}, {"_bvarPANELs_sample_m", (DL_FUNC) &_bvarPANELs_sample_m, 5}, {"_bvarPANELs_sample_w", (DL_FUNC) &_bvarPANELs_sample_w, 2}, diff --git a/src/forecast_panel.cpp b/src/forecast_panel.cpp index defc831..bdb673b 100644 --- a/src/forecast_panel.cpp +++ b/src/forecast_panel.cpp @@ -17,8 +17,8 @@ Rcpp::List forecast_bvarPANEL ( const int horizon ) { - const int N = posterior_A_c_cpp(0).n_cols; const int S = posterior_A_c_cpp.n_elem; + const int N = posterior_A_c_cpp(0).n_cols; const int K = posterior_A_c_cpp(0).n_rows; const int C = posterior_A_c_cpp(0).n_slices; diff --git a/src/panel_variance_decompositions.cpp b/src/panel_variance_decompositions.cpp new file mode 100644 index 0000000..5dc453c --- /dev/null +++ b/src/panel_variance_decompositions.cpp @@ -0,0 +1,104 @@ + +#include +#include + +using namespace Rcpp; +using namespace arma; + + +// [[Rcpp:interface(cpp,r)]] +// [[Rcpp::export]] +arma::cube Sigma2B_c ( + arma::cube& posterior_Sigma_c, // (N, N, S) + const bool lower = true +) { + + const int S = posterior_Sigma_c.n_slices; + const int N = posterior_Sigma_c.n_rows; + + cube posterior_B_c(N, N, S); + + for (int s = 0; s < S; s++) { + if ( lower ) { + posterior_B_c.slice(s) = chol(posterior_Sigma_c.slice(s), "lower"); + } else { + posterior_B_c.slice(s) = chol(posterior_Sigma_c.slice(s), "upper"); + } + } // END s loop + + return posterior_B_c; +} // END Sigma2B_c + + + +arma::cube flip_cube_rows_cols ( + arma::cube& x // (N, K, S) +) { + + const int N = x.n_rows; + const int K = x.n_cols; + const int S = x.n_slices; + + cube y(K, N, S); + + for (int s = 0; s < S; s++) { + y.slice(s) = trans(x.slice(s)); + } + return y; +} // END flip_cube_rows_cols + + + +// [[Rcpp:interface(cpp,r)]] +// [[Rcpp::export]] +arma::field panel_variance_decompositions ( + arma::field& posterior_Sigma, // (S)(N, N, C) + arma::field& posterior_A, // (S)(K, N, C) + arma::cube& global_Sigma, // (N, N, S) + arma::cube& global_A, // (K, N, S) + const int horizon, + const int p, + const bool lower = true +) { + + const int S = posterior_Sigma.n_rows; + const int C = posterior_Sigma(0).n_slices; + + field fevds(C + 1, S); + + for (int s=0; s posterior_irf_s = bsvars::bsvars_ir ( + posterior_B_s, + A_s_aperm, + horizon, + p, + false + ); + + field posterior_fevd_s = bsvars::bsvars_fevd_homosk ( posterior_irf_s ); + for (int c = 0; c < C; c++) { + fevds(c,s) = posterior_fevd_s(c); + } // END s loop + } // END c loop + + // Global fevds + cube global_B = Sigma2B_c ( global_Sigma, lower ); + cube global_A_aperm = flip_cube_rows_cols ( global_A ); + field global_irf = bsvars::bsvars_ir ( + global_B, + global_A_aperm, + horizon, + p, + true + ); + field global_fevd = bsvars::bsvars_fevd_homosk ( global_irf ); + + for (int s=0; s + +using namespace Rcpp; +using namespace arma; + + +arma::cube Sigma2B_c ( + arma::cube& posterior_Sigma_c, // (N, N, S) + const bool lower = true +); + + +arma::cube flip_cube_rows_cols ( + arma::cube& x // (N, K, S) +); + + +arma::field panel_variance_decompositions ( + arma::field& posterior_Sigma, // (C)(N, N, S) + arma::field& posterior_A, // (C)(K, N, S) + arma::cube& global_Sigma, // (N, N, S) + arma::cube& global_A, // (K, N, S) + const int horizon, + const int p, + const bool lower = true +); + + +#endif // _PANEL_VARIANCE_DECOMPOSITIONS_H_