included compute_historical_decompositions #25

+ cpp functions and wrappers + R method + tests + doc + removed bsvarTOOLs.cpp and bsvarTOOLs.h and adjusted other cppand h files for that
bsvars · Jul 16, 2024 · 6c8779b · 6c8779b
1 parent 5bb6248
commit 6c8779b
Show file tree

Hide file tree

Showing 20 changed files with 357 additions and 164 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(compute_fitted_values,PosteriorBSVARSIGN)
+S3method(compute_historical_decompositions,PosteriorBSVARSIGN)
 S3method(compute_impulse_responses,PosteriorBSVARSIGN)
 S3method(compute_structural_shocks,PosteriorBSVARSIGN)
 S3method(estimate,BSVARSIGN)

diff --git a/R/compute.R b/R/compute.R
@@ -124,14 +124,14 @@ compute_fitted_values.PosteriorBSVARSIGN <- function(posterior) {
 
 
 
+#' @method compute_impulse_responses PosteriorBSVARSIGN
+#' 
 #' @title Computes posterior draws of impulse responses 
 #'
 #' @description Each of the draws from the posterior estimation of models from 
 #' packages \pkg{bsvars} or \pkg{bsvarSIGNs} is transformed into
 #' a draw from the posterior distribution of the impulse responses. 
 #' 
-#' @method compute_impulse_responses PosteriorBSVARSIGN
-#' 
 #' @param posterior posterior estimation outcome - an object of class 
 #' \code{PosteriorBSVARSIGN} obtained by running the \code{estimate} function.
 #' @param horizon a positive integer number denoting the forecast horizon for the impulse responses computations.
@@ -191,3 +191,82 @@ compute_impulse_responses.PosteriorBSVARSIGN <- function(posterior, horizon, sta
 
   return(irfs)
 } # END compute_impulse_responses.PosteriorBSVARSIGN
+
+
+
+
+#' @method compute_historical_decompositions PosteriorBSVARSIGN
+#' 
+#' @title Computes posterior draws of historical decompositions
+#'
+#' @description Each of the draws from the posterior estimation of models from
+#' packages \pkg{bsvars} or \pkg{bsvarSIGNs} is transformed into
+#' a draw from the posterior distribution of the historical decompositions. 
+#' IMPORTANT! The historical decompositions are interpreted correctly for 
+#' covariance stationary data. Application to unit-root non-stationary data might
+#' result in non-interpretable outcomes.
+#' 
+#' @param posterior posterior estimation outcome - an object of class 
+#' \code{PosteriorBSVARSIGN} obtained by running the \code{estimate} function.
+#' @param show_progress a logical value, if \code{TRUE} the estimation progress bar is visible
+#' 
+#' @return An object of class \code{PosteriorHD}, that is, an \code{NxNxTxS} array 
+#' with attribute \code{PosteriorHD} containing \code{S} draws of the historical 
+#' decompositions.
+#'
+#' @seealso \code{\link{estimate}}, \code{\link{summary}}, \code{\link{plot}}
+#'
+#' @author Xiaolei Wang \email{[email protected]} and Tomasz Woźniak \email{[email protected]}
+#' 
+#' @references 
+#' Kilian, L., & Lütkepohl, H. (2017). Structural VAR Tools, Chapter 4, In: Structural vector autoregressive analysis. Cambridge University Press.
+#' 
+#' @examples
+#' # upload data
+#' data(oil)
+#' 
+#' # specify the model and set seed
+#' set.seed(123)
+#' sign_irf       = array(matrix(c(-1, -1, 1, rep(0, 6)), nrow = 3), dim = c(3, 3, 1))
+#' specification  = specify_bsvarSIGN$new(oil, sign_irf = sign_irf)
+#' 
+#' # run the burn-in
+#' posterior      = estimate(specification, 10)
+#' 
+#' # compute historical decompositions
+#' hd            = compute_historical_decompositions(posterior)
+#' 
+#' # workflow with the pipe |>
+#' ############################################################
+#' set.seed(123)
+#' oil |>
+#'   specify_bsvarSIGN$new(sign_irf = sign_irf) |> 
+#'   estimate(S = 10) |> 
+#'   compute_historical_decompositions() -> hd
+#'   
+#' @export
+compute_historical_decompositions.PosteriorBSVARSIGN <- function(posterior, show_progress = TRUE) {
+
+  posterior_Theta0  = posterior$posterior$Theta0
+  posterior_B       = posterior$posterior$B
+  posterior_A       = posterior$posterior$A
+  posterior_At      = aperm(posterior_A, c(2, 1, 3))
+
+  Y                 = posterior$last_draw$data_matrices$Y
+  X                 = posterior$last_draw$data_matrices$X
+
+  N                 = nrow(Y)
+  T                 = ncol(Y)
+  p                 = posterior$last_draw$p
+  S                 = dim(posterior_A)[3]
+
+  ss                = .Call(`_bsvarSIGNs_bsvarSIGNs_structural_shocks`, posterior_B, posterior_A, Y, X)
+  ir                = .Call(`_bsvarSIGNs_bsvarSIGNs_ir`, posterior_At, posterior_Theta0, T, p, TRUE)
+  qqq               = .Call(`_bsvarSIGNs_bsvarSIGNs_hd`, ir, ss, show_progress)
+
+  hd                = array(NA, c(N, N, T, S))
+  for (s in 1:S) hd[,,,s] = qqq[s][[1]]
+  class(hd)         = "PosteriorHD"
+
+  return(hd)
+} # END compute_historical_decompositions.PosteriorBSVARSIGN
diff --git a/inst/include/bsvarSIGNs_RcppExports.h b/inst/include/bsvarSIGNs_RcppExports.h
@@ -25,27 +25,6 @@ namespace bsvarSIGNs {
         }
     }
 
