diff --git a/DESCRIPTION b/DESCRIPTION index 1c5edaa..936fa64 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -16,7 +16,8 @@ LinkingTo: bsvars, Rcpp, RcppArmadillo, - RcppProgress + RcppProgress, + ramcmc Depends: R (>= 2.10), bsvars diff --git a/NAMESPACE b/NAMESPACE index 31d5bd7..96ea8ba 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -17,5 +17,4 @@ importFrom(bsvars,compute_structural_shocks) importFrom(bsvars,compute_variance_decompositions) importFrom(bsvars,estimate) importFrom(bsvars,forecast) -importFrom(stats,lm) useDynLib(bsvarSIGNs, .registration = TRUE) diff --git a/R/RcppExports.R b/R/RcppExports.R index 12cb357..105373a 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -9,20 +9,8 @@ riwish_cpp <- function(S, nu) { .Call(`_bsvarSIGNs_riwish_cpp`, S, nu) } -niw_cpp <- function(Y, X, prior) { - .Call(`_bsvarSIGNs_niw_cpp`, Y, X, prior) -} - -mn_prior <- function(p, lambda, psi) { - .Call(`_bsvarSIGNs_mn_prior`, p, lambda, psi) -} - -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) +niw_cpp <- function(Y, X, prior_B, prior_V, prior_S, prior_nu) { + .Call(`_bsvarSIGNs_niw_cpp`, Y, X, prior_B, prior_V, prior_S, prior_nu) } log_dgamma <- function(x, k, theta) { @@ -41,16 +29,28 @@ 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 <- function(b, Omega, Psi, d, Y, X) { + .Call(`_bsvarSIGNs_log_ml`, b, Omega, Psi, d, Y, X) +} + +log_ml_dummy <- function(hyper, model, Y, X, prior) { + .Call(`_bsvarSIGNs_log_ml_dummy`, hyper, model, Y, X, prior) +} + +log_posterior_hyper <- function(hyper, model, Y, X, prior) { + .Call(`_bsvarSIGNs_log_posterior_hyper`, hyper, model, Y, X, prior) +} + +extend_hyper <- function(init, model, hypers) { + .Call(`_bsvarSIGNs_extend_hyper`, init, model, hypers) } -log_ml_dummy <- function(p, hyper, model, Y, X, prior) { - .Call(`_bsvarSIGNs_log_ml_dummy`, p, hyper, model, Y, X, prior) +narrow_hyper <- function(model, hypers) { + .Call(`_bsvarSIGNs_narrow_hyper`, model, hypers) } -log_posterior_hyper <- function(p, hyper, model, Y, X, prior) { - .Call(`_bsvarSIGNs_log_posterior_hyper`, p, hyper, model, Y, X, prior) +sample_hyper <- function(S, start, init, model, Y, X, W, prior) { + .Call(`_bsvarSIGNs_sample_hyper`, S, start, init, model, Y, X, W, prior) } # Register entry points for exported C++ functions diff --git a/R/bsvarSIGNs-package.R b/R/bsvarSIGNs-package.R index 63fcf8e..4187203 100644 --- a/R/bsvarSIGNs-package.R +++ b/R/bsvarSIGNs-package.R @@ -16,7 +16,6 @@ #' @import bsvars #' @import RcppArmadillo #' @import RcppProgress -#' @importFrom stats lm #' @note This package is currently in active development. We give no #' warranty that anything here works. #' @author Xiaolei Wang \email{adamwang15@gmail.com} Tomasz Woźniak \email{wozniak.tom@pm.me} diff --git a/R/estimate.BSVARSIGN.R b/R/estimate.BSVARSIGN.R index 148800f..ab19e9d 100644 --- a/R/estimate.BSVARSIGN.R +++ b/R/estimate.BSVARSIGN.R @@ -92,7 +92,7 @@ #' posterior = estimate(specification, S = 10) #' #' @export -estimate.BSVARSIGN <- function(specification, S, thin = 1, show_progress = TRUE) { +estimate.BSVARSIGN = function(specification, S, thin = 1, show_progress = TRUE) { # get the inputs to estimation # prior = specification$last_draw$prior$get_prior() @@ -109,6 +109,7 @@ estimate.BSVARSIGN <- function(specification, S, thin = 1, show_progress = TRUE) identification$sign_irf, identification$sign_narrative, identification$sign_B, identification$zero_irf, 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/specify_bsvarSIGN.R b/R/specify_bsvarSIGN.R index ddc2362..933f4b6 100644 --- a/R/specify_bsvarSIGN.R +++ b/R/specify_bsvarSIGN.R @@ -1,5 +1,5 @@ -# construction Z_j matrices +# construct Z_j matrices get_Z = function(zero_irf) { if (sum(zero_irf) == 0) { return(NULL) @@ -85,28 +85,43 @@ verify_all = function(N, sign_irf, sign_narrative, sign_B) { # Minnesota prior of Normal-Inverse-Wishart form -niw_prior = function(Y, - p, - non_stationary, - lambda = 0.2) { - T = nrow(Y) - N = ncol(Y) - K = 1 + N * p - - B = matrix(0, K, N) - B[1:N, 1:N] = diag(non_stationary) - - sigma2 = sapply(1:N, \(i) summary(lm(Y[2:T, i] ~ Y[1:(T - 1), i]))$sigma^2) - V = matrix(0, K, K) - V[K, K] = 1e+6 - V[1:(K - 1), 1:(K - 1)] = diag(lambda^2 * kronecker((1:p)^-2, sigma2^-1)) - - S = diag(sigma2) - nu = N + 2 - - list(B = B, V = V, S = S, nu = nu) +# niw_prior = function(Y, +# p, +# non_stationary, +# lambda = 0.2) { +# T = nrow(Y) +# N = ncol(Y) +# K = 1 + N * p +# +# B = matrix(0, K, N) +# B[1:N, 1:N] = diag(non_stationary) +# +# sigma2 = sapply(1:N, \(i) summary(lm(Y[2:T, i] ~ Y[1:(T - 1), i]))$sigma^2) +# V = matrix(0, K, K) +# V[K, K] = 1e+6 +# V[1:(K - 1), 1:(K - 1)] = diag(lambda^2 * kronecker((1:p)^-2, sigma2^-1)) +# +# S = diag(sigma2) +# nu = N + 2 +# +# list(B = B, V = V, S = S, nu = nu) +# } + +gamma_scale = function(mode, variance) { + (2 * mode * variance) / (mode^2 + sqrt(mode^4 + 4 * variance * mode^2)) } +gamma_shape = function(mode, variance) { + (mode^2 + sqrt(mode^4 + 4 * (mode^2) * variance) + 2 * variance) / (2 * variance) +} + +igamma_scale = function(mode, variance) { + (0.5 * mode * (sqrt(variance * (variance + 24 * mode^2)) - 5 * variance)) / (mode^2 - variance) +} + +igamma_shape = function(mode, variance) { + 0.5 * (sqrt(variance * (variance + 24 * mode^2)) - 2 * mode^2 - 3 * variance ) / (mode^2 - variance) +} #' R6 Class Representing PriorBSVAR #' @@ -116,7 +131,7 @@ niw_prior = function(Y, #' @examples #' # a prior for 3-variable example with one lag #' data(oil) -#' prior = specify_prior_bsvarSIGN$new(oil, N = 3, p = 1) +#' prior = specify_prior_bsvarSIGN$new(oil, p = 1) #' prior$B # show autoregressive prior mean #' #' @export @@ -124,19 +139,81 @@ specify_prior_bsvarSIGN = R6::R6Class( "PriorBSVARSIGN", public = list( - #' @field hyper a \code{(N+3)xS} matrix of posterior draws of hyperparameters. + #' @field hyper a \code{(N+3)xS} matrix of hyper-parameters \eqn{\mu, \delta, \lambda, \psi}. hyper = matrix(), - #' @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(), + #' @field B a \code{KxN} normal prior mean matrix for the autoregressive + #' parameters. + B = 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(), + + #' @field nu a positive scalar greater than \code{N+1} - the shape of the + #' inverted-Wishart prior for error terms covariance matrix. + nu = NA, + + #' @field data an \code{TxN} matrix of observations. + data = matrix(), + + #' @field Y an \code{TxN} matrix of dependent variables. + Y = matrix(), + + #' @field X an \code{TxK} 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. + 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. + Xsur = matrix(), + + #' @field mu.scale a positive scalar - the shape of the gamma prior for \eqn{\mu}. + mu.scale = NA, + + #' @field mu.shape a positive scalar - the shape of the gamma prior for \eqn{\mu}. + mu.shape = NA, + + #' @field delta.scale a positive scalar - the shape of the gamma prior for \eqn{\delta}. + delta.scale = NA, + + #' @field delta.shape a positive scalar - the shape of the gamma prior for \eqn{\delta}. + delta.shape = NA, + + #' @field lambda.scale a positive scalar - the shape of the gamma prior for \eqn{\lambda}. + lambda.scale = NA, + + #' @field lambda.shape a positive scalar - the shape of the gamma prior for \eqn{\lambda}. + lambda.shape = NA, + + #' @field psi.scale a positive scalar - the shape of the inverted gamma prior for \eqn{\psi}. + psi.scale = NA, + + #' @field psi.shape a positive scalar - the shape of the inverted gamma prior for \eqn{\psi}. + psi.shape = NA, #' @description #' Create a new prior specification PriorBSVAR. - #' @param data the \code{TxN} data matrix. - #' @param N a positive integer - the number of dependent variables in the model. + #' @param data the \code{TxN} data matrix of observations. #' @param p a positive integer - the autoregressive lag order of the SVAR model. - #' @param d a positive integer - the number of \code{exogenous} variables in the model. + #' @param exogenous a \code{Txd} matrix of exogenous variables. #' @param stationary an \code{N} logical vector - its element set to \code{FALSE} sets #' the prior mean for the autoregressive parameters of the \code{N}th equation to the white noise process, #' otherwise to random walk. @@ -144,27 +221,89 @@ specify_prior_bsvarSIGN = R6::R6Class( #' @examples #' # a prior for 3-variable example with one lag and stationary data #' data(oil) - #' prior = specify_prior_bsvarSIGN$new(oil, N = 3, p = 1, stationary = rep(TRUE, 3)) + #' prior = specify_prior_bsvarSIGN$new(oil, p = 1) #' prior$B # show autoregressive prior mean #' - initialize = function(data, N, p, d = 0, stationary = rep(FALSE, N)){ - stopifnot("Argument N must be a positive integer number." = N > 0 & N %% 1 == 0) + initialize = function(data, p, exogenous = NULL, stationary = rep(FALSE, dim(data)[2])) { + stopifnot("Argument p must be a positive integer number." = p > 0 & p %% 1 == 0) - 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)) - Y = data - T = nrow(Y) - lambda = 0.2 - psi = sapply(1:N, \(i) summary(stats::lm(Y[2:T, i] ~ Y[1:(T - 1), i]))$sigma^2) + data_m = bsvars::specify_data_matrices$new(data, p, exogenous) + Y = t(data_m$Y) + X = t(data_m$X) + N = ncol(Y) + + stopifnot("Argument stationary must be a logical vector of length equal to the number of columns in data." = length(stationary) == N & is.logical(stationary)) + + d = 0 + if (!is.null(exogenous)) { + d = ncol(exogenous) + } + T = nrow(Y) + K = N * p + 1 + d + + 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)) + + s2.ols = rep(NA, N) + for (n in 1:N) { + y = as.matrix(Y[(p + 5 + 1):T, n]) + x = matrix(1, T - p - 5, 1) + for (i in 1:(p + 5)) { + x = cbind(x, Y[(p + 5 + 1):T - i, n]) + } + s2.ols[n] = sum(((diag(T - p - 5) - x %*% solve(t(x) %*% x) %*% t(x)) %*% y)^2) / + (T - p - 5) + } hyper = matrix(NA, N + 3, 1) - hyper[3, ] = lambda - hyper[4:(N + 3), ] = psi + hyper[1:3] = c(1, 1, 0.2) + hyper[4:(N + 3),] = s2.ols - self$model = c(FALSE, FALSE, FALSE, FALSE) - self$hyper = hyper + scale = gamma_scale(1, 1) + shape = gamma_shape(1, 1) + ybar = colMeans(matrix(Y[1:p,], ncol = N)) + Ysoc = diag(ybar) + Ysur = t(ybar) + Xsoc = cbind(kronecker(t(rep(1, p)), Ysoc), matrix(0, N, d + 1)) + Xsur = cbind(kronecker(t(rep(1, p)), Ysur), 1, matrix(0, 1, d)) + + # Ystar = rbind(diag(ybar), ybar) + # Xstar = Ystar + # if (p > 1) { + # for (i in 2:p) { + # Xstar = cbind(Xstar, Ystar) + # } + # } + # 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$hyper = hyper + self$B = 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$mu.scale = scale + self$mu.shape = shape + self$delta.scale = scale + self$delta.shape = shape + self$lambda.scale = gamma_scale(0.2, 0.4) + self$lambda.shape = gamma_shape(0.2, 0.4) + self$psi.scale = igamma_scale(0.02^2, 0.02^2) + self$psi.shape = igamma_shape(0.02^2, 0.02^2) }, # END initialize #' @description @@ -177,10 +316,86 @@ specify_prior_bsvarSIGN = R6::R6Class( #' get_prior = function(){ list( - model = self$model, - hyper = self$hyper + hyper = self$hyper, + B = self$B, + V = self$V, + Vp = self$Vp, + Vd = self$Vd, + S = self$S, + nu = self$nu, + Ysoc = self$Ysoc, + Xsoc = self$Xsoc, + Ysur = self$Ysur, + Xsur = self$Xsur, + mu.scale = self$mu.scale, + mu.shape = self$mu.shape, + delta.scale = self$delta.scale, + delta.shape = self$delta.shape, + lambda.scale = self$lambda.scale, + lambda.shape = self$lambda.shape, + psi.scale = self$psi.scale, + psi.shape = self$psi.shape ) - } # END get_prior + }, # END get_prior + + #' @description + #' Estimates hyper-parameters with adaptive Metropolis algorithm. + #' + #' @param mu whether to estimate the hyper-parameter in the + #' sum-of-coefficients dummy prior. + #' @param delta whether to estimate the hyper-parameter in the + #' single-unit-root dummy prior. + #' @param lambda whether to estimate the hyper-parameter of the + #' shrinkage in the Minnesota prior. + #' @param psi whether to estimate the hyper-parameter of the + #' variances in the Minnesota prior. + #' @param S number of MCMC draws. + #' @param burn_in number of burn-in draws. + #' + #' @examples + #' # a prior for 3-variable example with four lags + #' data(oil) + #' prior = specify_prior_bsvarSIGN$new(oil, p = 1) + #' prior$estimate_hyper(S = 5) + #' + estimate_hyper = function( + S = 10000, burn_in = S / 2, + mu = FALSE, delta = FALSE, lambda = TRUE, psi = FALSE + ) { + + model = c(mu, delta, lambda, psi) + + if (all(!model)) { + stop("At least one of the hyper-parameters must be estimated.") + } + + hyper = matrix(self$hyper[, ncol(self$hyper)]) + init = narrow_hyper(model, hyper) + prior = self$get_prior() + + result = stats::optim( + init, + \(x) -log_posterior_hyper(extend_hyper(hyper, model, matrix(x)), + model, self$Y, self$X, prior), + method = 'L-BFGS-B', + lower = rep(0, length(init)), + upper = init * 100, + hessian = TRUE + ) + + mode = extend_hyper(hyper, model, matrix(result$par)) + variance = result$hessian + + if (length(init) == 1){ + variance = 1 / variance + } else { + e = eigen(variance) + 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 = self$hyper[, -(1:burn_in)] + } # END estimate_hyper ) # END public ) # END specify_prior_bsvarSIGN @@ -509,7 +724,7 @@ specify_bsvarSIGN = R6::R6Class( sign_B, zero_irf, max_tries) - self$prior = specify_prior_bsvarSIGN$new(data, N, p, 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 diff --git a/R/utils.R b/R/utils.R index a084d9e..c3531dc 100644 --- a/R/utils.R +++ b/R/utils.R @@ -17,9 +17,7 @@ importance_sampling <- function(posterior) { posterior$posterior$Q = posterior$posterior$Q[, , indices] posterior$posterior$Sigma = posterior$posterior$Sigma[, , indices] posterior$posterior$Theta0 = posterior$posterior$Theta0[, , indices] - posterior$posterior$ess = sum(w)^2/sum(w^2) - cat("The effective sample size is:", round(posterior$posterior$ess, 0), "\n") return(posterior) } diff --git a/README.Rmd b/README.Rmd index a055c9b..85c0199 100644 --- a/README.Rmd +++ b/README.Rmd @@ -13,6 +13,10 @@ knitr::opts_chunk$set( ) ``` + +[![R-CMD-check](https://github.com/bsvars/bsvarSIGNs/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/bsvars/bsvarSIGNs/actions/workflows/R-CMD-check.yaml) + + # bsvarSIGNs Developing an R package for Bayesian Structural VARs identified by zero, sign and narrative restrictions. @@ -72,7 +76,3 @@ irf = compute_impulse_responses(posterior, horizon = 40) plot(irf, probability = 0.68) ``` - - -[![R-CMD-check](https://github.com/bsvars/bsvarSIGNs/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/bsvars/bsvarSIGNs/actions/workflows/R-CMD-check.yaml) - \ No newline at end of file diff --git a/README.md b/README.md index e544e12..07252d4 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,9 @@ + + +[![R-CMD-check](https://github.com/bsvars/bsvarSIGNs/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/bsvars/bsvarSIGNs/actions/workflows/R-CMD-check.yaml) + # bsvarSIGNs @@ -60,8 +64,3 @@ posterior = estimate(specification, S = 100) irf = compute_impulse_responses(posterior, horizon = 40) plot(irf, probability = 0.68) ``` - - - -[![R-CMD-check](https://github.com/bsvars/bsvarSIGNs/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/bsvars/bsvarSIGNs/actions/workflows/R-CMD-check.yaml) - diff --git a/inst/include/bsvarSIGNs_RcppExports.h b/inst/include/bsvarSIGNs_RcppExports.h index ab1d075..43061ec 100644 --- a/inst/include/bsvarSIGNs_RcppExports.h +++ b/inst/include/bsvarSIGNs_RcppExports.h @@ -319,11 +319,11 @@ namespace bsvarSIGNs { return Rcpp::as(rcpp_result_gen); } - inline Rcpp::List sample_Q(const int& p, const arma::mat& Y, const arma::mat& X, arma::mat& B, arma::mat& h_invp, arma::mat& chol_Sigma, const Rcpp::List& prior, const arma::field& VB, const arma::cube& sign_irf, const arma::mat& sign_narrative, const arma::mat& sign_B, const arma::field& Z, const int& max_tries) { + inline arma::field sample_Q(const int& p, const arma::mat& Y, const arma::mat& X, arma::mat& B, arma::mat& h_invp, arma::mat& chol_Sigma, const Rcpp::List& prior, const arma::field& VB, const arma::cube& sign_irf, const arma::mat& sign_narrative, const arma::mat& sign_B, const arma::field& Z, const int& max_tries) { typedef SEXP(*Ptr_sample_Q)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP); static Ptr_sample_Q p_sample_Q = NULL; if (p_sample_Q == NULL) { - validateSignature("Rcpp::List(*sample_Q)(const int&,const arma::mat&,const arma::mat&,arma::mat&,arma::mat&,arma::mat&,const Rcpp::List&,const arma::field&,const arma::cube&,const arma::mat&,const arma::mat&,const arma::field&,const int&)"); + validateSignature("arma::field(*sample_Q)(const int&,const arma::mat&,const arma::mat&,arma::mat&,arma::mat&,arma::mat&,const Rcpp::List&,const arma::field&,const arma::cube&,const arma::mat&,const arma::mat&,const arma::field&,const int&)"); p_sample_Q = (Ptr_sample_Q)R_GetCCallable("bsvarSIGNs", "_bsvarSIGNs_sample_Q"); } RObject rcpp_result_gen; @@ -337,7 +337,7 @@ namespace bsvarSIGNs { throw Rcpp::LongjumpException(rcpp_result_gen); if (rcpp_result_gen.inherits("try-error")) throw Rcpp::exception(Rcpp::as(rcpp_result_gen).c_str()); - return Rcpp::as(rcpp_result_gen); + return Rcpp::as >(rcpp_result_gen); } inline arma::mat qr_sign_cpp(const arma::mat& A) { diff --git a/inst/tinytest/test_estimate.R b/inst/tinytest/test_estimate.R index 78dc9a9..b25a88d 100644 --- a/inst/tinytest/test_estimate.R +++ b/inst/tinytest/test_estimate.R @@ -45,8 +45,8 @@ expect_identical(post1$posterior$B[, , 1], expect_identical(post2$posterior$B[, , 1], post3$posterior$B[, , 1],) -expect_identical(post1$last_draw$B[, , 1], - post2$last_draw$B[, , 1],) - -expect_identical(post2$last_draw$B[, , 1], - post3$last_draw$B[, , 1],) +# expect_identical(post1$last_draw$B[, , 1], +# post2$last_draw$B[, , 1],) +# +# expect_identical(post2$last_draw$B[, , 1], +# post3$last_draw$B[, , 1],) diff --git a/inst/tinytest/test_specify.R b/inst/tinytest/test_specify.R index 64525a8..d2b740f 100644 --- a/inst/tinytest/test_specify.R +++ b/inst/tinytest/test_specify.R @@ -14,14 +14,14 @@ expect_identical(dim(spec$data_matrices$Y)[1], 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$B)[1], dim(spec$data_matrices$X)[2]) expect_identical(class(spec$prior$get_prior()), "list") -# expect_true(det(spec$prior$V) > 0) +expect_true(det(spec$prior$V) > 0) -# expect_true(det(spec$prior$S) > 0) +expect_true(det(spec$prior$S) > 0) expect_identical(class(spec$identification)[1], "IdentificationBSVARSIGN") @@ -57,14 +57,14 @@ expect_identical(dim(spec$data_matrices$Y)[1], 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$B)[1], dim(spec$data_matrices$X)[2]) expect_identical(class(spec$prior$get_prior()), "list") -# expect_true(det(spec$prior$V) > 0) +expect_true(det(spec$prior$V) > 0) -# expect_true(det(spec$prior$S) > 0) +expect_true(det(spec$prior$S) > 0) expect_identical(class(spec$identification)[1], "IdentificationBSVARSIGN") diff --git a/inst/varia/nicetry.R b/inst/varia/nicetry.R new file mode 100644 index 0000000..6b428e9 --- /dev/null +++ b/inst/varia/nicetry.R @@ -0,0 +1,24 @@ +library(R.matlab) +data = readMat('../data/DataSW.mat') +name = sapply(1:7, \(i) data$ShortDescr[[i]][[1]]) +data = data$y +colnames(data) = name + +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) + +end_time = Sys.time() +end_time - start_time + +plot.ts(t(spec$prior$hyper)) + +# hist(spec$prior$hyper[1, ], breaks = 200) diff --git a/man/specify_prior_bsvarSIGN.Rd b/man/specify_prior_bsvarSIGN.Rd index 9a82b6d..4d32ae4 100644 --- a/man/specify_prior_bsvarSIGN.Rd +++ b/man/specify_prior_bsvarSIGN.Rd @@ -9,7 +9,7 @@ The class PriorBSVARSIGN presents a prior specification for the homoskedastic bs \examples{ # a prior for 3-variable example with one lag data(oil) -prior = specify_prior_bsvarSIGN$new(oil, N = 3, p = 1) +prior = specify_prior_bsvarSIGN$new(oil, p = 1) prior$B # show autoregressive prior mean @@ -19,7 +19,7 @@ prior$B # show autoregressive prior mean # a prior for 3-variable example with one lag and stationary data data(oil) -prior = specify_prior_bsvarSIGN$new(oil, N = 3, p = 1, stationary = rep(TRUE, 3)) +prior = specify_prior_bsvarSIGN$new(oil, p = 1) prior$B # show autoregressive prior mean @@ -31,14 +31,67 @@ prior$B # show autoregressive prior mean prior = specify_prior_bsvar$new(N = 3, p = 4) prior$get_prior() # show the prior as list + +## ------------------------------------------------ +## Method `specify_prior_bsvarSIGN$estimate_hyper` +## ------------------------------------------------ + +# a prior for 3-variable example with four lags +data(oil) +prior = specify_prior_bsvarSIGN$new(oil, p = 1) +prior$estimate_hyper(S = 5) + } \section{Public fields}{ \if{html}{\out{
}} \describe{ -\item{\code{hyper}}{a \code{(N+3)xS} matrix of posterior draws of hyperparameters.} +\item{\code{hyper}}{a \code{(N+3)xS} matrix of hyper-parameters \eqn{\mu, \delta, \lambda, \psi}.} + +\item{\code{B}}{a \code{KxN} normal prior mean matrix for the autoregressive +parameters.} + +\item{\code{V}}{a \code{KxK} matrix determining the normal prior column-specific +covariance for the autoregressive parameters.} + +\item{\code{Vp}}{a \code{px1} vector of lag shrinkage level.} + +\item{\code{Vd}}{a \code{(d+1)x1} vector of (large) variances for constants.} + +\item{\code{S}}{an \code{NxN} matrix determining the inverted-Wishart prior scale +of error terms covariance matrix.} + +\item{\code{nu}}{a positive scalar greater than \code{N+1} - the shape of the +inverted-Wishart prior for error terms covariance matrix.} + +\item{\code{data}}{an \code{TxN} matrix of observations.} + +\item{\code{Y}}{an \code{TxN} matrix of dependent variables.} + +\item{\code{X}}{an \code{TxK} matrix of independent variables.} + +\item{\code{Ysoc}}{an \code{NxN} matrix with the sum-of-coefficients dummy observations.} + +\item{\code{Xsoc}}{an \code{NxK} matrix with the sum-of-coefficients dummy observations.} + +\item{\code{Ysur}}{an \code{NxN} matrix with the single-unit-root dummy observations.} + +\item{\code{Xsur}}{an \code{NxK} matrix with the single-unit-root dummy observations.} + +\item{\code{mu.scale}}{a positive scalar - the shape of the gamma prior for \eqn{\mu}.} + +\item{\code{mu.shape}}{a positive scalar - the shape of the gamma prior for \eqn{\mu}.} + +\item{\code{delta.scale}}{a positive scalar - the shape of the gamma prior for \eqn{\delta}.} + +\item{\code{delta.shape}}{a positive scalar - the shape of the gamma prior for \eqn{\delta}.} + +\item{\code{lambda.scale}}{a positive scalar - the shape of the gamma prior for \eqn{\lambda}.} + +\item{\code{lambda.shape}}{a positive scalar - the shape of the gamma prior for \eqn{\lambda}.} + +\item{\code{psi.scale}}{a positive scalar - the shape of the inverted gamma prior for \eqn{\psi}.} -\item{\code{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.} +\item{\code{psi.shape}}{a positive scalar - the shape of the inverted gamma prior for \eqn{\psi}.} } \if{html}{\out{
}} } @@ -47,6 +100,7 @@ model[1] = dummy soc, model[2] = dummy sur, model[3] = mn lambda, model[4] = mn \itemize{ \item \href{#method-PriorBSVARSIGN-new}{\code{specify_prior_bsvarSIGN$new()}} \item \href{#method-PriorBSVARSIGN-get_prior}{\code{specify_prior_bsvarSIGN$get_prior()}} +\item \href{#method-PriorBSVARSIGN-estimate_hyper}{\code{specify_prior_bsvarSIGN$estimate_hyper()}} \item \href{#method-PriorBSVARSIGN-clone}{\code{specify_prior_bsvarSIGN$clone()}} } } @@ -56,19 +110,22 @@ model[1] = dummy soc, model[2] = dummy sur, model[3] = mn lambda, model[4] = mn \subsection{Method \code{new()}}{ Create a new prior specification PriorBSVAR. \subsection{Usage}{ -\if{html}{\out{
}}\preformatted{specify_prior_bsvarSIGN$new(data, N, p, d = 0, stationary = rep(FALSE, N))}\if{html}{\out{
}} +\if{html}{\out{
}}\preformatted{specify_prior_bsvarSIGN$new( + data, + p, + exogenous = NULL, + stationary = rep(FALSE, dim(data)[2]) +)}\if{html}{\out{
}} } \subsection{Arguments}{ \if{html}{\out{
}} \describe{ -\item{\code{data}}{the \code{TxN} data matrix.} - -\item{\code{N}}{a positive integer - the number of dependent variables in the model.} +\item{\code{data}}{the \code{TxN} data matrix of observations.} \item{\code{p}}{a positive integer - the autoregressive lag order of the SVAR model.} -\item{\code{d}}{a positive integer - the number of \code{exogenous} variables in the model.} +\item{\code{exogenous}}{a \code{Txd} matrix of exogenous variables.} \item{\code{stationary}}{an \code{N} logical vector - its element set to \code{FALSE} sets the prior mean for the autoregressive parameters of the \code{N}th equation to the white noise process, @@ -83,7 +140,7 @@ A new prior specification PriorBSVARSIGN. \if{html}{\out{
}} \preformatted{# a prior for 3-variable example with one lag and stationary data data(oil) -prior = specify_prior_bsvarSIGN$new(oil, N = 3, p = 1, stationary = rep(TRUE, 3)) +prior = specify_prior_bsvarSIGN$new(oil, p = 1) prior$B # show autoregressive prior mean } @@ -112,6 +169,56 @@ prior$get_prior() # show the prior as list } +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-PriorBSVARSIGN-estimate_hyper}{}}} +\subsection{Method \code{estimate_hyper()}}{ +Estimates hyper-parameters with adaptive Metropolis algorithm. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{specify_prior_bsvarSIGN$estimate_hyper( + S = 10000, + burn_in = S/2, + mu = FALSE, + delta = FALSE, + lambda = TRUE, + psi = FALSE +)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{S}}{number of MCMC draws.} + +\item{\code{burn_in}}{number of burn-in draws.} + +\item{\code{mu}}{whether to estimate the hyper-parameter in the +sum-of-coefficients dummy prior.} + +\item{\code{delta}}{whether to estimate the hyper-parameter in the +single-unit-root dummy prior.} + +\item{\code{lambda}}{whether to estimate the hyper-parameter of the +shrinkage in the Minnesota prior.} + +\item{\code{psi}}{whether to estimate the hyper-parameter of the +variances in the Minnesota prior.} +} +\if{html}{\out{
}} +} +\subsection{Examples}{ +\if{html}{\out{
}} +\preformatted{# a prior for 3-variable example with four lags +data(oil) +prior = specify_prior_bsvarSIGN$new(oil, p = 1) +prior$estimate_hyper(S = 5) + +} +\if{html}{\out{
}} + +} + } \if{html}{\out{
}} \if{html}{\out{}} diff --git a/src/Makevars b/src/Makevars index 22c7566..62a255e 100644 --- a/src/Makevars +++ b/src/Makevars @@ -17,3 +17,8 @@ PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) + +#PKG_CXXFLAGS += -Xclang -fopenmp +#PKG_FFLAGS += $(SHLIB_OPENMP_FFLAGS) +#PKG_CFLAGS += $(SHLIB_OPENMP_CFLAGS) +#PKG_LIBS += -lomp \ No newline at end of file diff --git a/src/Makevars.win b/src/Makevars.win index 22c7566..48f15f8 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -16,4 +16,4 @@ #CXX_STD = CXX11 PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) -PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) +PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) \ No newline at end of file diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 1dd9525..aef846b 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -520,15 +520,18 @@ BEGIN_RCPP END_RCPP } // niw_cpp -Rcpp::List niw_cpp(const arma::mat& Y, const arma::mat& X, const Rcpp::List prior); -RcppExport SEXP _bsvarSIGNs_niw_cpp(SEXP YSEXP, SEXP XSEXP, SEXP priorSEXP) { +arma::field niw_cpp(const arma::mat& Y, const arma::mat& X, const arma::mat& prior_B, const arma::mat& prior_V, const arma::mat& prior_S, const int& prior_nu); +RcppExport SEXP _bsvarSIGNs_niw_cpp(SEXP YSEXP, SEXP XSEXP, SEXP prior_BSEXP, SEXP prior_VSEXP, SEXP prior_SSEXP, SEXP prior_nuSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; 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(niw_cpp(Y, X, prior)); + Rcpp::traits::input_parameter< const arma::mat& >::type prior_B(prior_BSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type prior_V(prior_VSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type prior_S(prior_SSEXP); + Rcpp::traits::input_parameter< const int& >::type prior_nu(prior_nuSEXP); + rcpp_result_gen = Rcpp::wrap(niw_cpp(Y, X, prior_B, prior_V, prior_S, prior_nu)); return rcpp_result_gen; END_RCPP } @@ -569,7 +572,7 @@ RcppExport SEXP _bsvarSIGNs_match_sign_irf(SEXP QSEXP, SEXP sign_irfSEXP, SEXP i return rcpp_result_gen; } // sample_Q -Rcpp::List sample_Q(const int& p, const arma::mat& Y, const arma::mat& X, arma::mat& B, arma::mat& h_invp, arma::mat& chol_Sigma, const Rcpp::List& prior, const arma::field& VB, const arma::cube& sign_irf, const arma::mat& sign_narrative, const arma::mat& sign_B, const arma::field& Z, const int& max_tries); +arma::field sample_Q(const int& p, const arma::mat& Y, const arma::mat& X, arma::mat& B, arma::mat& h_invp, arma::mat& chol_Sigma, const Rcpp::List& prior, const arma::field& VB, const arma::cube& sign_irf, const arma::mat& sign_narrative, const arma::mat& sign_B, const arma::field& Z, const int& max_tries); static SEXP _bsvarSIGNs_sample_Q_try(SEXP pSEXP, SEXP YSEXP, SEXP XSEXP, SEXP BSEXP, SEXP h_invpSEXP, SEXP chol_SigmaSEXP, SEXP priorSEXP, SEXP VBSEXP, SEXP sign_irfSEXP, SEXP sign_narrativeSEXP, SEXP sign_BSEXP, SEXP ZSEXP, SEXP max_triesSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; @@ -614,48 +617,6 @@ 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 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, model, prior)); - return rcpp_result_gen; -END_RCPP -} -// extend_dummy -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, model, Y, X)); - return rcpp_result_gen; -END_RCPP -} // log_dgamma double log_dgamma(const double& x, const double& k, const double& theta); RcppExport SEXP _bsvarSIGNs_log_dgamma(SEXP xSEXP, SEXP kSEXP, SEXP thetaSEXP) { @@ -708,52 +669,91 @@ BEGIN_RCPP 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) { +double log_ml(const arma::mat& b, const arma::mat& Omega, const arma::mat& Psi, const int& d, const arma::mat& Y, const arma::mat& X); +RcppExport SEXP _bsvarSIGNs_log_ml(SEXP bSEXP, SEXP OmegaSEXP, SEXP PsiSEXP, SEXP dSEXP, 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)); + rcpp_result_gen = Rcpp::wrap(log_ml(b, Omega, Psi, d, Y, X)); return rcpp_result_gen; END_RCPP } // log_ml_dummy -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) { +double log_ml_dummy(const arma::vec& hyper, const arma::vec& model, const arma::mat& Y, const arma::mat& X, const Rcpp::List& prior); +RcppExport SEXP _bsvarSIGNs_log_ml_dummy(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, model, Y, X, prior)); + Rcpp::traits::input_parameter< const Rcpp::List& >::type prior(priorSEXP); + rcpp_result_gen = Rcpp::wrap(log_ml_dummy(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::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) { +double log_posterior_hyper(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 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, model, Y, X, prior)); + rcpp_result_gen = Rcpp::wrap(log_posterior_hyper(hyper, model, Y, X, prior)); + return rcpp_result_gen; +END_RCPP +} +// extend_hyper +arma::mat extend_hyper(const arma::vec& init, const arma::vec& model, const arma::mat& hypers); +RcppExport SEXP _bsvarSIGNs_extend_hyper(SEXP initSEXP, SEXP modelSEXP, SEXP hypersSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::vec& >::type init(initSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type model(modelSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type hypers(hypersSEXP); + rcpp_result_gen = Rcpp::wrap(extend_hyper(init, model, hypers)); + return rcpp_result_gen; +END_RCPP +} +// narrow_hyper +arma::mat narrow_hyper(const arma::vec& model, arma::mat hypers); +RcppExport SEXP _bsvarSIGNs_narrow_hyper(SEXP modelSEXP, SEXP hypersSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::vec& >::type model(modelSEXP); + Rcpp::traits::input_parameter< arma::mat >::type hypers(hypersSEXP); + rcpp_result_gen = Rcpp::wrap(narrow_hyper(model, hypers)); + return rcpp_result_gen; +END_RCPP +} +// sample_hyper +arma::mat sample_hyper(const int& S, const int& start, const arma::vec& init, const arma::vec& model, const arma::mat& Y, const arma::mat& X, const arma::mat& W, const Rcpp::List& prior); +RcppExport SEXP _bsvarSIGNs_sample_hyper(SEXP SSEXP, SEXP startSEXP, SEXP initSEXP, SEXP modelSEXP, SEXP YSEXP, SEXP XSEXP, SEXP WSEXP, 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 start(startSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type init(initSEXP); + 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 arma::mat& >::type W(WSEXP); + Rcpp::traits::input_parameter< const Rcpp::List& >::type prior(priorSEXP); + rcpp_result_gen = Rcpp::wrap(sample_hyper(S, start, init, model, Y, X, W, prior)); return rcpp_result_gen; END_RCPP } @@ -879,7 +879,7 @@ static int _bsvarSIGNs_RcppExport_validate(const char* sig) { signatures.insert("double(*weight_zero)(const arma::field&,const arma::mat&,const arma::mat&,const arma::mat&)"); signatures.insert("arma::mat(*rzeroQ)(const arma::field&,const arma::mat&)"); signatures.insert("bool(*match_sign_irf)(const arma::mat&,const arma::cube&,const arma::cube&)"); - signatures.insert("Rcpp::List(*sample_Q)(const int&,const arma::mat&,const arma::mat&,arma::mat&,arma::mat&,arma::mat&,const Rcpp::List&,const arma::field&,const arma::cube&,const arma::mat&,const arma::mat&,const arma::field&,const int&)"); + signatures.insert("arma::field(*sample_Q)(const int&,const arma::mat&,const arma::mat&,arma::mat&,arma::mat&,arma::mat&,const Rcpp::List&,const arma::field&,const arma::cube&,const arma::mat&,const arma::mat&,const arma::field&,const int&)"); signatures.insert("arma::mat(*qr_sign_cpp)(const arma::mat&)"); signatures.insert("arma::mat(*rortho_cpp)(const int&)"); signatures.insert("bool(*match_sign)(const arma::mat&,const arma::mat&)"); @@ -927,19 +927,19 @@ static const R_CallMethodDef CallEntries[] = { {"_bsvarSIGNs_rzeroQ", (DL_FUNC) &_bsvarSIGNs_rzeroQ, 2}, {"_bsvarSIGNs_rmatnorm_cpp", (DL_FUNC) &_bsvarSIGNs_rmatnorm_cpp, 3}, {"_bsvarSIGNs_riwish_cpp", (DL_FUNC) &_bsvarSIGNs_riwish_cpp, 2}, - {"_bsvarSIGNs_niw_cpp", (DL_FUNC) &_bsvarSIGNs_niw_cpp, 3}, + {"_bsvarSIGNs_niw_cpp", (DL_FUNC) &_bsvarSIGNs_niw_cpp, 6}, {"_bsvarSIGNs_match_sign_irf", (DL_FUNC) &_bsvarSIGNs_match_sign_irf, 3}, {"_bsvarSIGNs_sample_Q", (DL_FUNC) &_bsvarSIGNs_sample_Q, 13}, - {"_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, 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, 6}, - {"_bsvarSIGNs_log_posterior_hyper", (DL_FUNC) &_bsvarSIGNs_log_posterior_hyper, 6}, + {"_bsvarSIGNs_log_ml", (DL_FUNC) &_bsvarSIGNs_log_ml, 6}, + {"_bsvarSIGNs_log_ml_dummy", (DL_FUNC) &_bsvarSIGNs_log_ml_dummy, 5}, + {"_bsvarSIGNs_log_posterior_hyper", (DL_FUNC) &_bsvarSIGNs_log_posterior_hyper, 5}, + {"_bsvarSIGNs_extend_hyper", (DL_FUNC) &_bsvarSIGNs_extend_hyper, 3}, + {"_bsvarSIGNs_narrow_hyper", (DL_FUNC) &_bsvarSIGNs_narrow_hyper, 2}, + {"_bsvarSIGNs_sample_hyper", (DL_FUNC) &_bsvarSIGNs_sample_hyper, 8}, {"_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 825b856..ead6630 100644 --- a/src/bsvars_sign.cpp +++ b/src/bsvars_sign.cpp @@ -2,8 +2,8 @@ #include #include "progress.hpp" #include "Rcpp/Rmath.h" - #include +// #include #include "sample_hyper.h" #include "sample_Q.h" @@ -38,89 +38,123 @@ Rcpp::List bsvar_sign_cpp( } // Progress bar setup - vec prog_rep_points = arma::round(arma::linspace(0, S, 50)); + double num_threads = 1; + // #pragma omp parallel + // { + // num_threads = omp_get_num_threads(); + // } + vec prog_rep_points = arma::round(arma::linspace(0, S / num_threads, 50)); if (show_progress) { Rcout << "**************************************************|" << endl; - Rcout << "bsvars: Bayesian Structural Vector Autoregressions|" << endl; - Rcout << "**************************************************|" << endl; - Rcout << " Gibbs sampler for the SVAR model |" << endl; + Rcout << " bsvarSIGNs: Bayesian Structural VAR with zero, |" << endl; + Rcout << " sign and narrative restriction |" << endl; Rcout << "**************************************************|" << endl; - Rcout << " Progress of the MCMC simulation for " << S << " draws" << endl; - Rcout << " Every " << oo << "draw is saved via MCMC thinning" << endl; + // Rcout << " Gibbs sampler for the SVAR model |" << endl; + // Rcout << "**************************************************|" << endl; + Rcout << " Progress of simulation for " << S << " independent draws" << endl; + // Rcout << " Every " << oo << "draw is saved via MCMC thinning" << endl; Rcout << " Press Esc to interrupt the computations" << endl; Rcout << "**************************************************|" << endl; } Progress p(50, show_progress); - const int T = Y.n_rows; - const int N = Y.n_cols; - const int K = X.n_cols; + const int T = Y.n_rows; + const int N = Y.n_cols; + const int K = X.n_cols; - mat aux_B = as(starting_values["B"]); - mat aux_A = as(starting_values["A"]); - mat aux_hyper = as(starting_values["hyper"]); + vec posterior_w(S); + 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_shocks(N, T, S); - vec posterior_w(S); + mat hypers = as(prior["hyper"]); - vec model = as(prior["model"]); - mat sample_hyper = as(prior["hyper"]); + int S_hyper = hypers.n_cols - 1; + int prior_nu = as(prior["nu"]); + int post_nu = prior_nu + T; - 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_shocks(N, T, S); + double w, mu, delta, lambda; - 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; + vec hyper, psi; - while (s < S) { - - // Check for user interrupts - if (s % 200 == 0) checkUserInterrupt(); + mat B, Sigma, chol_Sigma, h_invp, Q, shocks; + mat prior_V, prior_S, post_B, post_V, post_S; + mat Ystar, Xstar, Yplus, Xplus; + mat prior_B = as(prior["B"]); + mat Ysoc = as(prior["Ysoc"]); + mat Xsoc = as(prior["Xsoc"]); + mat Ysur = as(prior["Ysur"]); + mat Xsur = as(prior["Xsur"]); + + field result; + + // #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 = sample_hyper.col(randi(distr_param(0, S_hyper))); - lambda = hyper(2); - psi = hyper.rows(3, N + 2); + 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(join_vert(lambda*lambda * kron(as(prior["Vp"]), 1 / psi), + as(prior["Vd"]))); + 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(post(3)); + + // sample reduced-form parameters + Sigma = iwishrnd(post_S, post_nu); + chol_Sigma = chol(Sigma, "lower"); + B = rmatnorm_cpp(post_B, post_V, Sigma); + h_invp = inv(trimatl(chol_Sigma)); // lower tri, h(Sigma) is upper tri + + result = sample_Q(lags, Y, X, B, h_invp, chol_Sigma, prior, + VB, sign_irf, sign_narrative, sign_B, Z, max_tries); + Q = result(0); + shocks = result(1); + w = as_scalar(result(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"]); + posterior_w(s) = w; + posterior_hyper.col(s) = hyper; + posterior_A.slice(s) = B.t(); + posterior_B.slice(s) = Q.t() * h_invp; + posterior_Sigma.slice(s) = Sigma; + posterior_Theta0.slice(s) = chol_Sigma * Q; + posterior_shocks.slice(s) = shocks; - Sigma = iwishrnd(post_S, post_nu); - chol_Sigma = chol(Sigma, "lower"); - B = rmatnorm_cpp(post_B, post_V, Sigma); - h_invp = inv(trimatl(chol_Sigma)); // lower triangular, h(Sigma) is upper triangular + // Check for user interrupts + if (s % 200 == 0) checkUserInterrupt(); - result = sample_Q(lags, Y, X, - B, h_invp, chol_Sigma, - prior, VB, - sign_irf, sign_narrative, sign_B, Z, - max_tries); - Q = as(result["Q"]); - w = as(result["w"]); - Epsilon = as(result["Epsilon"]); + // Increment progress bar + if (any(prog_rep_points == s)) p.increment(); + + // if (omp_get_thread_num() == 0) { + // // Check for user interrupts + // if (s % 10 == 0) checkUserInterrupt(); + // + // // Increment progress bar + // if (any(prog_rep_points == s)) p.increment(); + // } - if (w > 0) { - // Increment progress bar - if (any(prog_rep_points == s)) p.increment(); - - posterior_w(s) = w; - posterior_A.slice(s) = B.t(); - posterior_B.slice(s) = Q.t() * h_invp; - posterior_Sigma.slice(s) = Sigma; - posterior_Theta0.slice(s) = chol_Sigma * Q; - posterior_shocks.slice(s) = Epsilon; - posterior_hyper.col(s) = hyper; - - s++; - } } // END s loop return List::create( @@ -136,5 +170,3 @@ Rcpp::List bsvar_sign_cpp( ); } // END bsvar_sign_cpp - - diff --git a/src/mcmc.cpp b/src/mcmc.cpp new file mode 100644 index 0000000..decbbcc --- /dev/null +++ b/src/mcmc.cpp @@ -0,0 +1,70 @@ + +#include +#include +#include +#include "Rcpp/Rmath.h" +#include "progress.hpp" +#include "ramcmc.h" + +using namespace arma; + + +// adaptive Metropolis algorithm for strictly positive parameters +// [[Rcpp:interfaces(cpp)]] +arma::mat metropolis( + const int& T, + const int& t0, + arma::vec x, + arma::mat Sigma, + const std::function& log_target +) { + + int n = x.n_elem; + + x = log(x); + vec xbar = x; + double s = 2.38 / sqrt(n); + double d = log_target(x) + sum(x); + + double new_d, a; + vec new_x, diff; + + mat X(n, T); + X.col(0) = x; + + vec progress = round(linspace(0, T, 50)); + Progress p(50, true); + + for (int t = 1; t < T; t++) { + new_x = mvnrnd(x, s*s * Sigma); + new_d = log_target(exp(new_x)) + sum(new_x); + a = std::min(1.0, exp(new_d - d)); + + if (randu() < a) { + x = new_x; + d = new_d; + } + X.col(t) = x; + + if (t <= t0) { + diff = x - xbar; + s += pow(t, -0.6) * (a - 0.234); + xbar += diff / (t + 1); + Sigma = Sigma * t / (t + 1) + diff * diff.t() * t / (t + 1) / (t + 1); + } + + if (any(progress == t)) { + p.increment(); + } + } + + X = exp(X); + return X; +} + + + + + + + diff --git a/src/mcmc.h b/src/mcmc.h new file mode 100644 index 0000000..b76052a --- /dev/null +++ b/src/mcmc.h @@ -0,0 +1,23 @@ +#ifndef _MCMC_H_ +#define _MCMC_H_ + + +#include +#include +#include +#include "Rcpp/Rmath.h" +#include "progress.hpp" + +using namespace arma; + + +arma::mat metropolis( + const int& T, + const int& t0, + arma::vec x, + arma::mat Sigma, + const std::function& log_target +); + + +#endif // _MCMC_H_ \ No newline at end of file diff --git a/src/sample_NIW.cpp b/src/sample_NIW.cpp index 079a394..bea6cd9 100644 --- a/src/sample_NIW.cpp +++ b/src/sample_NIW.cpp @@ -46,19 +46,17 @@ arma::mat riwish_cpp ( // [[Rcpp:interface(cpp)]] // [[Rcpp::export]] -Rcpp::List niw_cpp( +arma::field niw_cpp( const arma::mat& Y, const arma::mat& X, - const Rcpp::List prior + const arma::mat& prior_B, + const arma::mat& prior_V, + const arma::mat& prior_S, + const int& prior_nu ) { const int T = Y.n_rows; - - mat prior_B = as(prior["B"]); - mat prior_V = as(prior["V"]); - mat prior_S = as(prior["S"]); - int prior_nu = as(prior["nu"]); - + // analytic solutions mat prior_V_inv = inv_sympd(prior_V); mat post_V_inv = prior_V_inv + X.t() * X; @@ -70,11 +68,12 @@ Rcpp::List niw_cpp( post_S = symmatu(post_S); int post_nu = prior_nu + T; - return List::create( - Named("B") = post_B, - Named("V") = post_V, - Named("S") = post_S, - Named("nu") = post_nu - ); + field post(4); + post(0) = post_B; + post(1) = post_V; + post(2) = post_S; + post(3) = mat(1, 1, fill::ones) * post_nu; + + return post; } diff --git a/src/sample_NIW.h b/src/sample_NIW.h index 39a194c..5e5832e 100644 --- a/src/sample_NIW.h +++ b/src/sample_NIW.h @@ -9,19 +9,24 @@ using namespace Rcpp; using namespace arma; -arma::mat rmatnorm_cpp(const arma::mat& M, - const arma::mat& U, - const arma::mat& V); +arma::mat rmatnorm_cpp( + const arma::mat& M, + const arma::mat& U, + const arma::mat& V +); arma::mat riwish_cpp ( const arma::mat& S, const double& nu ); -Rcpp::List niw_cpp( +arma::field niw_cpp( const arma::mat& Y, const arma::mat& X, - const Rcpp::List prior + const arma::mat& prior_B, + const arma::mat& prior_V, + const arma::mat& prior_S, + const int& prior_nu ); #endif // _SAMPLE_NIW_H_ diff --git a/src/sample_Q.cpp b/src/sample_Q.cpp index 65bb4b5..c82698d 100644 --- a/src/sample_Q.cpp +++ b/src/sample_Q.cpp @@ -32,7 +32,7 @@ bool match_sign_irf( // Sample rotation matrix Q and updates importance weight by reference // [[Rcpp::interfaces(cpp)]] // [[Rcpp::export]] -Rcpp::List sample_Q( +arma::field sample_Q( const int& p, const arma::mat& Y, const arma::mat& X, @@ -59,12 +59,12 @@ Rcpp::List sample_Q( bool has_zero = Z.n_elem > 0; mat Q(N, N); - mat lt_Epsilon = h_invp * (Y - X * B).t(); // lower triangular identified shocks + mat lt_shocks = h_invp * (Y - X * B).t(); // lower triangular identified shocks cube irf = ir1_cpp(B, chol_Sigma, h, p); // reduced-form irf bool success = false; - mat Epsilon; + mat shocks; while (n_tries < max_tries && !success) { if (has_zero) { Q = rzeroQ(Z, irf.slice(0)); @@ -72,13 +72,13 @@ Rcpp::List sample_Q( Q = rortho_cpp(N); } - Epsilon = Q.t() * lt_Epsilon; + shocks = Q.t() * lt_shocks; if (match_sign_irf(Q, sign_irf, irf) && match_sign(Q.t() * h_invp, sign_B)) { if (!has_narrative) { success = true; } else { - success = match_sign_narrative(Epsilon, sign_narrative, irf); + success = match_sign_narrative(shocks, sign_narrative, irf); } } @@ -97,9 +97,11 @@ Rcpp::List sample_Q( } } - return List::create( - Named("Q") = Q, - Named("w") = w, - Named("Epsilon") = Epsilon - ); -} \ No newline at end of file + field result(3); + result(0) = Q; + result(1) = shocks; + result(2) = w; + + return result; +} + diff --git a/src/sample_Q.h b/src/sample_Q.h index 651ec20..33487c9 100644 --- a/src/sample_Q.h +++ b/src/sample_Q.h @@ -8,7 +8,7 @@ using namespace Rcpp; using namespace arma; -Rcpp::List sample_Q( +arma::field sample_Q( const int& p, const arma::mat& Y, const arma::mat& X, diff --git a/src/sample_hyper.cpp b/src/sample_hyper.cpp index a48a54e..9e8ee58 100644 --- a/src/sample_hyper.cpp +++ b/src/sample_hyper.cpp @@ -1,137 +1,13 @@ #include +#include "utils.h" +#include "mcmc.h" + 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 -) { - - mat prior_B = as(prior["B"]); - int prior_nu = as(prior["nu"]); - - 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"]); - } - - double lambda; - if (model(2)) { - lambda = hyper(2); - } else { - lambda = 0.2; - } - prior_V *= lambda * lambda; - - return List::create( - _["B"] = prior_B, - _["V"] = prior_V, - _["S"] = prior_S, - _["nu"] = prior_nu - ); -} - - -// extend data with dummy observations -// [[Rcpp:interface(cpp)]] -// [[Rcpp::export]] -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; - int K = X.n_cols; - - mat ybar0 = mean(X.submat(0, 0, p - 1, N - 1), 0); - - 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 Yplus = join_vert(Ystar, Y); - mat Xplus = join_vert(Xstar, X); - - return List::create( - _["Yplus"] = Yplus, - _["Xplus"] = Xplus, - _["Ystar"] = Ystar, - _["Xstar"] = Xstar - ); -} - - // log density of gamma distribution // [[Rcpp:interface(cpp)]] // [[Rcpp::export]] @@ -167,25 +43,31 @@ double log_prior_hyper( const Rcpp::List& prior ) { - double log_prior = 0; - - mat prior_hyper = prior["hyper"]; + double log_prior = 0, shape, scale; if (model(0)) { - log_prior += log_dgamma(hyper(0), prior_hyper(0, 0), prior_hyper(0, 1)); + shape = as(prior["mu.shape"]); + scale = as(prior["mu.scale"]); + log_prior += log_dgamma(hyper(0), shape, scale); } if (model(1)) { - log_prior += log_dgamma(hyper(1), prior_hyper(1, 0), prior_hyper(1, 1)); + shape = as(prior["delta.shape"]); + scale = as(prior["delta.scale"]); + log_prior += log_dgamma(hyper(1), shape, scale); } if (model(2)) { - log_prior += log_dgamma(hyper(2), prior_hyper(2, 0), prior_hyper(2, 1)); + shape = as(prior["lambda.shape"]); + scale = as(prior["lambda.scale"]); + log_prior += log_dgamma(hyper(2), shape, scale); } if (model(3)) { + shape = as(prior["psi.shape"]); + scale = as(prior["psi.scale"]); for (int i = 3; i < hyper.n_elem; i++) { - log_prior += log_dinvgamma(hyper(i), prior_hyper(i, 0), prior_hyper(i, 1)); + log_prior += log_dinvgamma(hyper(i), shape, scale); } } @@ -210,17 +92,14 @@ double log_mvgamma( } -// log marginal likelihood -// notation as in Giannone, Lenza & Primiceri (2014) +// log marginal likelihood, notation as in Giannone, Lenza & Primiceri (2014) // [[Rcpp:interface(cpp)]] // [[Rcpp::export]] 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 ) { @@ -228,18 +107,25 @@ double log_ml( int T = Y.n_rows; int N = Y.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; + double log_ml = 0; + mat inv_Omega = diagmat(1 / Omega.diag()); - 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 + d) / 2 * log_det_sympd(Psi + A); + try { + mat Bhat = inv_sympd(X.t() * X + inv_Omega) * (X.t() * Y + inv_Omega * b); + mat ehat = Y - X * Bhat; + + log_ml += - N * T / 2.0 * log(M_PI); + log_ml += log_mvgamma(N, (T + d) / 2.0); + log_ml += -log_mvgamma(N, d / 2.0); + log_ml += - N / 2.0 * log_det_sympd(Omega); + log_ml += d / 2.0 * log_det_sympd(Psi); + log_ml += - N / 2.0 * log_det_sympd(X.t() * X + inv_Omega); + mat A = Psi + ehat.t() * ehat + (Bhat - b).t() * inv_Omega * (Bhat - b); + log_ml += - (T + d) / 2.0 * log_det_sympd(A); + + } catch(...) { + log_ml = -1e+10; + } return log_ml; } @@ -249,35 +135,39 @@ double log_ml( // [[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 + const Rcpp::List& prior ) { - int N = Y.n_cols; - int n_hyper = hyper.n_elem; - - 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, model, prior); - mat prior_B = as(prior["B"]); - mat prior_V = as(prior["V"]); - mat prior_S = as(prior["S"]); - int prior_nu = as(prior["nu"]); - - mat inv_V = diagmat(1 / prior_V.diag()); - - 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; + int N = Y.n_cols; + double mu = hyper(0); + double delta = hyper(1); + double lambda = hyper(2); + vec psi = hyper.rows(3, N + 2); + + // update Minnesota prior + mat prior_B = as(prior["B"]); + mat prior_V = diagmat(join_vert( + lambda*lambda * kron(as(prior["Vp"]), 1 / psi), + as(prior["Vd"]))); + mat prior_S = diagmat(psi); + int prior_nu = as(prior["nu"]); + + // update dummy observation prior + mat Ystar = join_vert(as(prior["Ysoc"]) / mu, + as(prior["Ysur"]) / delta); + mat Xstar = join_vert(as(prior["Xsoc"]) / mu, + as(prior["Xsur"]) / delta); + mat Yplus = join_vert(Ystar, Y); + mat Xplus = join_vert(Xstar, X); + + double log_ml_plus = log_ml(prior_B, prior_V, prior_S, prior_nu, Yplus, Xplus); + double log_ml_star = log_ml(prior_B, prior_V, prior_S, prior_nu, Ystar, Xstar); + + return log_ml_plus - log_ml_star; } @@ -285,7 +175,6 @@ double log_ml_dummy( // [[Rcpp:interface(cpp)]] // [[Rcpp::export]] double log_posterior_hyper( - const int& p, const arma::vec& hyper, const arma::vec& model, const arma::mat& Y, @@ -294,9 +183,110 @@ double log_posterior_hyper( ) { double log_prior = log_prior_hyper(hyper, model, prior); - double log_ml = log_ml_dummy(p, hyper, model, Y, X, prior); + double log_ml = log_ml_dummy(hyper, model, Y, X, prior); + double log_post = log_prior + log_ml; + + if (!std::isfinite(log_post)) { + log_post = -1e+10; + } + + return log_post; +} + + +// [[Rcpp:interface(cpp)]] +// [[Rcpp::export]] +arma::mat extend_hyper( + const arma::vec& init, + const arma::vec& model, + const arma::mat& hypers +) { + + int i = 0; + + mat extended = repmat(init, 1, hypers.n_cols); + + if (model(0)) { + extended.row(0) = hypers.row(i); + i++; + } + + if (model(1)) { + extended.row(1) = hypers.row(i); + i++; + } + + if (model(2)) { + extended.row(2) = hypers.row(i); + i++; + } + + if (model(3)) { + extended.rows(3, extended.n_rows - 1) = hypers.rows(i, hypers.n_rows - 1); + } + + return extended; +} + + +// [[Rcpp:interface(cpp)]] +// [[Rcpp::export]] +arma::mat narrow_hyper( + const arma::vec& model, + arma::mat hypers +) { + + uvec indices; + + if (!model(0)) { + indices = join_vert(indices, uvec({0})); + } - return log_prior + log_ml; + if (!model(1)) { + indices = join_vert(indices, uvec({1})); + } + + if (!model(2)) { + indices = join_vert(indices, uvec({2})); + } + + if (!model(3)) { + indices = join_vert(indices, regspace(3, hypers.n_rows - 1)); + } + + hypers.shed_rows(indices); + + return hypers; } +// sample hyper-parameters +// [[Rcpp:interface(cpp)]] +// [[Rcpp::export]] +arma::mat sample_hyper( + const int& S, + const int& start, + const arma::vec& init, + const arma::vec& model, + const arma::mat& Y, + const arma::mat& X, + const arma::mat& W, + const Rcpp::List& prior +) { + + mat hypers = metropolis( + + S, start, narrow_hyper(model, init), W, + + [init, model, Y, X, prior](const vec& x) { + + vec extended = extend_hyper(init, model, x); + return log_posterior_hyper(extended, model, Y, X, prior); + } + ); + + hypers = extend_hyper(init, model, hypers); + + return hypers; +} + diff --git a/src/sample_hyper.h b/src/sample_hyper.h index 2657272..ed0320e 100644 --- a/src/sample_hyper.h +++ b/src/sample_hyper.h @@ -7,36 +7,15 @@ 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 -); - - -Rcpp::List extend_dummy( - 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 int& start, + const arma::vec& init, + const arma::vec& model, const arma::mat& Y, const arma::mat& X, + const arma::mat& W, const Rcpp::List& prior ); - #endif // _SAMPLE_HYPER_H_ diff --git a/src/utils.cpp b/src/utils.cpp index 837d294..aaaad02 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -1,7 +1,6 @@ #include #include - #include using namespace arma; @@ -50,9 +49,9 @@ bool match_sign( // Numerical derivative of f: R^n -> R^m, at x (Rcpp::export is not possible) // [[Rcpp:interface(cpp)]] arma::mat Df( - const std::function& f, - const arma::colvec& x, - double h = 1e-10 + const std::function& f, + const arma::vec& x, + const double h = 1e-10 ) { colvec f_x = f(x); @@ -78,12 +77,3 @@ arma::mat Df( return result; } - - - - - - - - - diff --git a/src/utils.h b/src/utils.h index 8749870..1ef8311 100644 --- a/src/utils.h +++ b/src/utils.h @@ -19,4 +19,12 @@ arma::mat Df( double h = 1e-10 ); +arma::mat metropolis( + const int& T, + const int& t0, + arma::vec x, + arma::mat Sigma, + const std::function& log_target +); + #endif // _UTILS_H_