From 520dc878976701c52f2e7360a5e7f4a088fb1dea Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Tue, 16 Jul 2024 23:45:56 +1000 Subject: [PATCH 1/3] do not resample hyper when Q is invalid #31 --- src/bsvars_sign.cpp | 55 +++++++++++++++++----------------- src/restrictions_narrative.cpp | 2 +- 2 files changed, 29 insertions(+), 28 deletions(-) diff --git a/src/bsvars_sign.cpp b/src/bsvars_sign.cpp index c43cbde..1d23d4f 100644 --- a/src/bsvars_sign.cpp +++ b/src/bsvars_sign.cpp @@ -96,32 +96,36 @@ Rcpp::List bsvar_sign_cpp( // #pragma omp parallel for private(hyper, mu, delta, lambda, psi, prior_V, prior_S, Ystar, Xstar, Yplus, Xplus, result, post_B, post_V, post_S, Sigma, chol_Sigma, B, h_invp, Q, shocks, w) for (int s = 0; s < S; s++) { - w = 0; + hyper = hypers.col(randi(distr_param(0, S_hyper))); + mu = hyper(0); + delta = hyper(1); + lambda = hyper(2); + psi = hyper.rows(3, N + 2); + + // update Minnesota prior + prior_V = diagmat(prior_v % + join_vert(lambda * lambda * repmat(1 / psi, p, 1), + ones(K - N * p))); + prior_S = diagmat(psi); + + // update dummy observation prior + Ystar = join_vert(Ysoc / mu, Ysur / delta); + Xstar = join_vert(Xsoc / mu, Xsur / delta); + Yplus = join_vert(Ystar, Y); + Xplus = join_vert(Xstar, X); + + // posterior parameters + result = niw_cpp(Yplus, Xplus, prior_B, prior_V, prior_S, prior_nu); + post_B = result(0); + post_V = result(1); + post_S = result(2); + post_nu = as_scalar(result(3)); + + w = 0; + while (w == 0) { - hyper = hypers.col(randi(distr_param(0, S_hyper))); - mu = hyper(0); - delta = hyper(1); - lambda = hyper(2); - psi = hyper.rows(3, N + 2); - - // update Minnesota prior - prior_V = diagmat(prior_v % - join_vert(lambda * lambda * repmat(1 / psi, p, 1), - ones(K - N * p))); - prior_S = diagmat(psi); - // update dummy observation prior - Ystar = join_vert(Ysoc / mu, Ysur / delta); - Xstar = join_vert(Xsoc / mu, Xsur / delta); - Yplus = join_vert(Ystar, Y); - Xplus = join_vert(Xstar, X); - - // posterior parameters - result = niw_cpp(Yplus, Xplus, prior_B, prior_V, prior_S, prior_nu); - post_B = result(0); - post_V = result(1); - post_S = result(2); - post_nu = as_scalar(result(3)); + checkUserInterrupt(); // sample reduced-form parameters Sigma = iwishrnd(post_S, post_nu); @@ -145,9 +149,6 @@ Rcpp::List bsvar_sign_cpp( posterior_Theta0.slice(s) = chol_Sigma * Q; posterior_shocks.slice(s) = shocks; - // Check for user interrupts - if (s % 200 == 0) checkUserInterrupt(); - // Increment progress bar if (any(prog_rep_points == s)) bar.increment(); diff --git a/src/restrictions_narrative.cpp b/src/restrictions_narrative.cpp index 805bb66..3894f83 100644 --- a/src/restrictions_narrative.cpp +++ b/src/restrictions_narrative.cpp @@ -79,7 +79,7 @@ double weight_narrative( const arma::cube& irf ) { - const int M = 1e+04; // number of draws to approximate normal distribution + const int M = 1e+3; // number of draws to approximate normal distribution double n_success = 1.0e-15; From b26a2652255c2f60a0265ed266f489926714eadf Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Wed, 17 Jul 2024 00:00:59 +1000 Subject: [PATCH 2/3] update readme example --- README.Rmd | 49 ++++++++++++++++++++++--------------------------- README.md | 49 ++++++++++++++++++++++--------------------------- 2 files changed, 44 insertions(+), 54 deletions(-) diff --git a/README.Rmd b/README.Rmd index 85c0199..6054dbb 100644 --- a/README.Rmd +++ b/README.Rmd @@ -29,17 +29,35 @@ devtools::install_github("bsvars/bsvarSIGNs") # Example -A repliaction of Antolín-Díaz and Rubio-Ramírez (2018). +A replication of Arias, Rubio-Ramírez and Waggoner (2018). + +```r +data(optimism) + +# optimism shock +# no effect on productivity + positive effect on stock prices +sign_irf = matrix(c(0, 1, rep(NA, 23)), 5, 5) + +specification = specify_bsvarSIGN$new(optimism * 100, + p = 4, + sign_irf = sign_irf) + +posterior = estimate(specification, S = 100) +irf = compute_impulse_responses(posterior, horizon = 40) +plot(irf, probability = 0.68) +``` + +A replication of Antolín-Díaz and Rubio-Ramírez (2018). ```r data(monetary) # contractionary monetary policy shock -sign_irf = matrix(0, 6, 6) -sign_irf[, 1] = c(0, -1, -1, 0, -1, 1) +sign_irf = matrix(NA, 6, 6) +sign_irf[, 1] = c(NA, -1, -1, NA, -1, 1) sign_irf = array(sign_irf, dim = c(6, 6, 5)) -# in October 1979 the shock +# in October 1979 sign_narrative = rbind(c(1, 1, NA, 1, 166, 0), # is positive c(3, 1, 6, 1, 166, 0)) # greatest historical decomposition @@ -53,26 +71,3 @@ irf = compute_impulse_responses(posterior, horizon = 60) plot(irf, probability = 0.68) ``` -A repliaction of Arias, Rubio-Ramírez and Waggoner (2018). - -```r -data(optimism) - -# optimism shock -# no contemporaneous effect on productivity -zero_irf = matrix(0, nrow = 5, ncol = 5) -zero_irf[1, 1] = 1 -# positive contemporaneous effect on stock prices -sign_irf = array(0, dim = c(5, 5, 1)) -sign_irf[2, 1, 1] = 1 - -specification = specify_bsvarSIGN$new(optimism * 100, - p = 4, - sign_irf = sign_irf, - zero_irf = zero_irf) - -posterior = estimate(specification, S = 100) -irf = compute_impulse_responses(posterior, horizon = 40) -plot(irf, probability = 0.68) -``` - diff --git a/README.md b/README.md index 07252d4..4d23d3e 100644 --- a/README.md +++ b/README.md @@ -18,17 +18,35 @@ devtools::install_github("bsvars/bsvarSIGNs") # Example -A repliaction of Antolín-Díaz and Rubio-Ramírez (2018). +A replication of Arias, Rubio-Ramírez and Waggoner (2018). + +``` r +data(optimism) + +# optimism shock +# no effect on productivity + positive effect on stock prices +sign_irf = matrix(c(0, 1, rep(NA, 23)), 5, 5) + +specification = specify_bsvarSIGN$new(optimism * 100, + p = 4, + sign_irf = sign_irf) + +posterior = estimate(specification, S = 100) +irf = compute_impulse_responses(posterior, horizon = 40) +plot(irf, probability = 0.68) +``` + +A replication of Antolín-Díaz and Rubio-Ramírez (2018). ``` r data(monetary) # contractionary monetary policy shock -sign_irf = matrix(0, 6, 6) -sign_irf[, 1] = c(0, -1, -1, 0, -1, 1) +sign_irf = matrix(NA, 6, 6) +sign_irf[, 1] = c(NA, -1, -1, NA, -1, 1) sign_irf = array(sign_irf, dim = c(6, 6, 5)) -# in October 1979 the shock +# in October 1979 sign_narrative = rbind(c(1, 1, NA, 1, 166, 0), # is positive c(3, 1, 6, 1, 166, 0)) # greatest historical decomposition @@ -41,26 +59,3 @@ posterior = estimate(specification, S = 100) irf = compute_impulse_responses(posterior, horizon = 60) plot(irf, probability = 0.68) ``` - -A repliaction of Arias, Rubio-Ramírez and Waggoner (2018). - -``` r -data(optimism) - -# optimism shock -# no contemporaneous effect on productivity -zero_irf = matrix(0, nrow = 5, ncol = 5) -zero_irf[1, 1] = 1 -# positive contemporaneous effect on stock prices -sign_irf = array(0, dim = c(5, 5, 1)) -sign_irf[2, 1, 1] = 1 - -specification = specify_bsvarSIGN$new(optimism * 100, - p = 4, - sign_irf = sign_irf, - zero_irf = zero_irf) - -posterior = estimate(specification, S = 100) -irf = compute_impulse_responses(posterior, horizon = 40) -plot(irf, probability = 0.68) -``` From b14010078934886adad7699f53dafd5f403bd85c Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Wed, 17 Jul 2024 00:29:29 +1000 Subject: [PATCH 3/3] comment about notations --- README.Rmd | 2 +- README.md | 2 +- src/bsvars_sign.cpp | 10 +++++++++- 3 files changed, 11 insertions(+), 3 deletions(-) diff --git a/README.Rmd b/README.Rmd index 6054dbb..c18cfda 100644 --- a/README.Rmd +++ b/README.Rmd @@ -55,7 +55,7 @@ data(monetary) # contractionary monetary policy shock sign_irf = matrix(NA, 6, 6) sign_irf[, 1] = c(NA, -1, -1, NA, -1, 1) -sign_irf = array(sign_irf, dim = c(6, 6, 5)) +sign_irf = array(sign_irf, dim = c(6, 6, 6)) # in October 1979 sign_narrative = rbind(c(1, 1, NA, 1, 166, 0), # is positive diff --git a/README.md b/README.md index 4d23d3e..00ed0ed 100644 --- a/README.md +++ b/README.md @@ -44,7 +44,7 @@ data(monetary) # contractionary monetary policy shock sign_irf = matrix(NA, 6, 6) sign_irf[, 1] = c(NA, -1, -1, NA, -1, 1) -sign_irf = array(sign_irf, dim = c(6, 6, 5)) +sign_irf = array(sign_irf, dim = c(6, 6, 6)) # in October 1979 sign_narrative = rbind(c(1, 1, NA, 1, 166, 0), # is positive diff --git a/src/bsvars_sign.cpp b/src/bsvars_sign.cpp index 1d23d4f..1433621 100644 --- a/src/bsvars_sign.cpp +++ b/src/bsvars_sign.cpp @@ -13,6 +13,14 @@ using namespace Rcpp; using namespace arma; +/*** +All notations in the C++ code except for compute.cpp and forecast_bsvarSIGNs.cpp + are consistent with the notations in the papers: + Antolín-Díaz and Rubio-Ramírez (2018) and Arias, Rubio-Ramírez and Waggoner (2018) + which are different from the notations in the R code. +***/ + + // [[Rcpp::interfaces(cpp)]] // [[Rcpp::export]] Rcpp::List bsvar_sign_cpp( @@ -47,7 +55,7 @@ Rcpp::List bsvar_sign_cpp( if (show_progress) { Rcout << "**************************************************|" << endl; Rcout << " bsvarSIGNs: Bayesian Structural VAR with zero, |" << endl; - Rcout << " sign and narrative restrictions |" << endl; + Rcout << " sign and narrative restrictions |" << endl; Rcout << "**************************************************|" << endl; // Rcout << " Gibbs sampler for the SVAR model |" << endl; // Rcout << "**************************************************|" << endl;