-    inline arma::mat hd1_cpp(const int& var_i, const int& t, const int& h, const arma::mat& Epsilon, const arma::cube& irf) {
-        typedef SEXP(*Ptr_hd1_cpp)(SEXP,SEXP,SEXP,SEXP,SEXP);
-        static Ptr_hd1_cpp p_hd1_cpp = NULL;
-        if (p_hd1_cpp == NULL) {
-            validateSignature("arma::mat(*hd1_cpp)(const int&,const int&,const int&,const arma::mat&,const arma::cube&)");
-            p_hd1_cpp = (Ptr_hd1_cpp)R_GetCCallable("bsvarSIGNs", "_bsvarSIGNs_hd1_cpp");
-        }
-        RObject rcpp_result_gen;
-        {
-            RNGScope RCPP_rngScope_gen;
-            rcpp_result_gen = p_hd1_cpp(Shield<SEXP>(Rcpp::wrap(var_i)), Shield<SEXP>(Rcpp::wrap(t)), Shield<SEXP>(Rcpp::wrap(h)), Shield<SEXP>(Rcpp::wrap(Epsilon)), Shield<SEXP>(Rcpp::wrap(irf)));
-        }
-        if (rcpp_result_gen.inherits("interrupted-error"))
-            throw Rcpp::internal::InterruptedException();
-        if (Rcpp::internal::isLongjumpSentinel(rcpp_result_gen))
-            throw Rcpp::LongjumpException(rcpp_result_gen);
-        if (rcpp_result_gen.inherits("try-error"))
-            throw Rcpp::exception(Rcpp::as<std::string>(rcpp_result_gen).c_str());
-        return Rcpp::as<arma::mat >(rcpp_result_gen);
-    }
-
     inline Rcpp::List bsvar_sign_cpp(const int& S, const int& p, const arma::mat& Y, const arma::mat& X, const arma::field<arma::mat>& VB, const arma::cube& sign_irf, const arma::mat& sign_narrative, const arma::mat& sign_B, const arma::field<arma::mat>& Z, const Rcpp::List& prior, const Rcpp::List& starting_values, const bool show_progress = true, const int thin = 100, const int& max_tries = 10000) {
         typedef SEXP(*Ptr_bsvar_sign_cpp)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP);
         static Ptr_bsvar_sign_cpp p_bsvar_sign_cpp = NULL;
@@ -151,6 +130,48 @@ namespace bsvarSIGNs {
         return Rcpp::as<arma::field<arma::cube> >(rcpp_result_gen);
     }
 
