From f690c40eb53d6168b15fb3d80de459d98b705af2 Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Tue, 11 Jun 2024 13:51:08 +1000 Subject: [PATCH 1/4] mcmc scaling constant #13 the scaling constant parameter c in metropolis-hastings --- R/RcppExports.R | 4 ++-- src/RcppExports.cpp | 9 +++++---- src/sample_hyper.cpp | 25 +++++++++++++++---------- src/sample_hyper.h | 1 + 4 files changed, 23 insertions(+), 16 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index b08a07f..2ca23a1 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -41,8 +41,8 @@ log_posterior_hyper <- function(p, hyper, Y, X, prior) { .Call(`_bsvarSIGNs_log_posterior_hyper`, p, hyper, Y, X, prior) } -sample_hyper <- function(S, p, Y, X, prior) { - .Call(`_bsvarSIGNs_sample_hyper`, S, p, Y, X, prior) +sample_hyper <- function(S, p, c, Y, X, prior) { + .Call(`_bsvarSIGNs_sample_hyper`, S, p, c, Y, X, prior) } # Register entry points for exported C++ functions diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index fd3e3d2..ea3e5a8 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -710,17 +710,18 @@ BEGIN_RCPP END_RCPP } // sample_hyper -arma::mat sample_hyper(const int& S, const int& p, const arma::mat& Y, const arma::mat& X, const Rcpp::List& prior); -RcppExport SEXP _bsvarSIGNs_sample_hyper(SEXP SSEXP, SEXP pSEXP, SEXP YSEXP, SEXP XSEXP, SEXP priorSEXP) { +arma::mat sample_hyper(const int& S, const int& p, const double& c, const arma::mat& Y, const arma::mat& X, const Rcpp::List& prior); +RcppExport SEXP _bsvarSIGNs_sample_hyper(SEXP SSEXP, SEXP pSEXP, SEXP cSEXP, SEXP YSEXP, SEXP XSEXP, SEXP priorSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< const int& >::type S(SSEXP); Rcpp::traits::input_parameter< const int& >::type p(pSEXP); + Rcpp::traits::input_parameter< const double& >::type c(cSEXP); 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 Rcpp::List& >::type prior(priorSEXP); - rcpp_result_gen = Rcpp::wrap(sample_hyper(S, p, Y, X, prior)); + rcpp_result_gen = Rcpp::wrap(sample_hyper(S, p, c, Y, X, prior)); return rcpp_result_gen; END_RCPP } @@ -904,7 +905,7 @@ static const R_CallMethodDef CallEntries[] = { {"_bsvarSIGNs_log_prior_hyper", (DL_FUNC) &_bsvarSIGNs_log_prior_hyper, 2}, {"_bsvarSIGNs_log_ml_dummy", (DL_FUNC) &_bsvarSIGNs_log_ml_dummy, 5}, {"_bsvarSIGNs_log_posterior_hyper", (DL_FUNC) &_bsvarSIGNs_log_posterior_hyper, 5}, - {"_bsvarSIGNs_sample_hyper", (DL_FUNC) &_bsvarSIGNs_sample_hyper, 5}, + {"_bsvarSIGNs_sample_hyper", (DL_FUNC) &_bsvarSIGNs_sample_hyper, 6}, {"_bsvarSIGNs_qr_sign_cpp", (DL_FUNC) &_bsvarSIGNs_qr_sign_cpp, 1}, {"_bsvarSIGNs_rortho_cpp", (DL_FUNC) &_bsvarSIGNs_rortho_cpp, 1}, {"_bsvarSIGNs_match_sign", (DL_FUNC) &_bsvarSIGNs_match_sign, 2}, diff --git a/src/sample_hyper.cpp b/src/sample_hyper.cpp index fbdd587..128f4a9 100644 --- a/src/sample_hyper.cpp +++ b/src/sample_hyper.cpp @@ -20,8 +20,16 @@ Rcpp::List update_prior( vec phi = hyper.rows(3, n_hyper - 1); mat prior_B = as(prior["B"]); - mat prior_V = lambda * lambda * as(prior["V"]); - mat prior_S = diagmat(phi); + + int K = prior_B.n_rows; + vec v(K); + v(K - 1) = 1e+6; + v.rows(0, K - 2) = lambda * lambda * kron(pow(linspace(1, p, p), -2), 1 / phi); + mat prior_V = lambda * lambda * as(prior["V"]); + // mat prior_V = diagmat(v); + + mat prior_S = as(prior["S"]); + // mat prior_S = diagmat(phi); int prior_nu = as(prior["nu"]); return List::create( @@ -196,6 +204,9 @@ double log_ml_dummy( double log_ml_dummy = log_ml(p, prior_B, prior_V, prior_S, prior_nu, chol_V, cholinv_S, inv_V, Ystar, Xstar); + return log_ml(p, prior_B, prior_V, prior_S, prior_nu, + chol_V, cholinv_S, inv_V, Y, X); + return log_ml_extended - log_ml_dummy; } @@ -224,6 +235,7 @@ double log_posterior_hyper( arma::mat sample_hyper( const int& S, const int& p, + const double& c, const arma::mat& Y, const arma::mat& X, const Rcpp::List& prior @@ -236,8 +248,7 @@ arma::mat sample_hyper( double logp = log_posterior_hyper(p, hyper, Y, X, prior); - double c = 1; - mat W = as(prior["W"]); + mat W = c * as(prior["W"]); mat posterior_hyper(n_hyper, S); @@ -248,12 +259,6 @@ arma::mat sample_hyper( double logpp = log_posterior_hyper(p, hyperp, Y, X, prior); double r = exp(logpp - logp); - if (r < 0.2) { - c /= 10; - } else if (r > 0.5) { - c *= 2; - } - if (randu() < r) { hyper = hyperp; logp = logpp; diff --git a/src/sample_hyper.h b/src/sample_hyper.h index 1fa2382..9648922 100644 --- a/src/sample_hyper.h +++ b/src/sample_hyper.h @@ -25,6 +25,7 @@ Rcpp::List extend_dummy( arma::mat sample_hyper( const int& S, const int& p, + const double& c, const arma::mat& Y, const arma::mat& X, const Rcpp::List& prior From f884acb8ee214055f100dd73525ecf90437417f3 Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Wed, 12 Jun 2024 23:27:28 +1000 Subject: [PATCH 2/4] add constants --- R/RcppExports.R | 12 +++++- src/RcppExports.cpp | 42 ++++++++++++++++--- src/sample_hyper.cpp | 99 +++++++++++++++++++++++--------------------- 3 files changed, 98 insertions(+), 55 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 2ca23a1..96bb71d 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -21,8 +21,8 @@ extend_dummy <- function(p, hyper, Y, X) { .Call(`_bsvarSIGNs_extend_dummy`, p, hyper, Y, X) } -log_dgamma <- function(x, alpha, beta) { - .Call(`_bsvarSIGNs_log_dgamma`, x, alpha, beta) +log_dgamma <- function(x, k, theta) { + .Call(`_bsvarSIGNs_log_dgamma`, x, k, theta) } log_dinvgamma <- function(x, alpha, beta) { @@ -33,6 +33,14 @@ log_prior_hyper <- function(hyper, prior) { .Call(`_bsvarSIGNs_log_prior_hyper`, hyper, prior) } +log_mvgamma <- function(n, x) { + .Call(`_bsvarSIGNs_log_mvgamma`, n, x) +} + +log_ml <- function(p, b, Omega, Psi, d, inv_Omega, Y, X) { + .Call(`_bsvarSIGNs_log_ml`, p, b, Omega, Psi, d, inv_Omega, Y, X) +} + log_ml_dummy <- function(p, hyper, Y, X, prior) { .Call(`_bsvarSIGNs_log_ml_dummy`, p, hyper, Y, X, prior) } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index ea3e5a8..4c82396 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -642,15 +642,15 @@ BEGIN_RCPP END_RCPP } // log_dgamma -double log_dgamma(const double& x, const double& alpha, const double& beta); -RcppExport SEXP _bsvarSIGNs_log_dgamma(SEXP xSEXP, SEXP alphaSEXP, SEXP betaSEXP) { +double log_dgamma(const double& x, const double& k, const double& theta); +RcppExport SEXP _bsvarSIGNs_log_dgamma(SEXP xSEXP, SEXP kSEXP, SEXP thetaSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< const double& >::type x(xSEXP); - Rcpp::traits::input_parameter< const double& >::type alpha(alphaSEXP); - Rcpp::traits::input_parameter< const double& >::type beta(betaSEXP); - rcpp_result_gen = Rcpp::wrap(log_dgamma(x, alpha, beta)); + Rcpp::traits::input_parameter< const double& >::type k(kSEXP); + Rcpp::traits::input_parameter< const double& >::type theta(thetaSEXP); + rcpp_result_gen = Rcpp::wrap(log_dgamma(x, k, theta)); return rcpp_result_gen; END_RCPP } @@ -679,6 +679,36 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// log_mvgamma +double log_mvgamma(const int& n, const double& x); +RcppExport SEXP _bsvarSIGNs_log_mvgamma(SEXP nSEXP, SEXP xSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const int& >::type n(nSEXP); + Rcpp::traits::input_parameter< const double& >::type x(xSEXP); + rcpp_result_gen = Rcpp::wrap(log_mvgamma(n, x)); + return rcpp_result_gen; +END_RCPP +} +// log_ml +double log_ml(const int& p, const arma::mat& b, const arma::mat& Omega, const arma::mat& Psi, const int& d, const arma::mat& inv_Omega, const arma::mat& Y, const arma::mat& X); +RcppExport SEXP _bsvarSIGNs_log_ml(SEXP pSEXP, SEXP bSEXP, SEXP OmegaSEXP, SEXP PsiSEXP, SEXP dSEXP, SEXP inv_OmegaSEXP, SEXP YSEXP, SEXP XSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const int& >::type p(pSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type b(bSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type Omega(OmegaSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type Psi(PsiSEXP); + Rcpp::traits::input_parameter< const int& >::type d(dSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type inv_Omega(inv_OmegaSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP); + rcpp_result_gen = Rcpp::wrap(log_ml(p, b, Omega, Psi, d, inv_Omega, Y, X)); + return rcpp_result_gen; +END_RCPP +} // log_ml_dummy double log_ml_dummy(const int& p, const arma::vec& hyper, const arma::mat& Y, const arma::mat& X, Rcpp::List prior); RcppExport SEXP _bsvarSIGNs_log_ml_dummy(SEXP pSEXP, SEXP hyperSEXP, SEXP YSEXP, SEXP XSEXP, SEXP priorSEXP) { @@ -903,6 +933,8 @@ static const R_CallMethodDef CallEntries[] = { {"_bsvarSIGNs_log_dgamma", (DL_FUNC) &_bsvarSIGNs_log_dgamma, 3}, {"_bsvarSIGNs_log_dinvgamma", (DL_FUNC) &_bsvarSIGNs_log_dinvgamma, 3}, {"_bsvarSIGNs_log_prior_hyper", (DL_FUNC) &_bsvarSIGNs_log_prior_hyper, 2}, + {"_bsvarSIGNs_log_mvgamma", (DL_FUNC) &_bsvarSIGNs_log_mvgamma, 2}, + {"_bsvarSIGNs_log_ml", (DL_FUNC) &_bsvarSIGNs_log_ml, 8}, {"_bsvarSIGNs_log_ml_dummy", (DL_FUNC) &_bsvarSIGNs_log_ml_dummy, 5}, {"_bsvarSIGNs_log_posterior_hyper", (DL_FUNC) &_bsvarSIGNs_log_posterior_hyper, 5}, {"_bsvarSIGNs_sample_hyper", (DL_FUNC) &_bsvarSIGNs_sample_hyper, 6}, diff --git a/src/sample_hyper.cpp b/src/sample_hyper.cpp index 128f4a9..268ba00 100644 --- a/src/sample_hyper.cpp +++ b/src/sample_hyper.cpp @@ -17,19 +17,19 @@ Rcpp::List update_prior( int n_hyper = hyper.n_elem; double lambda = hyper(2); - vec phi = hyper.rows(3, n_hyper - 1); + // vec psi = hyper.rows(3, n_hyper - 1); mat prior_B = as(prior["B"]); int K = prior_B.n_rows; - vec v(K); - v(K - 1) = 1e+6; - v.rows(0, K - 2) = lambda * lambda * kron(pow(linspace(1, p, p), -2), 1 / phi); + // vec v(K); + // v(K - 1) = 1e+6; + // v.rows(0, K - 2) = lambda * lambda * kron(pow(linspace(1, p, p), -2), 1 / psi); mat prior_V = lambda * lambda * as(prior["V"]); // mat prior_V = diagmat(v); mat prior_S = as(prior["S"]); - // mat prior_S = diagmat(phi); + // mat prior_S = diagmat(psi); int prior_nu = as(prior["nu"]); return List::create( @@ -75,29 +75,29 @@ Rcpp::List extend_dummy( } -// log density of gamma distribution (up to a constant) +// log density of gamma distribution // [[Rcpp:interface(cpp)]] // [[Rcpp::export]] double log_dgamma( const double& x, - const double& alpha, // shape - const double& beta // rate + const double& k, + const double& theta ) { - return (alpha - 1) * log(x) - beta * x; + return (k - 1) * log(x) - x / theta - k * log(theta) - lgamma(k); } -// log density of inverse gamma distribution (up to a constant) +// log density of inverse gamma distribution // [[Rcpp:interface(cpp)]] // [[Rcpp::export]] double log_dinvgamma( const double& x, - const double& alpha, // shape - const double& beta // scale + const double& alpha, + const double& beta ) { - return -(alpha + 1) * log(x) - beta / x; + return alpha * log(beta) - (alpha + 1) * log(x) - beta / x - lgamma(alpha); } @@ -126,39 +126,53 @@ double log_prior_hyper( } -// log marginal likelihood (up to a constant) -// notation as in Domenico, Lenza & Primiceri (2014) +// log multivariate gamma function +// [[Rcpp:interface(cpp)]] +// [[Rcpp::export]] +double log_mvgamma( + const int& n, + const double& x +) { + + if (n == 1) { + return lgamma(x); + } + + double c = (n - 1) / 2; + return c * log(M_PI) + lgamma(x - c) + log_mvgamma(n - 1, x); +} + + +// log marginal likelihood +// notation as in Giannone, Lenza & Primiceri (2014) +// [[Rcpp:interface(cpp)]] +// [[Rcpp::export]] double log_ml( - const int& p, - const mat& b, - const mat& Omega, - const mat& Phi, - const int& d, - const mat& D_Omega, - const mat& D_Phi, - const mat& inv_Omega, - const mat& Y, - const mat& X + const int& p, + const arma::mat& b, + const arma::mat& Omega, + const arma::mat& Psi, + const int& d, + const arma::mat& inv_Omega, + const arma::mat& Y, + const arma::mat& X ) { - int T = Y.n_rows + p; + int T = Y.n_rows; int N = Y.n_cols; - int K = X.n_cols; mat Bhat = solve(X.t() * X + inv_Omega, X.t() * Y + inv_Omega * b); mat ehat = Y - X * Bhat; double log_ml = 0; - // cancels out in p(Yplus) / p(Ystar) - // log_ml -= (T - p) / 2 * log_det_sympd(Phi); - - log_ml -= N / 2 * log_det(D_Omega.t() * X.t() * X * D_Omega + eye(K, K)).real(); - // log_ml -= N / 2 * log_det(X.t() * X + inv_Omega).real(); - + log_ml += - N * T / 2 * log(M_PI); + log_ml += log_mvgamma(N, (T + d) / 2) - log_mvgamma(N, d / 2); + log_ml += - N / 2 * log_det_sympd(Omega); + log_ml += d / 2 * log_det_sympd(Psi); + log_ml += - N / 2 * log_det_sympd(X.t() * X + inv_Omega); mat A = ehat.t() * ehat + (Bhat - b).t() * inv_Omega * (Bhat - b); - log_ml -= (T - p + d) / 2 * log_det(D_Phi.t() * A * D_Phi + eye(N, N)).real(); - // log_ml -= (T - p + d) / 2 * log_det(Phi + A).real(); + log_ml += - (T + d) / 2 * log_det_sympd(Psi + A); return log_ml; } @@ -178,8 +192,6 @@ double log_ml_dummy( int N = Y.n_cols; int n_hyper = hyper.n_elem; - vec phi = hyper.rows(3, n_hyper - 1); - List extended = extend_dummy(p, hyper, Y, X); mat Yplus = as(extended["Y"]); mat Xplus = as(extended["X"]); @@ -190,22 +202,13 @@ double log_ml_dummy( mat prior_S = as(prior["S"]); int prior_nu = as(prior["nu"]); - // mat chol_V = chol(prior_V, "lower"); - // mat inv_V = inv_sympd(prior_V); - mat chol_V = diagmat(sqrt(prior_V.diag())); mat inv_V = diagmat(1 / prior_V.diag()); - mat cholinv_S = diagmat(sqrt(1 / phi)); mat Ystar = Yplus.rows(0, N); mat Xstar = Xplus.rows(0, N); - double log_ml_extended = log_ml(p, prior_B, prior_V, prior_S, prior_nu, - chol_V, cholinv_S, inv_V, Yplus, Xplus); - double log_ml_dummy = log_ml(p, prior_B, prior_V, prior_S, prior_nu, - chol_V, cholinv_S, inv_V, Ystar, Xstar); - - return log_ml(p, prior_B, prior_V, prior_S, prior_nu, - chol_V, cholinv_S, inv_V, Y, X); + double log_ml_extended = log_ml(p, prior_B, prior_V, prior_S, prior_nu, inv_V, Yplus, Xplus); + double log_ml_dummy = log_ml(p, prior_B, prior_V, prior_S, prior_nu, inv_V, Ystar, Xstar); return log_ml_extended - log_ml_dummy; } From 9f313be29ba7b4394c93c92cbed4561074319fcb Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Thu, 13 Jun 2024 23:42:36 +1000 Subject: [PATCH 3/4] sample lambda in R --- R/sample_hyper.R | 100 +++++++++++++++++++++++++++++++++++++++++++ src/sample_hyper.cpp | 4 ++ 2 files changed, 104 insertions(+) create mode 100644 R/sample_hyper.R diff --git a/R/sample_hyper.R b/R/sample_hyper.R new file mode 100644 index 0000000..bcdf196 --- /dev/null +++ b/R/sample_hyper.R @@ -0,0 +1,100 @@ + +log_det = function(A) { + determinant(A, logarithm = TRUE)$modulus +} + +dlambda = function(lambda, p, Y, X) { + + N = ncol(Y) + T = nrow(Y) + prior = niw_prior(Y, p, rep(1, N), lambda) + + b = prior$B + Omega = prior$V + Psi = prior$S + d = prior$nu + + inv_Omega = solve(Omega) + Bhat = solve(t(X) %*% X + inv_Omega) %*% (t(X) %*% Y + inv_Omega %*% b) + ehat = Y - X %*% Bhat + + lprior = dgamma( + lambda, + shape = (9 + sqrt(17)) / 8, + scale = 8 / (5 * (1 + sqrt(17))), + log = TRUE + ) + + llike = -N / 2 * log_det(Omega) + llike = llike - N / 2 * log_det(t(X) %*% X + inv_Omega) + A = Psi + t(ehat) %*% ehat + t(Bhat - b) %*% inv_Omega %*% (Bhat - b) + llike = llike - (T + d) / 2 * log_det(A) + + return(lprior + llike) +} + + +sample_lambda = function(S, p, Y, X) { + init = 0.2 + result = optim( + init, + \(lambda) { + nlogp = -dlambda(lambda, p, Y, X) + if (!is.finite(nlogp)) { + nlogp = 1e+6 + } + return(nlogp) + }, + method = "L-BFGS-B", + lower = init / 100, + upper = init * 100, + hessian = TRUE, + control = list(trace = 1, maxit = 1e4) + ) + + lambda = numeric(S) + lambda[1] = result$par + lambda[1] = 0.2 + dlambda = result$value + + c = 0.5 + W = c * 1 / result$hessian + + for (s in 2:S) { + new_lambda = rnorm(1, lambda[s - 1], sqrt(W)) + new_dlambda = dlambda(new_lambda, p, Y, X) + + if (runif(1) < exp(new_dlambda - dlambda)) { + lambda[s] = new_lambda + dlambda = new_dlambda + } else { + lambda[s] = lambda[s - 1] + } + } + + return(lambda) +} + + +library(R.matlab) +data = readMat('../data/DataSW.mat') +name = sapply(1:7, \(i) data$ShortDescr[[i]][[1]]) +data = data$y +colnames(data) = name +data = data[, c(1, 2, 7)] + +p = 4 +d = bsvars::specify_data_matrices$new(data, p, NULL) +Y = t(d$Y) +X = t(d$X) + +lambda = sample_lambda(1000, p, Y, X) + +par(mfrow = c(2, 1)) +plot(lambda, type = 'l') +hist(lambda, breaks = 100) + + + + + diff --git a/src/sample_hyper.cpp b/src/sample_hyper.cpp index 268ba00..73c9091 100644 --- a/src/sample_hyper.cpp +++ b/src/sample_hyper.cpp @@ -255,6 +255,7 @@ arma::mat sample_hyper( mat posterior_hyper(n_hyper, S); + int success = 0; for (int s = 0; s < S; s++) { // sample from proposal distribution @@ -265,11 +266,14 @@ arma::mat sample_hyper( if (randu() < r) { hyper = hyperp; logp = logpp; + success++; } posterior_hyper.col(s) = hyper; } + cout << "Acceptance rate: " << (double) success / S << endl; + return posterior_hyper; } From 8fc1b687554a0248229aabcf5daadbe1359b62ec Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Tue, 18 Jun 2024 09:06:03 +1000 Subject: [PATCH 4/4] hyperparameters in prior --- R/RcppExports.R | 28 ++-- R/estimate.BSVARSIGN.R | 2 +- R/sample_hyper.R | 100 ------------- R/specify_bsvarSIGN.R | 41 +++--- R/utils.R | 2 +- inst/include/bsvarSIGNs_RcppExports.h | 6 +- src/RcppExports.cpp | 90 ++++++------ src/bsvars_sign.cpp | 38 ++--- src/bsvars_sign.h | 4 +- src/sample_hyper.cpp | 199 ++++++++++++++------------ src/sample_hyper.h | 33 +++-- 11 files changed, 238 insertions(+), 305 deletions(-) delete mode 100644 R/sample_hyper.R diff --git a/R/RcppExports.R b/R/RcppExports.R index 96bb71d..12cb357 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -13,12 +13,16 @@ niw_cpp <- function(Y, X, prior) { .Call(`_bsvarSIGNs_niw_cpp`, Y, X, prior) } -update_prior <- function(p, hyper, prior) { - .Call(`_bsvarSIGNs_update_prior`, p, hyper, prior) +mn_prior <- function(p, lambda, psi) { + .Call(`_bsvarSIGNs_mn_prior`, p, lambda, psi) } -extend_dummy <- function(p, hyper, Y, X) { - .Call(`_bsvarSIGNs_extend_dummy`, p, hyper, Y, X) +update_prior <- function(p, hyper, model, prior) { + .Call(`_bsvarSIGNs_update_prior`, p, hyper, model, prior) +} + +extend_dummy <- function(p, hyper, model, Y, X) { + .Call(`_bsvarSIGNs_extend_dummy`, p, hyper, model, Y, X) } log_dgamma <- function(x, k, theta) { @@ -29,8 +33,8 @@ log_dinvgamma <- function(x, alpha, beta) { .Call(`_bsvarSIGNs_log_dinvgamma`, x, alpha, beta) } -log_prior_hyper <- function(hyper, prior) { - .Call(`_bsvarSIGNs_log_prior_hyper`, hyper, prior) +log_prior_hyper <- function(hyper, model, prior) { + .Call(`_bsvarSIGNs_log_prior_hyper`, hyper, model, prior) } log_mvgamma <- function(n, x) { @@ -41,16 +45,12 @@ log_ml <- function(p, b, Omega, Psi, d, inv_Omega, Y, X) { .Call(`_bsvarSIGNs_log_ml`, p, b, Omega, Psi, d, inv_Omega, Y, X) } -log_ml_dummy <- function(p, hyper, Y, X, prior) { - .Call(`_bsvarSIGNs_log_ml_dummy`, p, hyper, Y, X, prior) -} - -log_posterior_hyper <- function(p, hyper, Y, X, prior) { - .Call(`_bsvarSIGNs_log_posterior_hyper`, p, hyper, Y, X, prior) +log_ml_dummy <- function(p, hyper, model, Y, X, prior) { + .Call(`_bsvarSIGNs_log_ml_dummy`, p, hyper, model, Y, X, prior) } -sample_hyper <- function(S, p, c, Y, X, prior) { - .Call(`_bsvarSIGNs_sample_hyper`, S, p, c, Y, X, prior) +log_posterior_hyper <- function(p, hyper, model, Y, X, prior) { + .Call(`_bsvarSIGNs_log_posterior_hyper`, p, hyper, model, Y, X, prior) } # Register entry points for exported C++ functions diff --git a/R/estimate.BSVARSIGN.R b/R/estimate.BSVARSIGN.R index c76bed2..148800f 100644 --- a/R/estimate.BSVARSIGN.R +++ b/R/estimate.BSVARSIGN.R @@ -108,7 +108,7 @@ estimate.BSVARSIGN <- function(specification, S, thin = 1, show_progress = TRUE) data_matrices$Y, data_matrices$X, identification$VB, identification$sign_irf, identification$sign_narrative, identification$sign_B, identification$zero_irf, - prior, starting_values, thin, show_progress, max_tries) + prior, starting_values, show_progress, thin, max_tries) specification$starting_values$set_starting_values(qqq$last_draw) output = specify_posterior_bsvarSIGN$new(specification, qqq$posterior) diff --git a/R/sample_hyper.R b/R/sample_hyper.R deleted file mode 100644 index bcdf196..0000000 --- a/R/sample_hyper.R +++ /dev/null @@ -1,100 +0,0 @@ - -log_det = function(A) { - determinant(A, logarithm = TRUE)$modulus -} - -dlambda = function(lambda, p, Y, X) { - - N = ncol(Y) - T = nrow(Y) - prior = niw_prior(Y, p, rep(1, N), lambda) - - b = prior$B - Omega = prior$V - Psi = prior$S - d = prior$nu - - inv_Omega = solve(Omega) - Bhat = solve(t(X) %*% X + inv_Omega) %*% (t(X) %*% Y + inv_Omega %*% b) - ehat = Y - X %*% Bhat - - lprior = dgamma( - lambda, - shape = (9 + sqrt(17)) / 8, - scale = 8 / (5 * (1 + sqrt(17))), - log = TRUE - ) - - llike = -N / 2 * log_det(Omega) - llike = llike - N / 2 * log_det(t(X) %*% X + inv_Omega) - A = Psi + t(ehat) %*% ehat + t(Bhat - b) %*% inv_Omega %*% (Bhat - b) - llike = llike - (T + d) / 2 * log_det(A) - - return(lprior + llike) -} - - -sample_lambda = function(S, p, Y, X) { - init = 0.2 - result = optim( - init, - \(lambda) { - nlogp = -dlambda(lambda, p, Y, X) - if (!is.finite(nlogp)) { - nlogp = 1e+6 - } - return(nlogp) - }, - method = "L-BFGS-B", - lower = init / 100, - upper = init * 100, - hessian = TRUE, - control = list(trace = 1, maxit = 1e4) - ) - - lambda = numeric(S) - lambda[1] = result$par - lambda[1] = 0.2 - dlambda = result$value - - c = 0.5 - W = c * 1 / result$hessian - - for (s in 2:S) { - new_lambda = rnorm(1, lambda[s - 1], sqrt(W)) - new_dlambda = dlambda(new_lambda, p, Y, X) - - if (runif(1) < exp(new_dlambda - dlambda)) { - lambda[s] = new_lambda - dlambda = new_dlambda - } else { - lambda[s] = lambda[s - 1] - } - } - - return(lambda) -} - - -library(R.matlab) -data = readMat('../data/DataSW.mat') -name = sapply(1:7, \(i) data$ShortDescr[[i]][[1]]) -data = data$y -colnames(data) = name -data = data[, c(1, 2, 7)] - -p = 4 -d = bsvars::specify_data_matrices$new(data, p, NULL) -Y = t(d$Y) -X = t(d$X) - -lambda = sample_lambda(1000, p, Y, X) - -par(mfrow = c(2, 1)) -plot(lambda, type = 'l') -hist(lambda, breaks = 100) - - - - - diff --git a/R/specify_bsvarSIGN.R b/R/specify_bsvarSIGN.R index f72b646..ffff6e2 100644 --- a/R/specify_bsvarSIGN.R +++ b/R/specify_bsvarSIGN.R @@ -124,21 +124,12 @@ specify_prior_bsvarSIGN = R6::R6Class( "PriorBSVARSIGN", public = list( + #' @field hyper a \code{(N+3)xS} matrix of posterior draws of hyperparameters. + hyper = matrix(), - #' @field B an \code{KxN} matrix, the mean of the normal prior distribution for the parameter matrix \eqn{B}. - B = matrix(), - - #' @field V a \code{KxK} covariance matrix of the normal prior distribution for each of - #' the column of the parameter matrix \eqn{B}. This covariance matrix is equation invariant. - V = matrix(), - - #' @field S an \code{NxN} scale matrix of the inverse-Wishart prior distribution - #' for the covariance matrix \eqn{\Sigma}. This scale matrix is equation invariant. - S = matrix(), - - #' @field nu a positive real number greater of equal than \code{N}, a degree of freedom parameter - #' of the inverse-Wishart prior distribution for the covariance matrix \eqn{\Sigma}. - nu = NA, + #' @field model a \code{4x1} vector of Boolean values indicating prior specifications, + #' model[1] = dummy soc, model[2] = dummy sur, model[3] = mn lambda, model[4] = mn psi. + model = c(), #' @description #' Create a new prior specification PriorBSVAR. @@ -162,11 +153,17 @@ specify_prior_bsvarSIGN = R6::R6Class( stopifnot("Argument d must be a non-negative integer number." = d >= 0 & d %% 1 == 0) stopifnot("Argument stationary must be a logical vector of length N." = length(stationary) == N & is.logical(stationary)) - prior = niw_prior(data, p, !stationary) - self$B = prior$B - self$V = prior$V - self$S = prior$S - self$nu = prior$nu + Y = data + T = nrow(Y) + lambda = 0.2 + psi = sapply(1:N, \(i) summary(lm(Y[2:T, i] ~ Y[1:(T - 1), i]))$sigma^2) + + hyper = matrix(NA, N + 3, 1) + hyper[3, ] = lambda + hyper[4:(N + 3), ] = psi + + self$model = c(FALSE, FALSE, FALSE, FALSE) + self$hyper = hyper }, # END initialize @@ -180,10 +177,8 @@ specify_prior_bsvarSIGN = R6::R6Class( #' get_prior = function(){ list( - B = self$B, - V = self$V, - S = self$S, - nu = self$nu + model = self$model, + hyper = self$hyper ) } # END get_prior diff --git a/R/utils.R b/R/utils.R index b05d974..a084d9e 100644 --- a/R/utils.R +++ b/R/utils.R @@ -13,7 +13,7 @@ importance_sampling <- function(posterior) { posterior$posterior$A = posterior$posterior$A[, , indices] posterior$posterior$B = posterior$posterior$B[, , indices] - posterior$posterior$hyper = posterior$posterior$hyper[, , indices] + posterior$posterior$hyper = posterior$posterior$hyper[, indices] posterior$posterior$Q = posterior$posterior$Q[, , indices] posterior$posterior$Sigma = posterior$posterior$Sigma[, , indices] posterior$posterior$Theta0 = posterior$posterior$Theta0[, , indices] diff --git a/inst/include/bsvarSIGNs_RcppExports.h b/inst/include/bsvarSIGNs_RcppExports.h index c08c3b9..ab1d075 100644 --- a/inst/include/bsvarSIGNs_RcppExports.h +++ b/inst/include/bsvarSIGNs_RcppExports.h @@ -88,17 +88,17 @@ namespace bsvarSIGNs { return Rcpp::as >(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& VB, const arma::cube& sign_irf, const arma::mat& sign_narrative, const arma::mat& sign_B, const arma::field& Z, const Rcpp::List& prior, const Rcpp::List& starting_values, const int thin = 100, const bool show_progress = true, const int& max_tries = 10000) { + inline Rcpp::List bsvar_sign_cpp(const int& S, const int& lags, const arma::mat& Y, const arma::mat& X, const arma::field& VB, const arma::cube& sign_irf, const arma::mat& sign_narrative, const arma::mat& sign_B, const arma::field& 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&,const arma::cube&,const arma::mat&,const arma::mat&,const arma::field&,const Rcpp::List&,const Rcpp::List&,const int,const bool,const int&)"); + validateSignature("Rcpp::List(*bsvar_sign_cpp)(const int&,const int&,const arma::mat&,const arma::mat&,const arma::field&,const arma::cube&,const arma::mat&,const arma::mat&,const arma::field&,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_bsvar_sign_cpp(Shield(Rcpp::wrap(S)), Shield(Rcpp::wrap(lags)), Shield(Rcpp::wrap(Y)), Shield(Rcpp::wrap(X)), Shield(Rcpp::wrap(VB)), Shield(Rcpp::wrap(sign_irf)), Shield(Rcpp::wrap(sign_narrative)), Shield(Rcpp::wrap(sign_B)), Shield(Rcpp::wrap(Z)), Shield(Rcpp::wrap(prior)), Shield(Rcpp::wrap(starting_values)), Shield(Rcpp::wrap(thin)), Shield(Rcpp::wrap(show_progress)), Shield(Rcpp::wrap(max_tries))); + rcpp_result_gen = p_bsvar_sign_cpp(Shield(Rcpp::wrap(S)), Shield(Rcpp::wrap(lags)), Shield(Rcpp::wrap(Y)), Shield(Rcpp::wrap(X)), Shield(Rcpp::wrap(VB)), Shield(Rcpp::wrap(sign_irf)), Shield(Rcpp::wrap(sign_narrative)), Shield(Rcpp::wrap(sign_B)), Shield(Rcpp::wrap(Z)), Shield(Rcpp::wrap(prior)), Shield(Rcpp::wrap(starting_values)), Shield(Rcpp::wrap(show_progress)), Shield(Rcpp::wrap(thin)), Shield(Rcpp::wrap(max_tries))); } if (rcpp_result_gen.inherits("interrupted-error")) throw Rcpp::internal::InterruptedException(); diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 4c82396..1dd9525 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -127,8 +127,8 @@ 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& VB, const arma::cube& sign_irf, const arma::mat& sign_narrative, const arma::mat& sign_B, const arma::field& Z, const Rcpp::List& prior, const Rcpp::List& starting_values, const int thin, const bool show_progress, 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 thinSEXP, SEXP show_progressSEXP, SEXP max_triesSEXP) { +Rcpp::List bsvar_sign_cpp(const int& S, const int& lags, const arma::mat& Y, const arma::mat& X, const arma::field& VB, const arma::cube& sign_irf, const arma::mat& sign_narrative, const arma::mat& sign_B, const arma::field& 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) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::traits::input_parameter< const int& >::type S(SSEXP); @@ -142,18 +142,18 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const arma::field& >::type Z(ZSEXP); Rcpp::traits::input_parameter< const Rcpp::List& >::type prior(priorSEXP); Rcpp::traits::input_parameter< const Rcpp::List& >::type starting_values(starting_valuesSEXP); - Rcpp::traits::input_parameter< const int >::type thin(thinSEXP); 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, thin, show_progress, max_tries)); + 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)); 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 thinSEXP, SEXP show_progressSEXP, SEXP max_triesSEXP) { +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) { 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, thinSEXP, show_progressSEXP, max_triesSEXP)); + 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)); } Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error"); if (rcpp_isInterrupt_gen) { @@ -614,30 +614,45 @@ RcppExport SEXP _bsvarSIGNs_sample_Q(SEXP pSEXP, SEXP YSEXP, SEXP XSEXP, SEXP BS UNPROTECT(1); return rcpp_result_gen; } +// mn_prior +Rcpp::List mn_prior(const int& p, const double& lambda, const arma::vec& psi); +RcppExport SEXP _bsvarSIGNs_mn_prior(SEXP pSEXP, SEXP lambdaSEXP, SEXP psiSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const int& >::type p(pSEXP); + Rcpp::traits::input_parameter< const double& >::type lambda(lambdaSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type psi(psiSEXP); + rcpp_result_gen = Rcpp::wrap(mn_prior(p, lambda, psi)); + return rcpp_result_gen; +END_RCPP +} // update_prior -Rcpp::List update_prior(const int& p, const arma::vec& hyper, const Rcpp::List& prior); -RcppExport SEXP _bsvarSIGNs_update_prior(SEXP pSEXP, SEXP hyperSEXP, SEXP priorSEXP) { +Rcpp::List update_prior(const int& p, const arma::vec& hyper, const arma::vec& model, const Rcpp::List& prior); +RcppExport SEXP _bsvarSIGNs_update_prior(SEXP pSEXP, SEXP hyperSEXP, SEXP modelSEXP, SEXP priorSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< const int& >::type p(pSEXP); Rcpp::traits::input_parameter< const arma::vec& >::type hyper(hyperSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type model(modelSEXP); Rcpp::traits::input_parameter< const Rcpp::List& >::type prior(priorSEXP); - rcpp_result_gen = Rcpp::wrap(update_prior(p, hyper, prior)); + rcpp_result_gen = Rcpp::wrap(update_prior(p, hyper, model, prior)); return rcpp_result_gen; END_RCPP } // extend_dummy -Rcpp::List extend_dummy(const int& p, const arma::vec& hyper, const arma::mat& Y, const arma::mat& X); -RcppExport SEXP _bsvarSIGNs_extend_dummy(SEXP pSEXP, SEXP hyperSEXP, SEXP YSEXP, SEXP XSEXP) { +Rcpp::List extend_dummy(const int& p, const arma::vec& hyper, const arma::vec& model, const arma::mat& Y, const arma::mat& X); +RcppExport SEXP _bsvarSIGNs_extend_dummy(SEXP pSEXP, SEXP hyperSEXP, SEXP modelSEXP, SEXP YSEXP, SEXP XSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< const int& >::type p(pSEXP); Rcpp::traits::input_parameter< const arma::vec& >::type hyper(hyperSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type model(modelSEXP); Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP); Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP); - rcpp_result_gen = Rcpp::wrap(extend_dummy(p, hyper, Y, X)); + rcpp_result_gen = Rcpp::wrap(extend_dummy(p, hyper, model, Y, X)); return rcpp_result_gen; END_RCPP } @@ -668,14 +683,15 @@ BEGIN_RCPP END_RCPP } // log_prior_hyper -double log_prior_hyper(const arma::vec& hyper, const Rcpp::List& prior); -RcppExport SEXP _bsvarSIGNs_log_prior_hyper(SEXP hyperSEXP, SEXP priorSEXP) { +double log_prior_hyper(const arma::vec& hyper, const arma::vec& model, const Rcpp::List& prior); +RcppExport SEXP _bsvarSIGNs_log_prior_hyper(SEXP hyperSEXP, SEXP modelSEXP, SEXP priorSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< const arma::vec& >::type hyper(hyperSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type model(modelSEXP); Rcpp::traits::input_parameter< const Rcpp::List& >::type prior(priorSEXP); - rcpp_result_gen = Rcpp::wrap(log_prior_hyper(hyper, prior)); + rcpp_result_gen = Rcpp::wrap(log_prior_hyper(hyper, model, prior)); return rcpp_result_gen; END_RCPP } @@ -710,48 +726,34 @@ BEGIN_RCPP END_RCPP } // log_ml_dummy -double log_ml_dummy(const int& p, const arma::vec& hyper, const arma::mat& Y, const arma::mat& X, Rcpp::List prior); -RcppExport SEXP _bsvarSIGNs_log_ml_dummy(SEXP pSEXP, SEXP hyperSEXP, SEXP YSEXP, SEXP XSEXP, SEXP priorSEXP) { +double log_ml_dummy(const int& p, const arma::vec& hyper, const arma::vec& model, const arma::mat& Y, const arma::mat& X, Rcpp::List prior); +RcppExport SEXP _bsvarSIGNs_log_ml_dummy(SEXP pSEXP, SEXP hyperSEXP, SEXP modelSEXP, SEXP YSEXP, SEXP XSEXP, SEXP priorSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< const int& >::type p(pSEXP); Rcpp::traits::input_parameter< const arma::vec& >::type hyper(hyperSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type model(modelSEXP); Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP); Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP); Rcpp::traits::input_parameter< Rcpp::List >::type prior(priorSEXP); - rcpp_result_gen = Rcpp::wrap(log_ml_dummy(p, hyper, Y, X, prior)); + rcpp_result_gen = Rcpp::wrap(log_ml_dummy(p, hyper, model, Y, X, prior)); return rcpp_result_gen; END_RCPP } // log_posterior_hyper -double log_posterior_hyper(const int& p, const arma::vec& hyper, const arma::mat& Y, const arma::mat& X, const Rcpp::List& prior); -RcppExport SEXP _bsvarSIGNs_log_posterior_hyper(SEXP pSEXP, SEXP hyperSEXP, SEXP YSEXP, SEXP XSEXP, SEXP priorSEXP) { +double log_posterior_hyper(const int& p, const arma::vec& hyper, const arma::vec& model, const arma::mat& Y, const arma::mat& X, const Rcpp::List& prior); +RcppExport SEXP _bsvarSIGNs_log_posterior_hyper(SEXP pSEXP, SEXP hyperSEXP, SEXP modelSEXP, SEXP YSEXP, SEXP XSEXP, SEXP priorSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< const int& >::type p(pSEXP); Rcpp::traits::input_parameter< const arma::vec& >::type hyper(hyperSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type model(modelSEXP); 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 Rcpp::List& >::type prior(priorSEXP); - rcpp_result_gen = Rcpp::wrap(log_posterior_hyper(p, hyper, Y, X, prior)); - return rcpp_result_gen; -END_RCPP -} -// sample_hyper -arma::mat sample_hyper(const int& S, const int& p, const double& c, const arma::mat& Y, const arma::mat& X, const Rcpp::List& prior); -RcppExport SEXP _bsvarSIGNs_sample_hyper(SEXP SSEXP, SEXP pSEXP, SEXP cSEXP, SEXP YSEXP, SEXP XSEXP, SEXP priorSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< const int& >::type S(SSEXP); - Rcpp::traits::input_parameter< const int& >::type p(pSEXP); - Rcpp::traits::input_parameter< const double& >::type c(cSEXP); - 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 Rcpp::List& >::type prior(priorSEXP); - rcpp_result_gen = Rcpp::wrap(sample_hyper(S, p, c, Y, X, prior)); + rcpp_result_gen = Rcpp::wrap(log_posterior_hyper(p, hyper, model, Y, X, prior)); return rcpp_result_gen; END_RCPP } @@ -866,7 +868,7 @@ static int _bsvarSIGNs_RcppExport_validate(const char* sig) { signatures.insert("arma::mat(*hd1_cpp)(const int&,const int&,const int&,const arma::mat&,const arma::cube&)"); signatures.insert("arma::cube(*ir1_cpp)(const arma::mat&,const arma::mat&,int,const int&)"); signatures.insert("arma::field(*bsvarSIGNs_ir)(arma::cube&,arma::cube&,const int,const int)"); - signatures.insert("Rcpp::List(*bsvar_sign_cpp)(const int&,const int&,const arma::mat&,const arma::mat&,const arma::field&,const arma::cube&,const arma::mat&,const arma::mat&,const arma::field&,const Rcpp::List&,const Rcpp::List&,const int,const bool,const int&)"); + signatures.insert("Rcpp::List(*bsvar_sign_cpp)(const int&,const int&,const arma::mat&,const arma::mat&,const arma::field&,const arma::cube&,const arma::mat&,const arma::mat&,const arma::field&,const Rcpp::List&,const Rcpp::List&,const bool,const int,const int&)"); signatures.insert("bool(*match_sign_narrative)(const arma::mat&,const arma::mat&,const arma::cube&)"); signatures.insert("double(*weight_narrative)(const int&,arma::mat,const arma::cube&)"); signatures.insert("arma::field(*ZIRF)(const arma::field&,const arma::mat&)"); @@ -928,16 +930,16 @@ static const R_CallMethodDef CallEntries[] = { {"_bsvarSIGNs_niw_cpp", (DL_FUNC) &_bsvarSIGNs_niw_cpp, 3}, {"_bsvarSIGNs_match_sign_irf", (DL_FUNC) &_bsvarSIGNs_match_sign_irf, 3}, {"_bsvarSIGNs_sample_Q", (DL_FUNC) &_bsvarSIGNs_sample_Q, 13}, - {"_bsvarSIGNs_update_prior", (DL_FUNC) &_bsvarSIGNs_update_prior, 3}, - {"_bsvarSIGNs_extend_dummy", (DL_FUNC) &_bsvarSIGNs_extend_dummy, 4}, + {"_bsvarSIGNs_mn_prior", (DL_FUNC) &_bsvarSIGNs_mn_prior, 3}, + {"_bsvarSIGNs_update_prior", (DL_FUNC) &_bsvarSIGNs_update_prior, 4}, + {"_bsvarSIGNs_extend_dummy", (DL_FUNC) &_bsvarSIGNs_extend_dummy, 5}, {"_bsvarSIGNs_log_dgamma", (DL_FUNC) &_bsvarSIGNs_log_dgamma, 3}, {"_bsvarSIGNs_log_dinvgamma", (DL_FUNC) &_bsvarSIGNs_log_dinvgamma, 3}, - {"_bsvarSIGNs_log_prior_hyper", (DL_FUNC) &_bsvarSIGNs_log_prior_hyper, 2}, + {"_bsvarSIGNs_log_prior_hyper", (DL_FUNC) &_bsvarSIGNs_log_prior_hyper, 3}, {"_bsvarSIGNs_log_mvgamma", (DL_FUNC) &_bsvarSIGNs_log_mvgamma, 2}, {"_bsvarSIGNs_log_ml", (DL_FUNC) &_bsvarSIGNs_log_ml, 8}, - {"_bsvarSIGNs_log_ml_dummy", (DL_FUNC) &_bsvarSIGNs_log_ml_dummy, 5}, - {"_bsvarSIGNs_log_posterior_hyper", (DL_FUNC) &_bsvarSIGNs_log_posterior_hyper, 5}, - {"_bsvarSIGNs_sample_hyper", (DL_FUNC) &_bsvarSIGNs_sample_hyper, 6}, + {"_bsvarSIGNs_log_ml_dummy", (DL_FUNC) &_bsvarSIGNs_log_ml_dummy, 6}, + {"_bsvarSIGNs_log_posterior_hyper", (DL_FUNC) &_bsvarSIGNs_log_posterior_hyper, 6}, {"_bsvarSIGNs_qr_sign_cpp", (DL_FUNC) &_bsvarSIGNs_qr_sign_cpp, 1}, {"_bsvarSIGNs_rortho_cpp", (DL_FUNC) &_bsvarSIGNs_rortho_cpp, 1}, {"_bsvarSIGNs_match_sign", (DL_FUNC) &_bsvarSIGNs_match_sign, 2}, diff --git a/src/bsvars_sign.cpp b/src/bsvars_sign.cpp index 2d6e8e8..825b856 100644 --- a/src/bsvars_sign.cpp +++ b/src/bsvars_sign.cpp @@ -5,6 +5,7 @@ #include +#include "sample_hyper.h" #include "sample_Q.h" #include "sample_NIW.h" @@ -26,8 +27,8 @@ Rcpp::List bsvar_sign_cpp( const arma::field& Z, // a list of zero restrictions const Rcpp::List& prior, // a list of priors const Rcpp::List& starting_values, // a list of starting values - const int thin = 100, // introduce thinning const bool show_progress = true, + const int thin = 100, // introduce thinning const int& max_tries = 10000 // maximum tries for Q draw ) { @@ -51,9 +52,6 @@ Rcpp::List bsvar_sign_cpp( } Progress p(50, show_progress); - mat Yt = Y.t(); - mat Xt = X.t(); - const int T = Y.n_rows; const int N = Y.n_cols; const int K = X.n_cols; @@ -64,30 +62,37 @@ Rcpp::List bsvar_sign_cpp( vec posterior_w(S); + vec model = as(prior["model"]); + mat sample_hyper = as(prior["hyper"]); + + mat posterior_hyper(N + 3, S); cube posterior_A(N, K, S); cube posterior_B(N, N, S); cube posterior_Sigma(N, N, S); cube posterior_Theta0(N, N, S); - cube posterior_hyper(2 * N + 1, 2, S); cube posterior_shocks(N, T, S); - double w = 1; - + int s = 0, S_hyper = sample_hyper.n_cols - 1, post_nu; + double w = 1, lambda; + vec hyper, psi; + mat B, Sigma, chol_Sigma, h_invp, Q, Epsilon, post_B, post_V, post_S; + List result; - List result = niw_cpp(Y, X, prior); - mat post_B = as(result["B"]); - mat post_V = as(result["V"]); - mat post_S = as(result["S"]); - int post_nu = as(result["nu"]); - - mat B, Sigma, chol_Sigma, h_invp, Q, Epsilon; - - int s = 0; while (s < S) { // Check for user interrupts if (s % 200 == 0) checkUserInterrupt(); + hyper = sample_hyper.col(randi(distr_param(0, S_hyper))); + lambda = hyper(2); + psi = hyper.rows(3, N + 2); + + result = niw_cpp(Y, X, mn_prior(lags, lambda, psi)); + post_B = as(result["B"]); + post_V = as(result["V"]); + post_S = as(result["S"]); + post_nu = as(result["nu"]); + Sigma = iwishrnd(post_S, post_nu); chol_Sigma = chol(Sigma, "lower"); B = rmatnorm_cpp(post_B, post_V, Sigma); @@ -112,6 +117,7 @@ Rcpp::List bsvar_sign_cpp( posterior_Sigma.slice(s) = Sigma; posterior_Theta0.slice(s) = chol_Sigma * Q; posterior_shocks.slice(s) = Epsilon; + posterior_hyper.col(s) = hyper; s++; } diff --git a/src/bsvars_sign.h b/src/bsvars_sign.h index 65ab5a8..4a779c1 100644 --- a/src/bsvars_sign.h +++ b/src/bsvars_sign.h @@ -21,11 +21,11 @@ Rcpp::List bsvar_sign_cpp( const arma::cube& sign_irf, // NxNxh cube of signs for impulse response function const arma::mat& sign_narrative, // Mx6 matrix of signs for historical decomposition const arma::mat& sign_B, // Mx6 matrix of signs for B - const Rcpp::List& Z, // a list of zero restrictions + const arma::field& Z, // a list of zero restrictions const Rcpp::List& prior, // a list of priors const Rcpp::List& starting_values, // a list of starting values - const int thin = 100, // introduce thinning const bool show_progress = true, + const int thin = 100, // introduce thinning const int& max_tries = 10000 // maximum tries for Q draw ); diff --git a/src/sample_hyper.cpp b/src/sample_hyper.cpp index 73c9091..a48a54e 100644 --- a/src/sample_hyper.cpp +++ b/src/sample_hyper.cpp @@ -5,32 +5,74 @@ using namespace Rcpp; using namespace arma; +// get prior from hyper-parameters +// [[Rcpp:interface(cpp)]] +// [[Rcpp::export]] +Rcpp::List mn_prior( + const int& p, + const double& lambda, + const arma::vec& psi +) { + + int N = psi.n_elem; + int K = 1 + N * p; + + mat B(K, N, fill::eye); + + vec v(K); + v.rows(0, K - 2) = kron(pow(linspace(1, p, p), -2), 1 / psi); + mat V = lambda * lambda * diagmat(v); + V(K - 1, K - 1) = 1e+6; + + mat S = diagmat(psi); + int nu = N + 2; + + return List::create( + _["B"] = B, + _["V"] = V, + _["S"] = S, + _["nu"] = nu + ); +} + + // update prior with hyper-parameters // [[Rcpp:interface(cpp)]] // [[Rcpp::export]] Rcpp::List update_prior( const int& p, const arma::vec& hyper, + const arma::vec& model, const Rcpp::List& prior ) { - int n_hyper = hyper.n_elem; - - double lambda = hyper(2); - // vec psi = hyper.rows(3, n_hyper - 1); - mat prior_B = as(prior["B"]); + int prior_nu = as(prior["nu"]); - int K = prior_B.n_rows; - // vec v(K); - // v(K - 1) = 1e+6; - // v.rows(0, K - 2) = lambda * lambda * kron(pow(linspace(1, p, p), -2), 1 / psi); - mat prior_V = lambda * lambda * as(prior["V"]); - // mat prior_V = diagmat(v); + mat prior_V, prior_S; + if (model(3)) { + int K = prior_B.n_rows; + int n_hyper = hyper.n_elem; + + vec psi = hyper.rows(3, n_hyper - 1); + prior_S = diagmat(psi); + + vec v(K); + v(K - 1) = 1e+6; + v.rows(0, K - 2) = kron(pow(linspace(1, p, p), -2), 1 / psi); + prior_V = diagmat(v); + } else { + prior_V = as(prior["V"]); + prior_S = as(prior["S"]); + } - mat prior_S = as(prior["S"]); - // mat prior_S = diagmat(psi); - int prior_nu = as(prior["nu"]); + double lambda; + if (model(2)) { + lambda = hyper(2); + } else { + lambda = 0.2; + } + prior_V *= lambda * lambda; return List::create( _["B"] = prior_B, @@ -47,30 +89,45 @@ Rcpp::List update_prior( Rcpp::List extend_dummy( const int& p, const arma::vec& hyper, + const arma::vec& model, const arma::mat& Y, const arma::mat& X ) { int N = Y.n_cols; - - double mu = hyper(0); - double delta = hyper(1); + int K = X.n_cols; mat ybar0 = mean(X.submat(0, 0, p - 1, N - 1), 0); - mat yplus = diagmat(ybar0 / mu); - mat ypplus = ybar0 / delta; - mat Ystar = join_vert(ypplus, yplus); - mat Yplus = join_vert(Ystar, Y); + double mu, delta; + mat Ystar(0, N), Xstar(0, K); + + if (model(0)) { + mu = hyper(0); + mat yp = diagmat(ybar0 / mu); + mat xp = join_horiz(repmat(yp, 1, p), zeros(N, 1)); + + Ystar = join_vert(Ystar, yp); + Xstar = join_vert(Xstar, xp); + } + + if (model(1)) { + delta = hyper(1); + mat ypp = ybar0 / delta; + mat xpp = join_horiz(repmat(ypp, 1, p), mat(1, 1, fill::value(1 / delta))); + + Ystar = join_vert(Ystar, ypp); + Xstar = join_vert(Xstar, xpp); + } - mat Xstar = zeros(N + 1, 1); - Xstar(0, 0) = 1 / delta; - Xstar = join_horiz(repmat(Ystar, 1, p), Xstar); + mat Yplus = join_vert(Ystar, Y); mat Xplus = join_vert(Xstar, X); return List::create( - _["Y"] = Yplus, - _["X"] = Xplus + _["Yplus"] = Yplus, + _["Xplus"] = Xplus, + _["Ystar"] = Ystar, + _["Xstar"] = Xstar ); } @@ -101,11 +158,12 @@ double log_dinvgamma( } -// log prior density of hyper-parameters (up to a constant) +// log prior density of hyper-parameters // [[Rcpp:interface(cpp)]] // [[Rcpp::export]] double log_prior_hyper( const arma::vec& hyper, + const arma::vec& model, const Rcpp::List& prior ) { @@ -113,11 +171,20 @@ double log_prior_hyper( mat prior_hyper = prior["hyper"]; - for (int i = 0; i < hyper.n_elem; i++) { - // first 3 hyper-parameters ~ gamma - if (i < 3) { - log_prior += log_dgamma(hyper(i), prior_hyper(i, 0), prior_hyper(i, 1)); - } else { + if (model(0)) { + log_prior += log_dgamma(hyper(0), prior_hyper(0, 0), prior_hyper(0, 1)); + } + + if (model(1)) { + log_prior += log_dgamma(hyper(1), prior_hyper(1, 0), prior_hyper(1, 1)); + } + + if (model(2)) { + log_prior += log_dgamma(hyper(2), prior_hyper(2, 0), prior_hyper(2, 1)); + } + + if (model(3)) { + for (int i = 3; i < hyper.n_elem; i++) { log_prior += log_dinvgamma(hyper(i), prior_hyper(i, 0), prior_hyper(i, 1)); } } @@ -178,12 +245,13 @@ double log_ml( } -// log marginal likelihood with dummy observations (up to a constant) +// log marginal likelihood with dummy observations // [[Rcpp:interface(cpp)]] // [[Rcpp::export]] double log_ml_dummy( const int& p, const arma::vec& hyper, + const arma::vec& model, const arma::mat& Y, const arma::mat& X, Rcpp::List prior @@ -192,11 +260,13 @@ double log_ml_dummy( int N = Y.n_cols; int n_hyper = hyper.n_elem; - List extended = extend_dummy(p, hyper, Y, X); - mat Yplus = as(extended["Y"]); - mat Xplus = as(extended["X"]); + List extended = extend_dummy(p, model, hyper, Y, X); + mat Yplus = as(extended["Yplus"]); + mat Xplus = as(extended["Xplus"]); + mat Ystar = as(extended["Ystar"]); + mat Xstar = as(extended["Xstar"]); - prior = update_prior(p, hyper, prior); + prior = update_prior(p, hyper, model, prior); mat prior_B = as(prior["B"]); mat prior_V = as(prior["V"]); mat prior_S = as(prior["S"]); @@ -204,9 +274,6 @@ double log_ml_dummy( mat inv_V = diagmat(1 / prior_V.diag()); - mat Ystar = Yplus.rows(0, N); - mat Xstar = Xplus.rows(0, N); - double log_ml_extended = log_ml(p, prior_B, prior_V, prior_S, prior_nu, inv_V, Yplus, Xplus); double log_ml_dummy = log_ml(p, prior_B, prior_V, prior_S, prior_nu, inv_V, Ystar, Xstar); @@ -214,66 +281,22 @@ double log_ml_dummy( } -// log posterior of hyper-parameters +// log posterior of hyper-parameters (up to a constant) // [[Rcpp:interface(cpp)]] // [[Rcpp::export]] double log_posterior_hyper( const int& p, const arma::vec& hyper, + const arma::vec& model, const arma::mat& Y, const arma::mat& X, const Rcpp::List& prior ) { - double log_prior = log_prior_hyper(hyper, prior); - double log_ml = log_ml_dummy(p, hyper, Y, X, prior); + double log_prior = log_prior_hyper(hyper, model, prior); + double log_ml = log_ml_dummy(p, hyper, model, Y, X, prior); return log_prior + log_ml; } -// sample hyper-parameters with Metropolis-Hastings -// [[Rcpp:interface(cpp)]] -// [[Rcpp::export]] -arma::mat sample_hyper( - const int& S, - const int& p, - const double& c, - const arma::mat& Y, - const arma::mat& X, - const Rcpp::List& prior -) { - - int N = Y.n_cols; - int n_hyper = N + 3; - - vec hyper = as(prior["map"]); - - double logp = log_posterior_hyper(p, hyper, Y, X, prior); - - mat W = c * as(prior["W"]); - - mat posterior_hyper(n_hyper, S); - - int success = 0; - for (int s = 0; s < S; s++) { - - // sample from proposal distribution - vec hyperp = mvnrnd(hyper, c * W); - double logpp = log_posterior_hyper(p, hyperp, Y, X, prior); - double r = exp(logpp - logp); - - if (randu() < r) { - hyper = hyperp; - logp = logpp; - success++; - } - - posterior_hyper.col(s) = hyper; - } - - cout << "Acceptance rate: " << (double) success / S << endl; - - return posterior_hyper; -} - diff --git a/src/sample_hyper.h b/src/sample_hyper.h index 9648922..2657272 100644 --- a/src/sample_hyper.h +++ b/src/sample_hyper.h @@ -7,28 +7,35 @@ using namespace Rcpp; using namespace arma; +Rcpp::List mn_prior( + const int& p, + const double& lambda, + const arma::vec& psi +); + + Rcpp::List update_prior( - const int& p, - const arma::vec& hyper, - const Rcpp::List& prior + const int& p, + const arma::vec& hyper, + const Rcpp::List& prior ); Rcpp::List extend_dummy( - const int& p, - const arma::vec& hyper, - const arma::mat& Y, - const arma::mat& X + const int& p, + const arma::vec& hyper, + const arma::mat& Y, + const arma::mat& X ); arma::mat sample_hyper( - const int& S, - const int& p, - const double& c, - const arma::mat& Y, - const arma::mat& X, - const Rcpp::List& prior + const int& S, + const int& p, + const double& c, + const arma::mat& Y, + const arma::mat& X, + const Rcpp::List& prior );