Merge pull request #15 from bsvars/14-develop-method-compute_variance…

…_decompositionsbvarpanel

forecast error variance decompositions

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

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)

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)
 

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)
 }

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

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)
+}

R/estimate.bvarPANEL.R

@@ -83,7 +83,7 @@ estimate.BVARPANEL <- function(
 estimate.PosteriorBVARPANEL <- function(
     specification, 
     S, 
-    thin = 10, 
+    thin = 1, 
     show_progress = TRUE
 ) {
 

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."
+)
+

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<arma::cube> panel_variance_decompositions(arma::field<arma::cube>& posterior_Sigma, arma::field<arma::cube>& 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<arma::cube>& >::type posterior_Sigma(posterior_SigmaSEXP);
+    Rcpp::traits::input_parameter< arma::field<arma::cube>& >::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<arma::mat> 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},

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;