+    inline arma::field<arma::cube> bsvarSIGNs_hd(arma::field<arma::cube>& posterior_irf_T, arma::cube& structural_shocks, const bool show_progress = true) {
+        typedef SEXP(*Ptr_bsvarSIGNs_hd)(SEXP,SEXP,SEXP);
+        static Ptr_bsvarSIGNs_hd p_bsvarSIGNs_hd = NULL;
+        if (p_bsvarSIGNs_hd == NULL) {
+            validateSignature("arma::field<arma::cube>(*bsvarSIGNs_hd)(arma::field<arma::cube>&,arma::cube&,const bool)");
+            p_bsvarSIGNs_hd = (Ptr_bsvarSIGNs_hd)R_GetCCallable("bsvarSIGNs", "_bsvarSIGNs_bsvarSIGNs_hd");
+        }
+        RObject rcpp_result_gen;
+        {
+            RNGScope RCPP_rngScope_gen;
+            rcpp_result_gen = p_bsvarSIGNs_hd(Shield<SEXP>(Rcpp::wrap(posterior_irf_T)), Shield<SEXP>(Rcpp::wrap(structural_shocks)), Shield<SEXP>(Rcpp::wrap(show_progress)));
+        }
+        if (rcpp_result_gen.inherits("interrupted-error"))
+            throw Rcpp::internal::InterruptedException();
+        if (Rcpp::internal::isLongjumpSentinel(rcpp_result_gen))
+            throw Rcpp::LongjumpException(rcpp_result_gen);
+        if (rcpp_result_gen.inherits("try-error"))
+            throw Rcpp::exception(Rcpp::as<std::string>(rcpp_result_gen).c_str());
+        return Rcpp::as<arma::field<arma::cube> >(rcpp_result_gen);
+    }
+
+    inline arma::mat hd1_cpp(const int& var_i, const int& t, const int& h, const arma::mat& Epsilon, const arma::cube& irf) {
+        typedef SEXP(*Ptr_hd1_cpp)(SEXP,SEXP,SEXP,SEXP,SEXP);
+        static Ptr_hd1_cpp p_hd1_cpp = NULL;
+        if (p_hd1_cpp == NULL) {
+            validateSignature("arma::mat(*hd1_cpp)(const int&,const int&,const int&,const arma::mat&,const arma::cube&)");
+            p_hd1_cpp = (Ptr_hd1_cpp)R_GetCCallable("bsvarSIGNs", "_bsvarSIGNs_hd1_cpp");
+        }
+        RObject rcpp_result_gen;
+        {
+            RNGScope RCPP_rngScope_gen;
+            rcpp_result_gen = p_hd1_cpp(Shield<SEXP>(Rcpp::wrap(var_i)), Shield<SEXP>(Rcpp::wrap(t)), Shield<SEXP>(Rcpp::wrap(h)), Shield<SEXP>(Rcpp::wrap(Epsilon)), Shield<SEXP>(Rcpp::wrap(irf)));
+        }
+        if (rcpp_result_gen.inherits("interrupted-error"))
+            throw Rcpp::internal::InterruptedException();
+        if (Rcpp::internal::isLongjumpSentinel(rcpp_result_gen))
+            throw Rcpp::LongjumpException(rcpp_result_gen);
+        if (rcpp_result_gen.inherits("try-error"))
+            throw Rcpp::exception(Rcpp::as<std::string>(rcpp_result_gen).c_str());
+        return Rcpp::as<arma::mat >(rcpp_result_gen);
+    }
+
     inline arma::cube forecast_bsvarSIGNs(arma::cube& posterior_Sigma, arma::cube& posterior_A, arma::vec& X_T, arma::mat& exogenous_forecast, arma::mat& cond_forecast, const int& horizon) {
         typedef SEXP(*Ptr_forecast_bsvarSIGNs)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP);
         static Ptr_forecast_bsvarSIGNs p_forecast_bsvarSIGNs = NULL;

diff --git a/inst/tinytest/test_compute_historical_decompositions.R b/inst/tinytest/test_compute_historical_decompositions.R
@@ -0,0 +1,29 @@
+
+data(oil)
+
+set.seed(1)
+suppressMessages(
+  specification_no1 <- specify_bsvarSIGN$new(oil)
+)
+run_no1             <- estimate(specification_no1, 3, 1, show_progress = FALSE)
+hd                  <- compute_historical_decompositions(run_no1)
+
+set.seed(1)
+suppressMessages(
+  hd2               <- oil |>
+    specify_bsvarSIGN$new() |>
+    estimate(S = 3, thin = 1, show_progress = FALSE) |>
+    compute_historical_decompositions()
+)
+
+
+
+expect_equal(
+  length(dim(hd)), length(dim(hd2)),
+  info = "compute_historical_decompositions: same output dimensions for normal and pipe workflow."
+)
+
+expect_identical(
+  hd[3,3,3,3], hd2[3,3,3,3],
+  info = "compute_historical_decompositions: identical for normal and pipe workflow."
+)
diff --git a/man/compute_historical_decompositions.PosteriorBSVARSIGN.Rd b/man/compute_historical_decompositions.PosteriorBSVARSIGN.Rd