Skip to content

Commit

Permalink
included compute_impulse_responses #25
Browse files Browse the repository at this point in the history
+ updated cpp function and...
+ moved it from bsvarTOOLs
+ R method
+ docu
+ tests
  • Loading branch information
donotdespair committed Jul 16, 2024
1 parent b4faf0e commit 5bb6248
Show file tree
Hide file tree
Showing 11 changed files with 378 additions and 172 deletions.
1 change: 1 addition & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
# Generated by roxygen2: do not edit by hand

S3method(compute_fitted_values,PosteriorBSVARSIGN)
S3method(compute_impulse_responses,PosteriorBSVARSIGN)
S3method(compute_structural_shocks,PosteriorBSVARSIGN)
S3method(estimate,BSVARSIGN)
S3method(forecast,PosteriorBSVARSIGN)
Expand Down
74 changes: 73 additions & 1 deletion R/compute.R
Original file line number Diff line number Diff line change
Expand Up @@ -118,4 +118,76 @@ compute_fitted_values.PosteriorBSVARSIGN <- function(posterior) {
class(fv) = "PosteriorFitted"

return(fv)
}
} # END compute_fitted_values.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.
#' @param standardise a logical value. If \code{TRUE}, the impulse responses are standardised
#' so that the variables' own shocks at horizon 0 are equal to 1. Otherwise, the parameter estimates
#' determine this magnitude.
#'
#' @return An object of class PosteriorIR, that is, an \code{NxNx(horizon+1)xS} array with attribute PosteriorIR
#' containing \code{S} draws of the impulse responses.
#'
#' @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 impulse responses 2 years ahead
#' irf = compute_impulse_responses(posterior, horizon = 8)
#'
#' # workflow with the pipe |>
#' ############################################################
#' set.seed(123)
#' oil |>
#' specify_bsvarSIGN$new(sign_irf = sign_irf) |>
#' estimate(S = 10) |>
#' compute_impulse_responses(horizon = 8) -> ir
#'
#'
#' @export
compute_impulse_responses.PosteriorBSVARSIGN <- function(posterior, horizon, standardise = FALSE) {

posterior_Theta0 = posterior$posterior$Theta0
posterior_A = posterior$posterior$A
posterior_A = aperm(posterior_A, c(2, 1, 3))
N = dim(posterior_A)[2]
p = posterior$last_draw$p
S = dim(posterior_A)[3]

qqq = .Call(`_bsvarSIGNs_bsvarSIGNs_ir`, posterior_A, posterior_Theta0, horizon, p, standardise)

irfs = array(NA, c(N, N, horizon + 1, S))
for (s in 1:S) irfs[,,,s] = qqq[s][[1]]
class(irfs) = "PosteriorIR"

return(irfs)
} # END compute_impulse_responses.PosteriorBSVARSIGN
78 changes: 39 additions & 39 deletions inst/include/bsvarSIGNs_RcppExports.h
Original file line number Diff line number Diff line change
Expand Up @@ -46,80 +46,80 @@ namespace bsvarSIGNs {
return Rcpp::as<arma::mat >(rcpp_result_gen);
}

