diff --git a/R/estimate.R b/R/estimate.R index 8cca0fa..e58b3f9 100644 --- a/R/estimate.R +++ b/R/estimate.R @@ -109,10 +109,12 @@ estimate.BSVARSIGN = function(specification, S, thin = 1, show_progress = TRUE) Z = get_Z(identification$sign_irf) sign = identification$sign_irf sign[is.na(sign)] = 0 + narrative = identification$sign_narrative + narrative[, 5] = narrative[, 5] - p # estimation qqq = .Call(`_bsvarSIGNs_bsvar_sign_cpp`, S, p, Y, X, - identification$VB, sign, identification$sign_narrative, + identification$VB, sign, narrative, identification$sign_relation, Z, prior, starting_values, show_progress, thin, max_tries) diff --git a/R/specify.R b/R/specify.R index 2652dd1..cd592b9 100644 --- a/R/specify.R +++ b/R/specify.R @@ -468,7 +468,7 @@ specify_identification_bsvarSIGN = R6::R6Class( missing_all = FALSE } if (missing(sign_narrative)) { - sign_narrative = matrix(c(0, 1, 1, 1, 1, 0), ncol = 6, nrow = 1) + sign_narrative = matrix(c(0, 1, 1, 1, 1, 1), ncol = 6, nrow = 1) } else { missing_all = FALSE } @@ -558,7 +558,7 @@ specify_identification_bsvarSIGN = R6::R6Class( missing_all = FALSE } if (missing(sign_narrative)) { - sign_narrative = matrix(c(0, 1, 1, 1, 1, 0), ncol = 6, nrow = 1) + sign_narrative = matrix(c(0, 1, 1, 1, 1, 1), ncol = 6, nrow = 1) } else { missing_all = FALSE } @@ -683,7 +683,7 @@ specify_bsvarSIGN = R6::R6Class( missing_all = FALSE } if (missing(sign_narrative)) { - sign_narrative = matrix(c(0, 1, 1, 1, 1, 0), ncol = 6, nrow = 1) + sign_narrative = matrix(c(0, 1, 1, 1, 1, 1), ncol = 6, nrow = 1) } else { missing_all = FALSE } diff --git a/README.Rmd b/README.Rmd index fdd3f6d..7c0df8c 100644 --- a/README.Rmd +++ b/README.Rmd @@ -71,15 +71,16 @@ plot(irf, probability = 0.68) A replication of [Antolín-Díaz and Rubio-Ramírez (2018)](http://doi.org/10.1257/aer.20161852). + ```r 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, 6)) +sign_irf = array(sign_irf, dim = c(6, 6, 5)) -# in October 1979 +# in October 1979 the shock sign_narrative = rbind(c(1, 1, NA, 1, 166, 0), # is positive c(3, 1, 6, 1, 166, 0)) # greatest historical decomposition @@ -93,3 +94,26 @@ 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 f184045..75c620d 100644 --- a/README.md +++ b/README.md @@ -98,9 +98,9 @@ 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, 6)) +sign_irf = array(sign_irf, dim = c(6, 6, 5)) -# in October 1979 +# in October 1979 the shock sign_narrative = rbind(c(1, 1, NA, 1, 166, 0), # is positive c(3, 1, 6, 1, 166, 0)) # greatest historical decomposition @@ -113,3 +113,26 @@ 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) +``` diff --git a/src/bsvars_sign.cpp b/src/bsvars_sign.cpp index 1433621..33470ba 100644 --- a/src/bsvars_sign.cpp +++ b/src/bsvars_sign.cpp @@ -13,14 +13,6 @@ 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( @@ -55,7 +47,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; @@ -104,36 +96,32 @@ 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++) { - 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; - + 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); - checkUserInterrupt(); + // 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)); // sample reduced-form parameters Sigma = iwishrnd(post_S, post_nu); @@ -157,6 +145,9 @@ Rcpp::List bsvar_sign_cpp( posterior_Theta0.slice(s) = chol_Sigma * Q; posterior_shocks.slice(s) = shocks; + // Check for user interrupts + if (s % 1 == 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 3894f83..ea4a57a 100644 --- a/src/restrictions_narrative.cpp +++ b/src/restrictions_narrative.cpp @@ -81,7 +81,7 @@ double weight_narrative( const int M = 1e+3; // number of draws to approximate normal distribution - double n_success = 1.0e-15; + double n_success = 1.0e-6; cube Z(irf.n_rows, sign_narrative.col(5).max() + 1, M, fill::randn); // cube Z(irf.n_rows, T, M, fill::randn);