Skip to content

Commit

Permalink
Merge pull request #29 from bsvars/27-coherent-bsvars-notation
Browse files Browse the repository at this point in the history 
27 coherent bsvars notation
  • Loading branch information
@adamwang15
adamwang15 authored Jul 15, 2024
2 parents 93075c0 + ce7db40 commit 2787246
Show file tree
Hide file tree
Showing 10 changed files with 87 additions and 75 deletions.
13 changes: 10 additions & 3 deletions R/estimate.BSVARSIGN.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,16 +96,23 @@ estimate.BSVARSIGN = function(specification, S, thin = 1, show_progress = TRUE)

# get the inputs to estimation
# prior = specification$last_draw$prior$get_prior()
prior = specification$prior
prior = specification$prior$get_prior()
starting_values = specification$starting_values$get_starting_values()
identification = specification$identification$get_identification()
max_tries = identification$max_tries
data_matrices = specification$data_matrices$get_data_matrices()
p = specification$p

prior$B = t(prior$A)
prior$Ysoc = t(prior$Ysoc)
prior$Xsoc = t(prior$Xsoc)
prior$Ysur = t(prior$Ysur)
prior$Xsur = t(prior$Xsur)
Y = t(data_matrices$Y)
X = t(data_matrices$X)

# estimation
qqq = .Call(`_bsvarSIGNs_bsvar_sign_cpp`, S, p,
data_matrices$Y, data_matrices$X, identification$VB,
qqq = .Call(`_bsvarSIGNs_bsvar_sign_cpp`, S, p, Y, X, identification$VB,
identification$sign_irf, identification$sign_narrative,
identification$sign_B, identification$zero_irf,
prior, starting_values, show_progress, thin, max_tries)
Expand Down
66 changes: 32 additions & 34 deletions R/specify_bsvarSIGN.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,23 +139,20 @@ specify_prior_bsvarSIGN = R6::R6Class(
"PriorBSVARSIGN",

public = list(
#' @field p a positive integer - the number of lags.
p = 1,

#' @field hyper a \code{(N+3)xS} matrix of hyper-parameters \eqn{\mu, \delta, \lambda, \psi}.
hyper = matrix(),

#' @field B a \code{KxN} normal prior mean matrix for the autoregressive
#' @field A a \code{NxK} normal prior mean matrix for the autoregressive
#' parameters.
B = matrix(),
A = matrix(),

#' @field V a \code{KxK} matrix determining the normal prior column-specific
#' covariance for the autoregressive parameters.
V = matrix(),

#' @field Vp a \code{px1} vector of lag shrinkage level.
Vp = matrix(),

#' @field Vd a \code{(d+1)x1} vector of (large) variances for constants.
Vd = matrix(),

#' @field S an \code{NxN} matrix determining the inverted-Wishart prior scale
#' of error terms covariance matrix.
S = matrix(),
Expand All @@ -167,22 +164,22 @@ specify_prior_bsvarSIGN = R6::R6Class(
#' @field data an \code{TxN} matrix of observations.
data = matrix(),

#' @field Y an \code{TxN} matrix of dependent variables.
#' @field Y an \code{NxT} matrix of dependent variables.
Y = matrix(),

#' @field X an \code{TxK} matrix of independent variables.
#' @field X an \code{KxT} matrix of independent variables.
X = matrix(),

#' @field Ysoc an \code{NxN} matrix with the sum-of-coefficients dummy observations.
Ysoc = matrix(),

#' @field Xsoc an \code{NxK} matrix with the sum-of-coefficients dummy observations.
#' @field Xsoc an \code{KxN} matrix with the sum-of-coefficients dummy observations.
Xsoc = matrix(),

#' @field Ysur an \code{NxN} matrix with the single-unit-root dummy observations.
Ysur = matrix(),

#' @field Xsur an \code{NxK} matrix with the single-unit-root dummy observations.
#' @field Xsur an \code{KxN} matrix with the single-unit-root dummy observations.
Xsur = matrix(),

#' @field mu.scale a positive scalar - the shape of the gamma prior for \eqn{\mu}.
Expand Down Expand Up @@ -224,7 +221,7 @@ specify_prior_bsvarSIGN = R6::R6Class(
#' prior = specify_prior_bsvarSIGN$new(oil, p = 1)
#' prior$B # show autoregressive prior mean
#'
initialize = function(data, p, exogenous = NULL, stationary = rep(FALSE, dim(data)[2])) {
initialize = function(data, p, exogenous = NULL, stationary = rep(FALSE, ncol(data))) {

stopifnot("Argument p must be a positive integer number." = p > 0 & p %% 1 == 0)

Expand All @@ -245,9 +242,7 @@ specify_prior_bsvarSIGN = R6::R6Class(
B = matrix(0, K, N)
B[1:N,] = diag(!stationary)

Vp = (1:p)^-2
Vd = rep(100, 1 + d)
V = diag(c(kronecker(Vp, rep(1, N)), Vd))
V = diag(c(kronecker((1:p)^-2, rep(1, N)), rep(100, 1 + d)))

s2.ols = rep(NA, N)
for (n in 1:N) {
Expand Down Expand Up @@ -282,20 +277,18 @@ specify_prior_bsvarSIGN = R6::R6Class(
# }
# Xstar = cbind(Xstar, c(rep(0, N), 1), matrix(0, N + 1, d))

self$Y = Y
self$X = X
self$Vp = Vp
self$Vd = Vd

self$p = p
self$hyper = hyper
self$B = B
self$A = t(B)
self$V = V
self$S = diag(N)
self$nu = N + 2
self$Ysoc = Ysoc
self$Xsoc = Xsoc
self$Ysur = Ysur
self$Xsur = Xsur
self$Y = t(Y)
self$X = t(X)
self$Ysoc = t(Ysoc)
self$Xsoc = t(Xsoc)
self$Ysur = t(Ysur)
self$Xsur = t(Xsur)
self$mu.scale = scale
self$mu.shape = shape
self$delta.scale = scale
Expand All @@ -316,11 +309,10 @@ specify_prior_bsvarSIGN = R6::R6Class(
#'
get_prior = function(){
list(
p = self$p,
hyper = self$hyper,
B = self$B,
A = self$A,
V = self$V,
Vp = self$Vp,
Vd = self$Vd,
S = self$S,
nu = self$nu,
Ysoc = self$Ysoc,
Expand Down Expand Up @@ -373,10 +365,16 @@ specify_prior_bsvarSIGN = R6::R6Class(
init = narrow_hyper(model, hyper)
prior = self$get_prior()

prior$B = t(prior$A)
prior$Ysoc = t(prior$Ysoc)
prior$Xsoc = t(prior$Xsoc)
prior$Ysur = t(prior$Ysur)
prior$Xsur = t(prior$Xsur)

result = stats::optim(
init,
\(x) -log_posterior_hyper(extend_hyper(hyper, model, matrix(x)),
model, self$Y, self$X, prior),
model, t(self$Y), t(self$X), prior),
method = 'L-BFGS-B',
lower = rep(0, length(init)),
upper = init * 100,
Expand All @@ -393,7 +391,8 @@ specify_prior_bsvarSIGN = R6::R6Class(
variance = e$vectors %*% diag(as.vector(1 / abs(e$values))) %*% t(e$vectors)
}

self$hyper = sample_hyper(S, burn_in, mode, model, self$Y, self$X, variance, prior)
self$hyper = sample_hyper(S, burn_in, mode, model,
t(self$Y), t(self$X), variance, prior)
self$hyper = self$hyper[, -(1:burn_in)]
} # END estimate_hyper

Expand Down Expand Up @@ -716,15 +715,14 @@ specify_bsvarSIGN = R6::R6Class(
B[lower.tri(B, diag = TRUE)] = TRUE

self$data_matrices = bsvars::specify_data_matrices$new(data, p, exogenous)
self$data_matrices$Y = t(self$data_matrices$Y)
self$data_matrices$X = t(self$data_matrices$X)
self$identification = specify_identification_bsvarSIGN$new(N,
sign_irf,
sign_narrative,
sign_B,
zero_irf,
max_tries)
self$prior = specify_prior_bsvarSIGN$new(data, p, exogenous, stationary)
self$prior = specify_prior_bsvarSIGN$new(data, p, exogenous,
stationary)
self$starting_values = bsvars::specify_starting_values_bsvar$new(N, self$p, d)
}, # END initialize

Expand Down
Binary file modified data/oil.rda
View file
Binary file not shown.
4 changes: 2 additions & 2 deletions inst/include/bsvarSIGNs_RcppExports.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,7 @@ namespace bsvarSIGNs {
return Rcpp::as<arma::field<arma::cube> >(rcpp_result_gen);
}

inline Rcpp::List bsvar_sign_cpp(const int& S, const int& lags, 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) {
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) {
Expand All @@ -98,7 +98,7 @@ namespace bsvarSIGNs {
RObject rcpp_result_gen;
{
RNGScope RCPP_rngScope_gen;
rcpp_result_gen = p_bsvar_sign_cpp(Shield<SEXP>(Rcpp::wrap(S)), Shield<SEXP>(Rcpp::wrap(lags)), 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_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();
Expand Down
12 changes: 6 additions & 6 deletions inst/tinytest/test_specify.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,13 +8,13 @@ spec = specify_bsvarSIGN$new(oil)
expect_identical(class(spec)[1],
"BSVARSIGN")

expect_identical(dim(spec$data_matrices$Y)[1],
dim(spec$data_matrices$X)[1])
expect_identical(dim(spec$data_matrices$Y)[2],
dim(spec$data_matrices$X)[2])

expect_identical(length(spec$data_matrices$get_data_matrices()),
2L)

expect_identical(dim(spec$prior$B)[1], dim(spec$data_matrices$X)[2])
expect_identical(dim(spec$prior$A)[2], dim(spec$data_matrices$X)[1])

expect_identical(class(spec$prior$get_prior()),
"list")
Expand Down Expand Up @@ -51,13 +51,13 @@ spec = specify_bsvarSIGN$new(oil,
expect_identical(class(spec)[1],
"BSVARSIGN")

expect_identical(dim(spec$data_matrices$Y)[1],
dim(spec$data_matrices$X)[1])
expect_identical(dim(spec$data_matrices$Y)[2],
dim(spec$data_matrices$X)[2])

expect_identical(length(spec$data_matrices$get_data_matrices()),
2L)

expect_identical(dim(spec$prior$B)[1], dim(spec$data_matrices$X)[2])
expect_identical(dim(spec$prior$A)[2], dim(spec$data_matrices$X)[1])

expect_identical(class(spec$prior$get_prior()),
"list")
Expand Down
9 changes: 5 additions & 4 deletions inst/varia/nicetry.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,17 +4,18 @@ name = sapply(1:7, \(i) data$ShortDescr[[i]][[1]])
data = data$y
colnames(data) = name

data = optimism
# data = optimism

spec = specify_bsvarSIGN$new(data, p = 4)

start_time = Sys.time()

spec$prior$estimate_hyper(S = 10000, burn_in = 5000,
mu = TRUE, delta = TRUE, lambda = TRUE, psi = TRUE)
# post = estimate(spec, S = 1000)
# irf = compute_impulse_responses(post, horizon = 40)
# plot(irf, probability = 0.68)
post = estimate(spec, S = 100)
class(post) = "PosteriorBSVAR"
irf = compute_impulse_responses(post, horizon = 40)
plot(irf, probability = 0.68)

end_time = Sys.time()
end_time - start_time
Expand Down
18 changes: 8 additions & 10 deletions man/specify_prior_bsvarSIGN.Rd
View file

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

12 changes: 6 additions & 6 deletions src/RcppExports.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,12 +127,12 @@ RcppExport SEXP _bsvarSIGNs_bsvarSIGNs_ir(SEXP posterior_BSEXP, SEXP posterior_T
return rcpp_result_gen;
}
// bsvar_sign_cpp
Rcpp::List bsvar_sign_cpp(const int& S, const int& lags, 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, const int thin, const int& max_tries);
static SEXP _bsvarSIGNs_bsvar_sign_cpp_try(SEXP SSEXP, SEXP lagsSEXP, SEXP YSEXP, SEXP XSEXP, SEXP VBSEXP, SEXP sign_irfSEXP, SEXP sign_narrativeSEXP, SEXP sign_BSEXP, SEXP ZSEXP, SEXP priorSEXP, SEXP starting_valuesSEXP, SEXP show_progressSEXP, SEXP thinSEXP, SEXP max_triesSEXP) {
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, const int thin, const int& max_tries);
static SEXP _bsvarSIGNs_bsvar_sign_cpp_try(SEXP SSEXP, SEXP pSEXP, SEXP YSEXP, SEXP XSEXP, SEXP VBSEXP, SEXP sign_irfSEXP, SEXP sign_narrativeSEXP, SEXP sign_BSEXP, SEXP ZSEXP, SEXP priorSEXP, SEXP starting_valuesSEXP, SEXP show_progressSEXP, SEXP thinSEXP, SEXP max_triesSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::traits::input_parameter< const int& >::type S(SSEXP);
Rcpp::traits::input_parameter< const int& >::type lags(lagsSEXP);
Rcpp::traits::input_parameter< const int& >::type p(pSEXP);
Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP);
Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP);
Rcpp::traits::input_parameter< const arma::field<arma::mat>& >::type VB(VBSEXP);
Expand All @@ -145,15 +145,15 @@ BEGIN_RCPP
Rcpp::traits::input_parameter< const bool >::type show_progress(show_progressSEXP);
Rcpp::traits::input_parameter< const int >::type thin(thinSEXP);
Rcpp::traits::input_parameter< const int& >::type max_tries(max_triesSEXP);
rcpp_result_gen = Rcpp::wrap(bsvar_sign_cpp(S, lags, Y, X, VB, sign_irf, sign_narrative, sign_B, Z, prior, starting_values, show_progress, thin, max_tries));
rcpp_result_gen = Rcpp::wrap(bsvar_sign_cpp(S, p, Y, X, VB, sign_irf, sign_narrative, sign_B, Z, prior, starting_values, show_progress, thin, max_tries));
return rcpp_result_gen;
END_RCPP_RETURN_ERROR
}
RcppExport SEXP _bsvarSIGNs_bsvar_sign_cpp(SEXP SSEXP, SEXP lagsSEXP, SEXP YSEXP, SEXP XSEXP, SEXP VBSEXP, SEXP sign_irfSEXP, SEXP sign_narrativeSEXP, SEXP sign_BSEXP, SEXP ZSEXP, SEXP priorSEXP, SEXP starting_valuesSEXP, SEXP show_progressSEXP, SEXP thinSEXP, SEXP max_triesSEXP) {
RcppExport SEXP _bsvarSIGNs_bsvar_sign_cpp(SEXP SSEXP, SEXP pSEXP, SEXP YSEXP, SEXP XSEXP, SEXP VBSEXP, SEXP sign_irfSEXP, SEXP sign_narrativeSEXP, SEXP sign_BSEXP, SEXP ZSEXP, SEXP priorSEXP, SEXP starting_valuesSEXP, SEXP show_progressSEXP, SEXP thinSEXP, SEXP max_triesSEXP) {
SEXP rcpp_result_gen;
{
Rcpp::RNGScope rcpp_rngScope_gen;
rcpp_result_gen = PROTECT(_bsvarSIGNs_bsvar_sign_cpp_try(SSEXP, lagsSEXP, YSEXP, XSEXP, VBSEXP, sign_irfSEXP, sign_narrativeSEXP, sign_BSEXP, ZSEXP, priorSEXP, starting_valuesSEXP, show_progressSEXP, thinSEXP, max_triesSEXP));
rcpp_result_gen = PROTECT(_bsvarSIGNs_bsvar_sign_cpp_try(SSEXP, pSEXP, YSEXP, XSEXP, VBSEXP, sign_irfSEXP, sign_narrativeSEXP, sign_BSEXP, ZSEXP, priorSEXP, starting_valuesSEXP, show_progressSEXP, thinSEXP, max_triesSEXP));
}
Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error");
if (rcpp_isInterrupt_gen) {
Expand Down
Loading

0 comments on commit 2787246

Please sign in to comment.