inline arma::cube ir1_cpp(const arma::mat& B, const arma::mat& Theta0, int horizon, const int& p) {
typedef SEXP(*Ptr_ir1_cpp)(SEXP,SEXP,SEXP,SEXP);
static Ptr_ir1_cpp p_ir1_cpp = NULL;
if (p_ir1_cpp == NULL) {
validateSignature("arma::cube(*ir1_cpp)(const arma::mat&,const arma::mat&,int,const int&)");
p_ir1_cpp = (Ptr_ir1_cpp)R_GetCCallable("bsvarSIGNs", "_bsvarSIGNs_ir1_cpp");
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;
if (p_bsvar_sign_cpp == NULL) {
validateSignature("Rcpp::List(*bsvar_sign_cpp)(const int&,const int&,const arma::mat&,const arma::mat&,const arma::field<arma::mat>&,const arma::cube&,const arma::mat&,const arma::mat&,const arma::field<arma::mat>&,const Rcpp::List&,const Rcpp::List&,const bool,const int,const int&)");
p_bsvar_sign_cpp = (Ptr_bsvar_sign_cpp)R_GetCCallable("bsvarSIGNs", "_bsvarSIGNs_bsvar_sign_cpp");
}
RObject rcpp_result_gen;
{
RNGScope RCPP_rngScope_gen;
rcpp_result_gen = p_ir1_cpp(Shield<SEXP>(Rcpp::wrap(B)), Shield<SEXP>(Rcpp::wrap(Theta0)), Shield<SEXP>(Rcpp::wrap(horizon)), Shield<SEXP>(Rcpp::wrap(p)));
rcpp_result_gen = p_bsvar_sign_cpp(Shield<SEXP>(Rcpp::wrap(S)), Shield<SEXP>(Rcpp::wrap(p)), Shield<SEXP>(Rcpp::wrap(Y)), Shield<SEXP>(Rcpp::wrap(X)), Shield<SEXP>(Rcpp::wrap(VB)), Shield<SEXP>(Rcpp::wrap(sign_irf)), Shield<SEXP>(Rcpp::wrap(sign_narrative)), Shield<SEXP>(Rcpp::wrap(sign_B)), Shield<SEXP>(Rcpp::wrap(Z)), Shield<SEXP>(Rcpp::wrap(prior)), Shield<SEXP>(Rcpp::wrap(starting_values)), Shield<SEXP>(Rcpp::wrap(show_progress)), Shield<SEXP>(Rcpp::wrap(thin)), Shield<SEXP>(Rcpp::wrap(max_tries)));
}
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::cube >(rcpp_result_gen);
return Rcpp::as<Rcpp::List >(rcpp_result_gen);
}

inline arma::field<arma::cube> bsvarSIGNs_ir(arma::cube& posterior_B, arma::cube& posterior_Theta0, const int horizon, const int p) {
typedef SEXP(*Ptr_bsvarSIGNs_ir)(SEXP,SEXP,SEXP,SEXP);
static Ptr_bsvarSIGNs_ir p_bsvarSIGNs_ir = NULL;
if (p_bsvarSIGNs_ir == NULL) {
validateSignature("arma::field<arma::cube>(*bsvarSIGNs_ir)(arma::cube&,arma::cube&,const int,const int)");
p_bsvarSIGNs_ir = (Ptr_bsvarSIGNs_ir)R_GetCCallable("bsvarSIGNs", "_bsvarSIGNs_bsvarSIGNs_ir");
inline arma::cube bsvarSIGNs_structural_shocks(const arma::cube& posterior_B, const arma::cube& posterior_A, const arma::mat& Y, const arma::mat& X) {
typedef SEXP(*Ptr_bsvarSIGNs_structural_shocks)(SEXP,SEXP,SEXP,SEXP);
static Ptr_bsvarSIGNs_structural_shocks p_bsvarSIGNs_structural_shocks = NULL;
if (p_bsvarSIGNs_structural_shocks == NULL) {
validateSignature("arma::cube(*bsvarSIGNs_structural_shocks)(const arma::cube&,const arma::cube&,const arma::mat&,const arma::mat&)");
p_bsvarSIGNs_structural_shocks = (Ptr_bsvarSIGNs_structural_shocks)R_GetCCallable("bsvarSIGNs", "_bsvarSIGNs_bsvarSIGNs_structural_shocks");
}
RObject rcpp_result_gen;
{
RNGScope RCPP_rngScope_gen;
rcpp_result_gen = p_bsvarSIGNs_ir(Shield<SEXP>(Rcpp::wrap(posterior_B)), Shield<SEXP>(Rcpp::wrap(posterior_Theta0)), Shield<SEXP>(Rcpp::wrap(horizon)), Shield<SEXP>(Rcpp::wrap(p)));
rcpp_result_gen = p_bsvarSIGNs_structural_shocks(Shield<SEXP>(Rcpp::wrap(posterior_B)), Shield<SEXP>(Rcpp::wrap(posterior_A)), Shield<SEXP>(Rcpp::wrap(Y)), Shield<SEXP>(Rcpp::wrap(X)));
}
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);
return Rcpp::as<arma::cube >(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;
if (p_bsvar_sign_cpp == NULL) {
validateSignature("Rcpp::List(*bsvar_sign_cpp)(const int&,const int&,const arma::mat&,const arma::mat&,const arma::field<arma::mat>&,const arma::cube&,const arma::mat&,const arma::mat&,const arma::field<arma::mat>&,const Rcpp::List&,const Rcpp::List&,const bool,const int,const int&)");
p_bsvar_sign_cpp = (Ptr_bsvar_sign_cpp)R_GetCCallable("bsvarSIGNs", "_bsvarSIGNs_bsvar_sign_cpp");
inline arma::cube bsvarSIGNs_fitted_values(arma::cube& posterior_A, arma::cube& posterior_B, arma::cube& posterior_sigma, arma::mat& X) {
typedef SEXP(*Ptr_bsvarSIGNs_fitted_values)(SEXP,SEXP,SEXP,SEXP);
static Ptr_bsvarSIGNs_fitted_values p_bsvarSIGNs_fitted_values = NULL;
if (p_bsvarSIGNs_fitted_values == NULL) {
validateSignature("arma::cube(*bsvarSIGNs_fitted_values)(arma::cube&,arma::cube&,arma::cube&,arma::mat&)");
p_bsvarSIGNs_fitted_values = (Ptr_bsvarSIGNs_fitted_values)R_GetCCallable("bsvarSIGNs", "_bsvarSIGNs_bsvarSIGNs_fitted_values");
}
RObject rcpp_result_gen;
{
RNGScope RCPP_rngScope_gen;
rcpp_result_gen = p_bsvar_sign_cpp(Shield<SEXP>(Rcpp::wrap(S)), Shield<SEXP>(Rcpp::wrap(p)), Shield<SEXP>(Rcpp::wrap(Y)), Shield<SEXP>(Rcpp::wrap(X)), Shield<SEXP>(Rcpp::wrap(VB)), Shield<SEXP>(Rcpp::wrap(sign_irf)), Shield<SEXP>(Rcpp::wrap(sign_narrative)), Shield<SEXP>(Rcpp::wrap(sign_B)), Shield<SEXP>(Rcpp::wrap(Z)), Shield<SEXP>(Rcpp::wrap(prior)), Shield<SEXP>(Rcpp::wrap(starting_values)), Shield<SEXP>(Rcpp::wrap(show_progress)), Shield<SEXP>(Rcpp::wrap(thin)), Shield<SEXP>(Rcpp::wrap(max_tries)));
rcpp_result_gen = p_bsvarSIGNs_fitted_values(Shield<SEXP>(Rcpp::wrap(posterior_A)), Shield<SEXP>(Rcpp::wrap(posterior_B)), Shield<SEXP>(Rcpp::wrap(posterior_sigma)), Shield<SEXP>(Rcpp::wrap(X)));
}
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<Rcpp::List >(rcpp_result_gen);
return Rcpp::as<arma::cube >(rcpp_result_gen);
}

inline arma::cube bsvarSIGNs_structural_shocks(const arma::cube& posterior_B, const arma::cube& posterior_A, const arma::mat& Y, const arma::mat& X) {
typedef SEXP(*Ptr_bsvarSIGNs_structural_shocks)(SEXP,SEXP,SEXP,SEXP);
static Ptr_bsvarSIGNs_structural_shocks p_bsvarSIGNs_structural_shocks = NULL;
if (p_bsvarSIGNs_structural_shocks == NULL) {
validateSignature("arma::cube(*bsvarSIGNs_structural_shocks)(const arma::cube&,const arma::cube&,const arma::mat&,const arma::mat&)");
p_bsvarSIGNs_structural_shocks = (Ptr_bsvarSIGNs_structural_shocks)R_GetCCallable("bsvarSIGNs", "_bsvarSIGNs_bsvarSIGNs_structural_shocks");
inline arma::cube ir1_cpp(const arma::mat& B, const arma::mat& Theta0, int horizon, const int& p) {
typedef SEXP(*Ptr_ir1_cpp)(SEXP,SEXP,SEXP,SEXP);
static Ptr_ir1_cpp p_ir1_cpp = NULL;
if (p_ir1_cpp == NULL) {
validateSignature("arma::cube(*ir1_cpp)(const arma::mat&,const arma::mat&,int,const int&)");
p_ir1_cpp = (Ptr_ir1_cpp)R_GetCCallable("bsvarSIGNs", "_bsvarSIGNs_ir1_cpp");
}
RObject rcpp_result_gen;
{
RNGScope RCPP_rngScope_gen;
rcpp_result_gen = p_bsvarSIGNs_structural_shocks(Shield<SEXP>(Rcpp::wrap(posterior_B)), Shield<SEXP>(Rcpp::wrap(posterior_A)), Shield<SEXP>(Rcpp::wrap(Y)), Shield<SEXP>(Rcpp::wrap(X)));
rcpp_result_gen = p_ir1_cpp(Shield<SEXP>(Rcpp::wrap(B)), Shield<SEXP>(Rcpp::wrap(Theta0)), Shield<SEXP>(Rcpp::wrap(horizon)), Shield<SEXP>(Rcpp::wrap(p)));
}
if (rcpp_result_gen.inherits("interrupted-error"))
throw Rcpp::internal::InterruptedException();
Expand All @@ -130,25 +130,25 @@ namespace bsvarSIGNs {
return Rcpp::as<arma::cube >(rcpp_result_gen);
}

inline arma::cube bsvarSIGNs_fitted_values(arma::cube& posterior_A, arma::cube& posterior_B, arma::cube& posterior_sigma, arma::mat& X) {
typedef SEXP(*Ptr_bsvarSIGNs_fitted_values)(SEXP,SEXP,SEXP,SEXP);
static Ptr_bsvarSIGNs_fitted_values p_bsvarSIGNs_fitted_values = NULL;
if (p_bsvarSIGNs_fitted_values == NULL) {
validateSignature("arma::cube(*bsvarSIGNs_fitted_values)(arma::cube&,arma::cube&,arma::cube&,arma::mat&)");
p_bsvarSIGNs_fitted_values = (Ptr_bsvarSIGNs_fitted_values)R_GetCCallable("bsvarSIGNs", "_bsvarSIGNs_bsvarSIGNs_fitted_values");
inline arma::field<arma::cube> bsvarSIGNs_ir(arma::cube& posterior_B, arma::cube& posterior_Theta0, const int horizon, const int p, const bool standardise = false) {
typedef SEXP(*Ptr_bsvarSIGNs_ir)(SEXP,SEXP,SEXP,SEXP,SEXP);
static Ptr_bsvarSIGNs_ir p_bsvarSIGNs_ir = NULL;
if (p_bsvarSIGNs_ir == NULL) {
validateSignature("arma::field<arma::cube>(*bsvarSIGNs_ir)(arma::cube&,arma::cube&,const int,const int,const bool)");
p_bsvarSIGNs_ir = (Ptr_bsvarSIGNs_ir)R_GetCCallable("bsvarSIGNs", "_bsvarSIGNs_bsvarSIGNs_ir");
}
RObject rcpp_result_gen;
{
RNGScope RCPP_rngScope_gen;
rcpp_result_gen = p_bsvarSIGNs_fitted_values(Shield<SEXP>(Rcpp::wrap(posterior_A)), Shield<SEXP>(Rcpp::wrap(posterior_B)), Shield<SEXP>(Rcpp::wrap(posterior_sigma)), Shield<SEXP>(Rcpp::wrap(X)));
rcpp_result_gen = p_bsvarSIGNs_ir(Shield<SEXP>(Rcpp::wrap(posterior_B)), Shield<SEXP>(Rcpp::wrap(posterior_Theta0)), Shield<SEXP>(Rcpp::wrap(horizon)), Shield<SEXP>(Rcpp::wrap(p)), Shield<SEXP>(Rcpp::wrap(standardise)));
}
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::cube >(rcpp_result_gen);
return Rcpp::as<arma::field<arma::cube> >(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) {
Expand Down
62 changes: 62 additions & 0 deletions inst/tinytest/test_compute_impulse_responses.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,62 @@

data(oil)

set.seed(1)
suppressMessages(
specification_no1 <- specify_bsvarSIGN$new(oil)
)
run_no1 <- estimate(specification_no1, 3, 1, show_progress = FALSE)
irf <- compute_impulse_responses(run_no1, horizon = 2)

set.seed(1)
suppressMessages(
irf2 <- oil |>
specify_bsvarSIGN$new() |>
estimate(S = 3, thin = 1, show_progress = FALSE) |>
compute_impulse_responses(horizon = 2)
)


expect_error(
compute_impulse_responses(run_no1),
info = "compute_impulse_responses: specify horizon."
)

expect_identical(
irf[3,3,3,3], irf2[3,3,3,3],
info = "compute_impulse_responses: identical for normal and pipe workflow."
)


set.seed(1)
suppressMessages(
specification_no1 <- specify_bsvarSIGN$new(oil)
)
run_no1 <- estimate(specification_no1, 3, 1, show_progress = FALSE)
irf <- compute_impulse_responses(run_no1, horizon = 2, standardise = TRUE)

set.seed(1)
suppressMessages(
irf2 <- oil |>
specify_bsvarSIGN$new() |>
estimate(S = 3, thin = 1, show_progress = FALSE) |>
compute_impulse_responses(horizon = 2, standardise = TRUE)
)



expect_equal(
irf[1,1,1,1], 1,
info = "compute_impulse_responses: unit own shock at 0 horizon."
)

expect_error(
compute_impulse_responses(run_no1),
info = "compute_impulse_responses: specify horizon."
)

expect_identical(
irf[3,3,3,3], irf2[3,3,3,3],
info = "compute_impulse_responses: identical for normal and pipe workflow."
)

61 changes: 61 additions & 0 deletions man/compute_impulse_responses.PosteriorBSVARSIGN.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

Loading

0 comments on commit 5bb6248

Please sign in to comment.