From a5e2ed12517c96b49ad8a627ae0361df0f4287a3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Mon, 24 Jun 2024 15:28:47 +1000 Subject: [PATCH 01/22] get rid of the lm --- NAMESPACE | 1 - R/bsvarSIGNs-package.R | 1 - 2 files changed, 2 deletions(-) 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/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} From c79468376ca457c750c0bba7215fdd60d65c7639 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Mon, 24 Jun 2024 15:29:04 +1000 Subject: [PATCH 02/22] specify_priors fully developed --- R/specify_bsvarSIGN.R | 192 ++++++++++++++++++++++++++------- man/specify_prior_bsvarSIGN.Rd | 55 ++++++++-- 2 files changed, 198 insertions(+), 49 deletions(-) diff --git a/R/specify_bsvarSIGN.R b/R/specify_bsvarSIGN.R index ddc2362..b648f88 100644 --- a/R/specify_bsvarSIGN.R +++ b/R/specify_bsvarSIGN.R @@ -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,12 +139,57 @@ 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 $\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 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 Yplus an \code{(N+1)xN} matrix with the sum-of-coefficients and + #' dummy-initial-observation prior LHS matrix. + Yplus = matrix(), + + #' @field Xplus an \code{(N+1)xK} matrix with the sum-of-coefficients and + #' dummy-initial-observation prior LHS matrix. + Xplus = matrix(), + + #' @field mu.scale a positive scalar - the shape of the gamma prior for $\mu$. + mu.scale = NA, + + #' @field mu.shape a positive scalar - the shape of the gamma prior for $\mu$. + mu.shape = NA, + + #' @field delta.scale a positive scalar - the shape of the gamma prior for $\delta$. + delta.scale = NA, + + #' @field delta.shape a positive scalar - the shape of the gamma prior for $\delta$. + delta.shape = NA, + + #' @field lambda.scale a positive scalar - the shape of the gamma prior for $\lambda$. + lambda.scale = NA, + + #' @field lambda.shape a positive scalar - the shape of the gamma prior for $\lambda$. + lambda.shape = NA, + + #' @field psi.scale a positive scalar - the shape of the inverted gamma prior for $\psi$. + psi.scale = NA, + + #' @field psi.shape a positive scalar - the shape of the inverted gamma prior for $\psi$. + psi.shape = NA, #' @description #' Create a new prior specification PriorBSVAR. @@ -144,27 +204,66 @@ 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, stationary = rep(TRUE, 3)) #' 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, d = 0, stationary = rep(FALSE, ncol(data))){ 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 + 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)) + + T = nrow(Y) + K = N * p + 1 + d - 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) + B = matrix(0, K, N) + B[1:N,] = diag(!stationary) + + V = diag(c(kronecker((1:p)^-2, rep(1, N)), rep(100, 1 + d))) + + 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(Y[1:p,]) + Yplus = rbind(diag(ybar), ybar) + Xplus = Yplus + if (p > 1) { + for (i in 2:p) { + Xplus = cbind(Xplus, Yplus) + } + } + Xplus = cbind(Xplus, c(rep(0, N), 1), matrix(0, N + 1, d)) + + self$hyper = hyper + self$B = B + self$V = V + self$S = diag(N) + self$nu = N + 2 + self$Yplus = Yplus + self$Xplus = Xplus + 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,8 +276,21 @@ specify_prior_bsvarSIGN = R6::R6Class( #' get_prior = function(){ list( - model = self$model, - hyper = self$hyper + hyper = self$hyper, + B = self$B, + V = self$V, + S = self$S, + nu = self$nu, + Yplus = self$Yplus, + Xplus = self$Xplus, + 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 @@ -509,7 +621,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, d, stationary) self$starting_values = bsvars::specify_starting_values_bsvar$new(N, self$p, d) }, # END initialize diff --git a/man/specify_prior_bsvarSIGN.Rd b/man/specify_prior_bsvarSIGN.Rd index 9a82b6d..02b2b03 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, stationary = rep(TRUE, 3)) prior$B # show autoregressive prior mean @@ -35,10 +35,42 @@ prior$get_prior() # show the prior as list \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 $\mu$, $\delta$, +$\lambda$, $\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{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{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{Yplus}}{an \code{(N+1)xN} matrix with the sum-of-coefficients and +dummy-initial-observation prior LHS matrix.} + +\item{\code{Xplus}}{an \code{(N+1)xK} matrix with the sum-of-coefficients and +dummy-initial-observation prior LHS matrix.} + +\item{\code{mu.scale}}{a positive scalar - the shape of the gamma prior for $\mu$.} + +\item{\code{mu.shape}}{a positive scalar - the shape of the gamma prior for $\mu$.} + +\item{\code{delta.scale}}{a positive scalar - the shape of the gamma prior for $\delta$.} + +\item{\code{delta.shape}}{a positive scalar - the shape of the gamma prior for $\delta$.} + +\item{\code{lambda.scale}}{a positive scalar - the shape of the gamma prior for $\lambda$.} + +\item{\code{lambda.shape}}{a positive scalar - the shape of the gamma prior for $\lambda$.} + +\item{\code{psi.scale}}{a positive scalar - the shape of the inverted gamma prior for $\psi$.} + +\item{\code{psi.shape}}{a positive scalar - the shape of the inverted gamma prior for $\psi$.} } \if{html}{\out{
}} } @@ -56,7 +88,12 @@ 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, + d = 0, + stationary = rep(FALSE, ncol(data)) +)}\if{html}{\out{
}} } \subsection{Arguments}{ @@ -64,8 +101,6 @@ Create a new prior specification PriorBSVAR. \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{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.} @@ -73,6 +108,8 @@ Create a new prior specification PriorBSVAR. \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, otherwise to random walk.} + +\item{\code{N}}{a positive integer - the number of dependent variables in the model.} } \if{html}{\out{}} } @@ -83,7 +120,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, stationary = rep(TRUE, 3)) prior$B # show autoregressive prior mean } From 17dce6f303df174e3bd676fc8f91bd9e33c53a33 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Mon, 24 Jun 2024 15:36:02 +1000 Subject: [PATCH 03/22] now it'll work with p=1 --- R/specify_bsvarSIGN.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/specify_bsvarSIGN.R b/R/specify_bsvarSIGN.R index b648f88..129b00f 100644 --- a/R/specify_bsvarSIGN.R +++ b/R/specify_bsvarSIGN.R @@ -239,7 +239,7 @@ specify_prior_bsvarSIGN = R6::R6Class( scale = gamma_scale(1, 1) shape = gamma_shape(1, 1) - ybar = colMeans(Y[1:p,]) + ybar = colMeans(matrix(Y[1:p,], ncol = N)) Yplus = rbind(diag(ybar), ybar) Xplus = Yplus if (p > 1) { From b4286327bb06c2012957e7d85aac45ecf22f689d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Mon, 24 Jun 2024 15:36:20 +1000 Subject: [PATCH 04/22] specify_prior passes all the tests: old and new :) --- inst/tinytest/test_specify.R | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) 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") From bd8d05a111e806c89eac7bed0f1664a1579374b5 Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Mon, 24 Jun 2024 19:47:29 +1000 Subject: [PATCH 05/22] remove unused variable --- src/bsvars_sign.cpp | 20 +++++++------------- 1 file changed, 7 insertions(+), 13 deletions(-) diff --git a/src/bsvars_sign.cpp b/src/bsvars_sign.cpp index 825b856..c64b874 100644 --- a/src/bsvars_sign.cpp +++ b/src/bsvars_sign.cpp @@ -52,19 +52,11 @@ Rcpp::List bsvar_sign_cpp( } Progress p(50, show_progress); - 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"]); + const int T = Y.n_rows; + const int N = Y.n_cols; + const int K = X.n_cols; 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); @@ -72,7 +64,9 @@ Rcpp::List bsvar_sign_cpp( cube posterior_Theta0(N, N, S); cube posterior_shocks(N, T, S); - int s = 0, S_hyper = sample_hyper.n_cols - 1, post_nu; + mat hypers = as(prior["hyper"]); + + int s = 0, S_hyper = hypers.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; @@ -83,7 +77,7 @@ Rcpp::List bsvar_sign_cpp( // Check for user interrupts if (s % 200 == 0) checkUserInterrupt(); - hyper = sample_hyper.col(randi(distr_param(0, S_hyper))); + hyper = hypers.col(randi(distr_param(0, S_hyper))); lambda = hyper(2); psi = hyper.rows(3, N + 2); From f0e1743bd36e84b7343a2948ba181e6349a99c6f Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Sat, 29 Jun 2024 21:59:36 +1000 Subject: [PATCH 06/22] GLP draft code --- inst/varia/nicetry.R | 166 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 166 insertions(+) create mode 100644 inst/varia/nicetry.R diff --git a/inst/varia/nicetry.R b/inst/varia/nicetry.R new file mode 100644 index 0000000..05319ed --- /dev/null +++ b/inst/varia/nicetry.R @@ -0,0 +1,166 @@ +rm(list = ls()) + +niw = function(Y, + p, + non_stationary, + lambda, + psi) { + T = nrow(Y) + N = ncol(Y) + K = 1 + N * p + + B = matrix(0, K, N) + B[1:N, 1:N] = diag(non_stationary) + + V = matrix(0, K, K) + V[K, K] = 1e+6 + V[1:(K - 1), 1:(K - 1)] = diag(lambda^2 * kronecker((1:p)^-2, psi^-1)) + + S = diag(psi) + nu = N + 2 + + list(B = B, V = V, S = S, nu = nu) +} + +log_det = function(A) { + determinant(A, logarithm = TRUE)$modulus +} + +dhyper = function(p, hyper, model, Y, X) { + N = ncol(Y) + + if (model[3]) { + lambda = hyper[3] + } else { + lambda = 0.2 + } + + if (model[4]) { + psi = hyper[4:length(hyper)] + } else { + psi = sapply(1:N, \(i) summary(lm(Y[2:T, i] ~ Y[1:(T - 1), i]))$sigma^2) + } + + prior = niw(Y, p, rep(1, N), lambda, psi) + + extended = extend_dummy(p, hyper, model, Y, X) + Ystar = extended$Ystar + Xstar = extended$Xstar + Yplus = extended$Yplus + Xplus = extended$Xplus + + logml_plus = logml(prior, p, Yplus, Xplus) + if (dim(Ystar)[1] > 0) { + logml_star = logml(prior, p, Ystar, Xstar) + } else { + logml_star = 0 + } + + as.numeric(logml_plus - logml_star) +} + +logml = function(prior, p, Y, X) { + + N = ncol(Y) + T = nrow(Y) + + b = prior$B + Omega = prior$V + Psi = prior$S + d = prior$nu + + llike = -1e10 + tryCatch({ + inv_Omega = solve(Omega) + Bhat = solve(t(X) %*% X + inv_Omega, t(X) %*% Y + inv_Omega %*% b) + ehat = Y - X %*% Bhat + + llike = N * T / 2 * log(1 / pi) + log_mvgamma(N, (T + d) / 2) - log_mvgamma(N, d / 2) + llike = llike - N / 2 * log_det(Omega) + llike = llike + d / 2 * log_det(Psi) + 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) + }, error = function(e) { + }) + + if (!is.finite(llike)) { + llike = -1e10 + } + + return(llike) +} + +set.seed(123) + +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) +T = nrow(Y) +N = ncol(Y) + +init = c(1, 1, 0.2) +model = c(TRUE, TRUE, TRUE, FALSE) + +init = c(init, rep(0.02^2, N)) +model[4] = TRUE + +result = optim(init, + \(x) -dhyper(p, x, model, Y, X), + method = 'L-BFGS-B', + control = list(trace = 1), + lower = rep(0, length(init)), + upper = init * 100, + hessian = TRUE + ) + +h = result$par +d = dhyper(p, h, model, Y, X) +h = log(h) + +c = 1 +W = result$hessian +e = eigen(W) +W = e$vectors %*% diag(1 / abs(e$values)) %*% t(e$vectors) + +S = 10000 +start = 1000 +hyper = matrix(0, S, length(init)) +hyper[1, ] = h + +for (s in 2:S) { + newh = mvtnorm::rmvnorm(1, h, c^2 * W) + newd = dhyper(p, exp(newh), model, Y, X) + a = min(1, exp(newd - d + sum(h) - sum(newh))) + + if (runif(1) < a) { + h = newh + d = newd + } + + if (s > start) { + W = var(hyper[1:s, ]) + c = c + (s - start) ^ -0.6 * (a - 0.234) + } + + if (s %% 100 == 0) { + cat(s, "\r") + } + + hyper[s, ] = h +} + +burn = floor(S / 2) +plot.ts(exp(hyper[burn:S, ])) +# hist(exp(hyper[burn:S, 3]), breaks = 100) + + + From e1018e6cace478a0118309bf6dec44095d732b78 Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Mon, 1 Jul 2024 01:01:03 +1000 Subject: [PATCH 07/22] log ml --- R/RcppExports.R | 16 +++--- inst/varia/nicetry.R | 26 ++++++---- src/RcppExports.cpp | 42 +++++----------- src/sample_hyper.cpp | 117 ++++++++++++++----------------------------- src/sample_hyper.h | 17 ------- 5 files changed, 71 insertions(+), 147 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 12cb357..dab637d 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -17,10 +17,6 @@ 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) } @@ -41,16 +37,16 @@ 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(p, b, Omega, Psi, d, Y, X) { + .Call(`_bsvarSIGNs_log_ml`, p, b, Omega, Psi, d, Y, X) } -log_ml_dummy <- function(p, hyper, model, Y, X, prior) { - .Call(`_bsvarSIGNs_log_ml_dummy`, p, hyper, model, Y, X, prior) +log_ml_dummy <- function(p, hyper, model, Y, X) { + .Call(`_bsvarSIGNs_log_ml_dummy`, p, hyper, model, Y, X) } -log_posterior_hyper <- function(p, hyper, model, Y, X, prior) { - .Call(`_bsvarSIGNs_log_posterior_hyper`, p, hyper, model, Y, X, prior) +log_posterior_hyper <- function(p, hyper, model, Y, X) { + .Call(`_bsvarSIGNs_log_posterior_hyper`, p, hyper, model, Y, X) } # Register entry points for exported C++ functions diff --git a/inst/varia/nicetry.R b/inst/varia/nicetry.R index 05319ed..08671e8 100644 --- a/inst/varia/nicetry.R +++ b/inst/varia/nicetry.R @@ -75,7 +75,8 @@ logml = function(prior, p, Y, X) { Bhat = solve(t(X) %*% X + inv_Omega, t(X) %*% Y + inv_Omega %*% b) ehat = Y - X %*% Bhat - llike = N * T / 2 * log(1 / pi) + log_mvgamma(N, (T + d) / 2) - log_mvgamma(N, d / 2) + llike = - N * T / 2 * log(pi) + llike = llike + log_mvgamma(N, (T + d) / 2) - log_mvgamma(N, d / 2) llike = llike - N / 2 * log_det(Omega) llike = llike + d / 2 * log_det(Psi) llike = llike - N / 2 * log_det(t(X) %*% X + inv_Omega) @@ -85,7 +86,7 @@ logml = function(prior, p, Y, X) { }) if (!is.finite(llike)) { - llike = -1e10 + llike = -1e+10 } return(llike) @@ -93,6 +94,8 @@ logml = function(prior, p, Y, X) { set.seed(123) +dhyper = log_posterior_hyper + library(R.matlab) data = readMat('../data/DataSW.mat') name = sapply(1:7, \(i) data$ShortDescr[[i]][[1]]) @@ -116,7 +119,7 @@ model[4] = TRUE result = optim(init, \(x) -dhyper(p, x, model, Y, X), method = 'L-BFGS-B', - control = list(trace = 1), + control = list(trace = 1, maxit = 1e5), lower = rep(0, length(init)), upper = init * 100, hessian = TRUE @@ -130,12 +133,15 @@ c = 1 W = result$hessian e = eigen(W) W = e$vectors %*% diag(1 / abs(e$values)) %*% t(e$vectors) +W = MASS::ginv(result$hessian) +W = (W + t(W)) / 2 S = 10000 -start = 1000 +start = S / 4 hyper = matrix(0, S, length(init)) hyper[1, ] = h +start_time = Sys.time() for (s in 2:S) { newh = mvtnorm::rmvnorm(1, h, c^2 * W) newd = dhyper(p, exp(newh), model, Y, X) @@ -148,7 +154,7 @@ for (s in 2:S) { if (s > start) { W = var(hyper[1:s, ]) - c = c + (s - start) ^ -0.6 * (a - 0.234) + c = c + s ^ -0.6 * (a - 0.234) } if (s %% 100 == 0) { @@ -157,10 +163,8 @@ for (s in 2:S) { hyper[s, ] = h } +end_time = Sys.time() +end_time - start_time -burn = floor(S / 2) -plot.ts(exp(hyper[burn:S, ])) -# hist(exp(hyper[burn:S, 3]), breaks = 100) - - - +hyper = exp(hyper) +plot.ts(hyper) diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 1dd9525..348cb15 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -627,20 +627,6 @@ BEGIN_RCPP 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) { @@ -708,8 +694,8 @@ 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 int& p, 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 pSEXP, SEXP bSEXP, SEXP OmegaSEXP, SEXP PsiSEXP, SEXP dSEXP, SEXP YSEXP, SEXP XSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -718,16 +704,15 @@ BEGIN_RCPP 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(p, 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 int& p, const arma::vec& hyper, const arma::vec& model, const arma::mat& Y, const arma::mat& X); +RcppExport SEXP _bsvarSIGNs_log_ml_dummy(SEXP pSEXP, SEXP hyperSEXP, SEXP modelSEXP, SEXP YSEXP, SEXP XSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -736,14 +721,13 @@ BEGIN_RCPP 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_result_gen = Rcpp::wrap(log_ml_dummy(p, hyper, model, Y, X)); 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 int& p, const arma::vec& hyper, const arma::vec& model, const arma::mat& Y, const arma::mat& X); +RcppExport SEXP _bsvarSIGNs_log_posterior_hyper(SEXP pSEXP, SEXP hyperSEXP, SEXP modelSEXP, SEXP YSEXP, SEXP XSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -752,8 +736,7 @@ BEGIN_RCPP 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(p, hyper, model, Y, X)); return rcpp_result_gen; END_RCPP } @@ -931,15 +914,14 @@ static const R_CallMethodDef CallEntries[] = { {"_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, 7}, + {"_bsvarSIGNs_log_ml_dummy", (DL_FUNC) &_bsvarSIGNs_log_ml_dummy, 5}, + {"_bsvarSIGNs_log_posterior_hyper", (DL_FUNC) &_bsvarSIGNs_log_posterior_hyper, 5}, {"_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 a48a54e..d508dd3 100644 --- a/src/sample_hyper.cpp +++ b/src/sample_hyper.cpp @@ -36,53 +36,6 @@ Rcpp::List mn_prior( } -// 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]] @@ -210,8 +163,7 @@ 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( @@ -220,7 +172,6 @@ double log_ml( 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 +179,28 @@ 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_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); + + if (!std::isfinite(log_ml)) { + log_ml = -1e+10; + } + + } catch(...) { + log_ml = -1e+10; + } return log_ml; } @@ -253,31 +214,29 @@ double log_ml_dummy( const arma::vec& hyper, const arma::vec& model, const arma::mat& Y, - const arma::mat& X, - Rcpp::List prior + const arma::mat& X ) { int N = Y.n_cols; - int n_hyper = hyper.n_elem; - List extended = extend_dummy(p, model, hyper, Y, X); + List extended = extend_dummy(p, hyper, model, 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 lambda = hyper(2); + vec psi = hyper.rows(3, N + 2); + List prior = mn_prior(p, lambda, psi); + mat prior_B = as(prior["B"]); + mat prior_V = as(prior["V"]); + mat prior_S = as(prior["S"]); + int prior_nu = as(prior["nu"]); - 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); + double log_ml_plus = log_ml(p, prior_B, prior_V, prior_S, prior_nu, Yplus, Xplus); + double log_ml_star = log_ml(p, prior_B, prior_V, prior_S, prior_nu, Ystar, Xstar); - return log_ml_extended - log_ml_dummy; + return log_ml_plus - log_ml_star; } @@ -289,12 +248,12 @@ double log_posterior_hyper( const arma::vec& hyper, const arma::vec& model, const arma::mat& Y, - const arma::mat& X, - const Rcpp::List& prior + const arma::mat& X ) { - double log_prior = log_prior_hyper(hyper, model, prior); - double log_ml = log_ml_dummy(p, hyper, model, Y, X, prior); + // double log_prior = log_prior_hyper(hyper, model, prior); + double log_prior = 0; + double log_ml = log_ml_dummy(p, hyper, model, Y, X); return log_prior + log_ml; } diff --git a/src/sample_hyper.h b/src/sample_hyper.h index 2657272..1417503 100644 --- a/src/sample_hyper.h +++ b/src/sample_hyper.h @@ -14,13 +14,6 @@ Rcpp::List mn_prior( ); -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, @@ -29,14 +22,4 @@ 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 -); - - #endif // _SAMPLE_HYPER_H_ From 7223e363f673ac78dc129f40855930aa69d5bd99 Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Mon, 1 Jul 2024 15:37:08 +1000 Subject: [PATCH 08/22] use last acceptance rate --- inst/varia/nicetry.R | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/inst/varia/nicetry.R b/inst/varia/nicetry.R index 08671e8..f81bb50 100644 --- a/inst/varia/nicetry.R +++ b/inst/varia/nicetry.R @@ -143,11 +143,16 @@ hyper[1, ] = h start_time = Sys.time() for (s in 2:S) { + + if (s %% 100 == 0) { + cat(s, "\r") + } + newh = mvtnorm::rmvnorm(1, h, c^2 * W) newd = dhyper(p, exp(newh), model, Y, X) - a = min(1, exp(newd - d + sum(h) - sum(newh))) + newa = min(1, exp(newd - d + sum(h) - sum(newh))) - if (runif(1) < a) { + if (runif(1) < newa) { h = newh d = newd } @@ -157,10 +162,7 @@ for (s in 2:S) { c = c + s ^ -0.6 * (a - 0.234) } - if (s %% 100 == 0) { - cat(s, "\r") - } - + a = newa hyper[s, ] = h } end_time = Sys.time() From b90cc44cea2fb35b8e30bc85d43ae42e476d230b Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Fri, 12 Jul 2024 14:16:52 +1000 Subject: [PATCH 09/22] greek letters --- R/specify_bsvarSIGN.R | 21 ++++++++++----------- man/specify_prior_bsvarSIGN.Rd | 19 +++++++++---------- 2 files changed, 19 insertions(+), 21 deletions(-) diff --git a/R/specify_bsvarSIGN.R b/R/specify_bsvarSIGN.R index 129b00f..150c2a6 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) @@ -139,8 +139,7 @@ specify_prior_bsvarSIGN = R6::R6Class( "PriorBSVARSIGN", public = list( - #' @field hyper a \code{(N+3)xS} matrix of hyper-parameters $\mu$, $\delta$, - #' $\lambda$, $\psi$. + #' @field hyper a \code{(N+3)xS} matrix of hyper-parameters \eqn{\mu, \delta, \lambda, \psi}. hyper = matrix(), #' @field B a \code{KxN} normal prior mean matrix for the autoregressive @@ -167,28 +166,28 @@ specify_prior_bsvarSIGN = R6::R6Class( #' dummy-initial-observation prior LHS matrix. Xplus = matrix(), - #' @field mu.scale a positive scalar - the shape of the gamma prior for $\mu$. + #' @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 $\mu$. + #' @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 $\delta$. + #' @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 $\delta$. + #' @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 $\lambda$. + #' @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 $\lambda$. + #' @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 $\psi$. + #' @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 $\psi$. + #' @field psi.shape a positive scalar - the shape of the inverted gamma prior for \eqn{\psi}. psi.shape = NA, #' @description diff --git a/man/specify_prior_bsvarSIGN.Rd b/man/specify_prior_bsvarSIGN.Rd index 02b2b03..d425c66 100644 --- a/man/specify_prior_bsvarSIGN.Rd +++ b/man/specify_prior_bsvarSIGN.Rd @@ -35,8 +35,7 @@ prior$get_prior() # show the prior as list \section{Public fields}{ \if{html}{\out{
}} \describe{ -\item{\code{hyper}}{a \code{(N+3)xS} matrix of hyper-parameters $\mu$, $\delta$, -$\lambda$, $\psi$.} +\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.} @@ -56,21 +55,21 @@ dummy-initial-observation prior LHS matrix.} \item{\code{Xplus}}{an \code{(N+1)xK} matrix with the sum-of-coefficients and dummy-initial-observation prior LHS matrix.} -\item{\code{mu.scale}}{a positive scalar - the shape of the gamma prior for $\mu$.} +\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 $\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 $\delta$.} +\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 $\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 $\lambda$.} +\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 $\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 $\psi$.} +\item{\code{psi.scale}}{a positive scalar - the shape of the inverted gamma prior for \eqn{\psi}.} -\item{\code{psi.shape}}{a positive scalar - the shape of the inverted gamma prior for $\psi$.} +\item{\code{psi.shape}}{a positive scalar - the shape of the inverted gamma prior for \eqn{\psi}.} } \if{html}{\out{
}} } From f8dfc5ff18bee1e65f386843f4dd38bdb5a82ed5 Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Fri, 12 Jul 2024 14:18:39 +1000 Subject: [PATCH 10/22] adaptive metropolis --- R/RcppExports.R | 12 +++++ inst/varia/nicetry.R | 67 +++++++++++----------------- src/RcppExports.cpp | 46 +++++++++++++++++++ src/sample_hyper.cpp | 103 +++++++++++++++++++++++++++++++++++++++++++ src/sample_hyper.h | 11 +++++ src/utils.cpp | 49 ++++++++++++++++++-- src/utils.h | 8 ++++ 7 files changed, 251 insertions(+), 45 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index dab637d..2079232 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -49,6 +49,18 @@ log_posterior_hyper <- function(p, hyper, model, Y, X) { .Call(`_bsvarSIGNs_log_posterior_hyper`, p, hyper, model, Y, X) } +extend_hyper <- function(init, model, hypers) { + .Call(`_bsvarSIGNs_extend_hyper`, init, model, hypers) +} + +narrow_hyper <- function(model, hypers) { + .Call(`_bsvarSIGNs_narrow_hyper`, model, hypers) +} + +sample_hyper <- function(S, start, p, init, model, Y, X, W) { + .Call(`_bsvarSIGNs_sample_hyper`, S, start, p, init, model, Y, X, W) +} + # Register entry points for exported C++ functions methods::setLoadAction(function(ns) { .Call(`_bsvarSIGNs_RcppExport_registerCCallable`) diff --git a/inst/varia/nicetry.R b/inst/varia/nicetry.R index f81bb50..1c76b6d 100644 --- a/inst/varia/nicetry.R +++ b/inst/varia/nicetry.R @@ -111,13 +111,12 @@ T = nrow(Y) N = ncol(Y) init = c(1, 1, 0.2) -model = c(TRUE, TRUE, TRUE, FALSE) +init = c(init, rep(0.02^2, N)) +model = c(FALSE, FALSE, TRUE, TRUE) +# model[4] = TRUE -init = c(init, rep(0.02^2, N)) -model[4] = TRUE - -result = optim(init, - \(x) -dhyper(p, x, model, Y, X), +result = optim(narrow_hyper(model, matrix(init)), + \(x) -dhyper(p, extend_hyper(init, model, matrix(x)), model, Y, X), method = 'L-BFGS-B', control = list(trace = 1, maxit = 1e5), lower = rep(0, length(init)), @@ -126,47 +125,33 @@ result = optim(init, ) h = result$par -d = dhyper(p, h, model, Y, X) +d = dhyper(p, extend_hyper(init, model, matrix(h)), model, Y, X) h = log(h) + c = 1 W = result$hessian -e = eigen(W) -W = e$vectors %*% diag(1 / abs(e$values)) %*% t(e$vectors) -W = MASS::ginv(result$hessian) -W = (W + t(W)) / 2 +if (length(h) == 1){ + W = 1 / W +} else { + e = eigen(W) + W = e$vectors %*% diag(as.vector(1 / abs(e$values))) %*% t(e$vectors) +} -S = 10000 -start = S / 4 -hyper = matrix(0, S, length(init)) -hyper[1, ] = h +S = 10000 +start = S / 10 start_time = Sys.time() -for (s in 2:S) { - - if (s %% 100 == 0) { - cat(s, "\r") - } - - newh = mvtnorm::rmvnorm(1, h, c^2 * W) - newd = dhyper(p, exp(newh), model, Y, X) - newa = min(1, exp(newd - d + sum(h) - sum(newh))) - - if (runif(1) < newa) { - h = newh - d = newd - } - - if (s > start) { - W = var(hyper[1:s, ]) - c = c + s ^ -0.6 * (a - 0.234) - } - - a = newa - hyper[s, ] = h -} -end_time = Sys.time() -end_time - start_time -hyper = exp(hyper) + +hyper = sample_hyper(S, start, p, + extend_hyper(init, model, matrix(result$par)), + model, Y, X, W) + + +hyper = t(hyper) plot.ts(hyper) + + +Sys.time() - start_time + diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 348cb15..6ed481a 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -740,6 +740,49 @@ BEGIN_RCPP 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 int& p, const arma::vec& init, const arma::vec& model, const arma::mat& Y, const arma::mat& X, const arma::mat& W); +RcppExport SEXP _bsvarSIGNs_sample_hyper(SEXP SSEXP, SEXP startSEXP, SEXP pSEXP, SEXP initSEXP, SEXP modelSEXP, SEXP YSEXP, SEXP XSEXP, SEXP WSEXP) { +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 int& >::type p(pSEXP); + 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_result_gen = Rcpp::wrap(sample_hyper(S, start, p, init, model, Y, X, W)); + return rcpp_result_gen; +END_RCPP +} // qr_sign_cpp arma::mat qr_sign_cpp(const arma::mat& A); static SEXP _bsvarSIGNs_qr_sign_cpp_try(SEXP ASEXP) { @@ -922,6 +965,9 @@ static const R_CallMethodDef CallEntries[] = { {"_bsvarSIGNs_log_ml", (DL_FUNC) &_bsvarSIGNs_log_ml, 7}, {"_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/sample_hyper.cpp b/src/sample_hyper.cpp index d508dd3..dbaa4f7 100644 --- a/src/sample_hyper.cpp +++ b/src/sample_hyper.cpp @@ -1,6 +1,8 @@ #include +#include "utils.h" + using namespace Rcpp; using namespace arma; @@ -259,3 +261,104 @@ double log_posterior_hyper( } +// [[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})); + } + + 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 int& p, + const arma::vec& init, + const arma::vec& model, + const arma::mat& Y, + const arma::mat& X, + const arma::mat& W +) { + + mat hypers = metropolis( + S, start, narrow_hyper(model, init), W, + [p, init, model, Y, X](const vec& x) { + vec extended = extend_hyper(init, model, x); + return log_posterior_hyper(p, extended, model, Y, X); + } + ); + + hypers = extend_hyper(init, model, hypers); + + return hypers; +} + + + + + + + + + diff --git a/src/sample_hyper.h b/src/sample_hyper.h index 1417503..399fbb1 100644 --- a/src/sample_hyper.h +++ b/src/sample_hyper.h @@ -22,4 +22,15 @@ Rcpp::List extend_dummy( ); +arma::mat sample_hyper( + const int& S, + const int& start, + const int& p, + const arma::vec& init, + const arma::vec& model, + const arma::mat& Y, + const arma::mat& X, + const arma::mat& W +); + #endif // _SAMPLE_HYPER_H_ diff --git a/src/utils.cpp b/src/utils.cpp index 837d294..e0e9a3d 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -50,9 +50,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); @@ -79,7 +79,48 @@ arma::mat Df( } - +// 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; + + double s = 1.0; + double d = log_target(x); + + double new_d, a; + vec new_x; + + mat X(n, T); + x = log(x); + X.col(0) = x; + + for (int t = 1; t < T; t++) { + new_x = mvnrnd(x, s*s * Sigma); + new_d = log_target(exp(new_x)); + a = std::min(1.0, exp(new_d - d + sum(x - new_x))); + + if (randu() < a) { + x = new_x; + d = new_d; + } + X.col(t) = x; + + if (t >= t0) { + Sigma = cov(X.cols(0, t).t()); + s += pow(t, -0.6) * (a - 0.234); + } + } + + X = exp(X); + return X; +} 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_ From 5c680332235067e1ace76911374ee02d610c270d Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Fri, 12 Jul 2024 23:38:39 +1000 Subject: [PATCH 11/22] return field instead of list --- src/sample_NIW.cpp | 27 +++++++++++++-------------- src/sample_NIW.h | 15 ++++++++++----- 2 files changed, 23 insertions(+), 19 deletions(-) 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_ From 5955289a7beb9aa63fa5a30be5767e97d095c3f8 Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Fri, 12 Jul 2024 23:39:07 +1000 Subject: [PATCH 12/22] log prior --- R/RcppExports.R | 20 ++-- R/specify_bsvarSIGN.R | 167 +++++++++++++++++++++++++++------ inst/varia/nicetry.R | 165 +++++--------------------------- man/specify_prior_bsvarSIGN.Rd | 90 ++++++++++++++++-- src/RcppExports.cpp | 46 ++++----- src/bsvars_sign.cpp | 40 ++++++-- src/sample_hyper.cpp | 90 +++++++++++------- src/sample_hyper.h | 4 +- 8 files changed, 367 insertions(+), 255 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 2079232..79b1fd5 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -9,8 +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) +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) } mn_prior <- function(p, lambda, psi) { @@ -37,16 +37,16 @@ log_mvgamma <- function(n, x) { .Call(`_bsvarSIGNs_log_mvgamma`, n, x) } -log_ml <- function(p, b, Omega, Psi, d, Y, X) { - .Call(`_bsvarSIGNs_log_ml`, p, b, Omega, Psi, d, 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(p, hyper, model, Y, X) { - .Call(`_bsvarSIGNs_log_ml_dummy`, p, hyper, model, 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(p, hyper, model, Y, X) { - .Call(`_bsvarSIGNs_log_posterior_hyper`, p, hyper, model, Y, X) +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) { @@ -57,8 +57,8 @@ narrow_hyper <- function(model, hypers) { .Call(`_bsvarSIGNs_narrow_hyper`, model, hypers) } -sample_hyper <- function(S, start, p, init, model, Y, X, W) { - .Call(`_bsvarSIGNs_sample_hyper`, S, start, p, init, model, Y, X, W) +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/specify_bsvarSIGN.R b/R/specify_bsvarSIGN.R index 150c2a6..76c6ac5 100644 --- a/R/specify_bsvarSIGN.R +++ b/R/specify_bsvarSIGN.R @@ -131,7 +131,7 @@ igamma_shape = function(mode, variance) { #' @examples #' # a prior for 3-variable example with one lag #' data(oil) -#' prior = specify_prior_bsvarSIGN$new(oil, p = 1) +#' prior = specify_prior_bsvarSIGN$new(oil, p = 1, stationary = rep(TRUE, 3)) #' prior$B # show autoregressive prior mean #' #' @export @@ -150,6 +150,12 @@ specify_prior_bsvarSIGN = R6::R6Class( #' 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(), @@ -158,13 +164,26 @@ specify_prior_bsvarSIGN = R6::R6Class( #' inverted-Wishart prior for error terms covariance matrix. nu = NA, - #' @field Yplus an \code{(N+1)xN} matrix with the sum-of-coefficients and - #' dummy-initial-observation prior LHS matrix. - Yplus = matrix(), + #' @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 Xplus an \code{(N+1)xK} matrix with the sum-of-coefficients and - #' dummy-initial-observation prior LHS matrix. - Xplus = 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, @@ -192,10 +211,9 @@ specify_prior_bsvarSIGN = R6::R6Class( #' @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. @@ -206,20 +224,30 @@ specify_prior_bsvarSIGN = R6::R6Class( #' prior = specify_prior_bsvarSIGN$new(oil, p = 1, stationary = rep(TRUE, 3)) #' prior$B # show autoregressive prior mean #' - initialize = function(data, p, d = 0, stationary = rep(FALSE, ncol(data))){ + initialize = function(data, p, exogenous = NULL, stationary = rep(FALSE, ncol(data))) { + stopifnot("Argument p must be a positive integer number." = p > 0 & p %% 1 == 0) - stopifnot("Argument d must be a non-negative integer number." = d >= 0 & d %% 1 == 0) - Y = data + + data = bsvars::specify_data_matrices$new(data, p, exogenous) + Y = t(data$Y) + X = t(data$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) - V = diag(c(kronecker((1:p)^-2, rep(1, N)), rep(100, 1 + d))) + 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) { @@ -228,7 +256,8 @@ specify_prior_bsvarSIGN = R6::R6Class( 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) + 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) @@ -239,22 +268,35 @@ specify_prior_bsvarSIGN = R6::R6Class( shape = gamma_shape(1, 1) ybar = colMeans(matrix(Y[1:p,], ncol = N)) - Yplus = rbind(diag(ybar), ybar) - Xplus = Yplus - if (p > 1) { - for (i in 2:p) { - Xplus = cbind(Xplus, Yplus) - } - } - Xplus = cbind(Xplus, c(rep(0, N), 1), matrix(0, N + 1, d)) + 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$Yplus = Yplus - self$Xplus = Xplus + self$Ysoc = Ysoc + self$Xsoc = Xsoc + self$Ysur = Ysur + self$Xsur = Xsur self$mu.scale = scale self$mu.shape = shape self$delta.scale = scale @@ -278,10 +320,14 @@ specify_prior_bsvarSIGN = R6::R6Class( hyper = self$hyper, B = self$B, V = self$V, + Vp = self$Vp, + Vd = self$Vd, S = self$S, nu = self$nu, - Yplus = self$Yplus, - Xplus = self$Xplus, + 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, @@ -291,7 +337,70 @@ specify_prior_bsvarSIGN = R6::R6Class( 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 shrinkage in the Minnesota prior. + #' @param psi whether to estimate the hyper-parameter of variances in the Minnesota prior. + #' @param S number of draws. + #' @param start starting point of the adaptive Metropolis algorithm. + #' @param burn 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, stationary = rep(TRUE, 3)) + #' prior$estimate_hyper(S = 5) + #' + estimate_hyper = function(mu = FALSE, delta = FALSE, lambda = TRUE, psi = FALSE, + S = 10000, start = 1000, burn = 5000){ + + 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() + + if (!mu) { + prior$Ysoc = matrix(0, 0, ncol(self$Ysoc)) + prior$Xsoc = matrix(0, 0, ncol(self$Xsoc)) + } + if (!delta) { + prior$Ysur = matrix(0, 0, ncol(self$Ysur)) + prior$Xsur = matrix(0, 0, ncol(self$Xsur)) + } + + 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)) + W = result$hessian + + if (length(init) == 1){ + W = 1 / W + } else { + e = eigen(W) + W = e$vectors %*% diag(as.vector(1 / abs(e$values))) %*% t(e$vectors) + } + + self$hyper = sample_hyper(S, start, mode, model, self$Y, self$X, W, prior) + self$hyper = self$hyper[, -(1:burn)] + } # END estimate_hyper ) # END public ) # END specify_prior_bsvarSIGN @@ -620,7 +729,7 @@ specify_bsvarSIGN = R6::R6Class( sign_B, zero_irf, max_tries) - self$prior = specify_prior_bsvarSIGN$new(data, p, d, 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/inst/varia/nicetry.R b/inst/varia/nicetry.R index 1c76b6d..b23a652 100644 --- a/inst/varia/nicetry.R +++ b/inst/varia/nicetry.R @@ -1,157 +1,44 @@ -rm(list = ls()) - -niw = function(Y, - p, - non_stationary, - lambda, - psi) { - T = nrow(Y) - N = ncol(Y) - K = 1 + N * p - - B = matrix(0, K, N) - B[1:N, 1:N] = diag(non_stationary) - - V = matrix(0, K, K) - V[K, K] = 1e+6 - V[1:(K - 1), 1:(K - 1)] = diag(lambda^2 * kronecker((1:p)^-2, psi^-1)) - - S = diag(psi) - nu = N + 2 - - list(B = B, V = V, S = S, nu = nu) -} - -log_det = function(A) { - determinant(A, logarithm = TRUE)$modulus -} - -dhyper = function(p, hyper, model, Y, X) { - N = ncol(Y) - - if (model[3]) { - lambda = hyper[3] - } else { - lambda = 0.2 - } - - if (model[4]) { - psi = hyper[4:length(hyper)] - } else { - psi = sapply(1:N, \(i) summary(lm(Y[2:T, i] ~ Y[1:(T - 1), i]))$sigma^2) - } - - prior = niw(Y, p, rep(1, N), lambda, psi) - - extended = extend_dummy(p, hyper, model, Y, X) - Ystar = extended$Ystar - Xstar = extended$Xstar - Yplus = extended$Yplus - Xplus = extended$Xplus - - logml_plus = logml(prior, p, Yplus, Xplus) - if (dim(Ystar)[1] > 0) { - logml_star = logml(prior, p, Ystar, Xstar) - } else { - logml_star = 0 - } - - as.numeric(logml_plus - logml_star) -} - -logml = function(prior, p, Y, X) { - - N = ncol(Y) - T = nrow(Y) - - b = prior$B - Omega = prior$V - Psi = prior$S - d = prior$nu - - llike = -1e10 - tryCatch({ - inv_Omega = solve(Omega) - Bhat = solve(t(X) %*% X + inv_Omega, t(X) %*% Y + inv_Omega %*% b) - ehat = Y - X %*% Bhat - - llike = - N * T / 2 * log(pi) - llike = llike + log_mvgamma(N, (T + d) / 2) - log_mvgamma(N, d / 2) - llike = llike - N / 2 * log_det(Omega) - llike = llike + d / 2 * log_det(Psi) - 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) - }, error = function(e) { - }) - - if (!is.finite(llike)) { - llike = -1e+10 - } - - return(llike) -} - -set.seed(123) - -dhyper = log_posterior_hyper - 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) -T = nrow(Y) -N = ncol(Y) - -init = c(1, 1, 0.2) -init = c(init, rep(0.02^2, N)) -model = c(FALSE, FALSE, TRUE, TRUE) -# model[4] = TRUE -result = optim(narrow_hyper(model, matrix(init)), - \(x) -dhyper(p, extend_hyper(init, model, matrix(x)), model, Y, X), - method = 'L-BFGS-B', - control = list(trace = 1, maxit = 1e5), - lower = rep(0, length(init)), - upper = init * 100, - hessian = TRUE - ) +# data = optimism -h = result$par -d = dhyper(p, extend_hyper(init, model, matrix(h)), model, Y, X) -h = log(h) +spec = specify_bsvarSIGN$new(data, p = 4) +start_time = Sys.time() -c = 1 -W = result$hessian -if (length(h) == 1){ - W = 1 / W -} else { - e = eigen(W) - W = e$vectors %*% diag(as.vector(1 / abs(e$values))) %*% t(e$vectors) -} +spec$prior$estimate_hyper(TRUE, TRUE, TRUE, TRUE, S = 20000, burn = 10000, start = 2000) +# 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 -S = 10000 -start = S / 10 -start_time = Sys.time() +plot.ts(t(spec$prior$hyper)) -hyper = sample_hyper(S, start, p, - extend_hyper(init, model, matrix(result$par)), - model, Y, X, W) -hyper = t(hyper) -plot.ts(hyper) +data(optimism) +# optimism shock +# no contemporaneous effect on productivity +zero_irf = matrix(0, nrow = 5, ncol = 5) +zero_irf[1, 1] = 1 +# positive contemporaneous effect on stock prices +sign_irf = array(0, dim = c(5, 5, 1)) +sign_irf[2, 1, 1] = 1 -Sys.time() - start_time +specification = specify_bsvarSIGN$new(optimism, + p = 4, + sign_irf = sign_irf, + zero_irf = zero_irf) +specification$prior$estimate_hyper(FALSE, TRUE, TRUE, TRUE) +posterior = estimate(specification, S = 10000) +irf = compute_impulse_responses(posterior, horizon = 40) +plot(irf, probability = 0.68) \ No newline at end of file diff --git a/man/specify_prior_bsvarSIGN.Rd b/man/specify_prior_bsvarSIGN.Rd index d425c66..d9140e7 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, p = 1) +prior = specify_prior_bsvarSIGN$new(oil, p = 1, stationary = rep(TRUE, 3)) prior$B # show autoregressive prior mean @@ -31,6 +31,16 @@ 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, stationary = rep(TRUE, 3)) +prior$estimate_hyper(S = 5) + } \section{Public fields}{ \if{html}{\out{
}} @@ -43,17 +53,29 @@ 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{Yplus}}{an \code{(N+1)xN} matrix with the sum-of-coefficients and -dummy-initial-observation prior LHS matrix.} +\item{\code{data}}{an \code{TxN} matrix of observations.} + +\item{\code{Y}}{an \code{TxN} matrix of dependent variables.} -\item{\code{Xplus}}{an \code{(N+1)xK} matrix with the sum-of-coefficients and -dummy-initial-observation prior LHS matrix.} +\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}.} @@ -78,6 +100,7 @@ dummy-initial-observation prior LHS matrix.} \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()}} } } @@ -90,7 +113,7 @@ Create a new prior specification PriorBSVAR. \if{html}{\out{
}}\preformatted{specify_prior_bsvarSIGN$new( data, p, - d = 0, + exogenous = NULL, stationary = rep(FALSE, ncol(data)) )}\if{html}{\out{
}} } @@ -98,17 +121,15 @@ Create a new prior specification PriorBSVAR. \subsection{Arguments}{ \if{html}{\out{
}} \describe{ -\item{\code{data}}{the \code{TxN} data matrix.} +\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, otherwise to random walk.} - -\item{\code{N}}{a positive integer - the number of dependent variables in the model.} } \if{html}{\out{
}} } @@ -148,6 +169,55 @@ 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( + mu = FALSE, + delta = FALSE, + lambda = TRUE, + psi = FALSE, + S = 10000, + start = S/10, + burn = 5000 +)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\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 shrinkage in the Minnesota prior.} + +\item{\code{psi}}{whether to estimate the hyper-parameter of variances in the Minnesota prior.} + +\item{\code{S}}{number of draws.} + +\item{\code{start}}{starting point of the adaptive Metropolis algorithm.} + +\item{\code{burn}}{number of burn-in draws.} +} +\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, stationary = rep(TRUE, 3)) +prior$estimate_hyper(S = 5) + +} +\if{html}{\out{
}} + +} + } \if{html}{\out{
}} \if{html}{\out{}} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 6ed481a..e5a54f3 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 } @@ -694,49 +697,48 @@ 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& Y, const arma::mat& X); -RcppExport SEXP _bsvarSIGNs_log_ml(SEXP pSEXP, SEXP bSEXP, SEXP OmegaSEXP, SEXP PsiSEXP, SEXP dSEXP, 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 Y(YSEXP); Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP); - rcpp_result_gen = Rcpp::wrap(log_ml(p, b, Omega, Psi, d, 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); -RcppExport SEXP _bsvarSIGNs_log_ml_dummy(SEXP pSEXP, SEXP hyperSEXP, SEXP modelSEXP, SEXP YSEXP, SEXP XSEXP) { +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_result_gen = Rcpp::wrap(log_ml_dummy(p, hyper, model, Y, X)); + 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); -RcppExport SEXP _bsvarSIGNs_log_posterior_hyper(SEXP pSEXP, SEXP hyperSEXP, SEXP modelSEXP, SEXP YSEXP, SEXP XSEXP) { +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_result_gen = Rcpp::wrap(log_posterior_hyper(p, hyper, model, Y, X)); + Rcpp::traits::input_parameter< const Rcpp::List& >::type prior(priorSEXP); + rcpp_result_gen = Rcpp::wrap(log_posterior_hyper(hyper, model, Y, X, prior)); return rcpp_result_gen; END_RCPP } @@ -766,20 +768,20 @@ BEGIN_RCPP END_RCPP } // sample_hyper -arma::mat sample_hyper(const int& S, const int& start, const int& p, const arma::vec& init, const arma::vec& model, const arma::mat& Y, const arma::mat& X, const arma::mat& W); -RcppExport SEXP _bsvarSIGNs_sample_hyper(SEXP SSEXP, SEXP startSEXP, SEXP pSEXP, SEXP initSEXP, SEXP modelSEXP, SEXP YSEXP, SEXP XSEXP, SEXP WSEXP) { +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 int& >::type p(pSEXP); 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_result_gen = Rcpp::wrap(sample_hyper(S, start, p, init, model, Y, X, W)); + 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 } @@ -953,7 +955,7 @@ 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}, @@ -962,7 +964,7 @@ static const R_CallMethodDef CallEntries[] = { {"_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, 7}, + {"_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}, diff --git a/src/bsvars_sign.cpp b/src/bsvars_sign.cpp index c64b874..677bdf1 100644 --- a/src/bsvars_sign.cpp +++ b/src/bsvars_sign.cpp @@ -72,28 +72,54 @@ Rcpp::List bsvar_sign_cpp( mat B, Sigma, chol_Sigma, h_invp, Q, Epsilon, post_B, post_V, post_S; List result; + mat prior_B = as(prior["B"]), prior_V, prior_S; + int prior_nu = as(prior["nu"]); + field post; + while (s < S) { // Check for user interrupts if (s % 200 == 0) checkUserInterrupt(); + // update Minnesota prior hyper = hypers.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"]); + prior_V = diagmat(join_vert( + lambda*lambda * kron(as(prior["Vp"]), 1 / psi), + as(prior["Vd"]) + )); + prior_S = diagmat(psi); + + // update dummy observation prior + mat Ystar(0, N), Xstar(0, K); + double mu = hyper(0); + Ystar = join_vert(Ystar, as(prior["Ysoc"]) / mu); + Xstar = join_vert(Xstar, as(prior["Xsoc"]) / mu); + + double delta = hyper(1); + Ystar = join_vert(Ystar, as(prior["Ysur"]) / delta); + Xstar = join_vert(Xstar, as(prior["Xsur"]) / delta); + + mat Yplus = join_vert(Ystar, Y); + mat Xplus = join_vert(Xstar, X); + + // posterior parameters + post = niw_cpp(Yplus, Xplus, prior_B, prior_V, prior_S, prior_nu); + post_B = post(0); + post_V = post(1); + post_S = post(2); + post_nu = as_scalar(post(3)); + // sample reduced-form parameters (B, Sigma, Q) 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 + h_invp = inv(trimatl(chol_Sigma)); // lower tri, h(Sigma) is upper tri result = sample_Q(lags, Y, X, - B, h_invp, chol_Sigma, + B, h_invp, chol_Sigma, prior, VB, sign_irf, sign_narrative, sign_B, Z, max_tries); diff --git a/src/sample_hyper.cpp b/src/sample_hyper.cpp index dbaa4f7..fdac00f 100644 --- a/src/sample_hyper.cpp +++ b/src/sample_hyper.cpp @@ -122,25 +122,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); } } @@ -169,7 +175,6 @@ double log_mvgamma( // [[Rcpp:interface(cpp)]] // [[Rcpp::export]] double log_ml( - const int& p, const arma::mat& b, const arma::mat& Omega, const arma::mat& Psi, @@ -196,10 +201,6 @@ double log_ml( mat A = Psi + ehat.t() * ehat + (Bhat - b).t() * inv_Omega * (Bhat - b); log_ml += - (T + d) / 2.0 * log_det_sympd(A); - if (!std::isfinite(log_ml)) { - log_ml = -1e+10; - } - } catch(...) { log_ml = -1e+10; } @@ -212,31 +213,41 @@ 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 + const arma::mat& X, + const Rcpp::List& prior ) { - int N = Y.n_cols; - - List extended = extend_dummy(p, hyper, model, Y, X); - mat Yplus = as(extended["Yplus"]); - mat Xplus = as(extended["Xplus"]); - mat Ystar = as(extended["Ystar"]); - mat Xstar = as(extended["Xstar"]); + int N = Y.n_cols; + int K = X.n_cols; double lambda = hyper(2); vec psi = hyper.rows(3, N + 2); - List prior = mn_prior(p, lambda, psi); + mat prior_B = as(prior["B"]); - mat prior_V = as(prior["V"]); - mat prior_S = as(prior["S"]); + 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"]); - double log_ml_plus = log_ml(p, prior_B, prior_V, prior_S, prior_nu, Yplus, Xplus); - double log_ml_star = log_ml(p, prior_B, prior_V, prior_S, prior_nu, Ystar, Xstar); + mat Ystar(0, N), Xstar(0, K); + double mu = hyper(0); + Ystar = join_vert(Ystar, as(prior["Ysoc"]) / mu); + Xstar = join_vert(Xstar, as(prior["Xsoc"]) / mu); + + double delta = hyper(1); + Ystar = join_vert(Ystar, as(prior["Ysur"]) / delta); + Xstar = join_vert(Xstar, 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; } @@ -246,18 +257,22 @@ 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, - const arma::mat& X + const arma::mat& X, + const Rcpp::List& prior ) { - // double log_prior = log_prior_hyper(hyper, model, prior); - double log_prior = 0; - double log_ml = log_ml_dummy(p, hyper, model, Y, X); + double log_prior = log_prior_hyper(hyper, model, prior); + double log_ml = log_ml_dummy(hyper, model, Y, X, prior); + double log_post = log_prior + log_ml; - return log_prior + log_ml; + if (!std::isfinite(log_post)) { + log_post = -1e+10; + } + + return log_post; } @@ -333,19 +348,22 @@ arma::mat narrow_hyper( arma::mat sample_hyper( const int& S, const int& start, - const int& p, const arma::vec& init, const arma::vec& model, const arma::mat& Y, const arma::mat& X, - const arma::mat& W + const arma::mat& W, + const Rcpp::List& prior ) { mat hypers = metropolis( + S, start, narrow_hyper(model, init), W, - [p, init, model, Y, X](const vec& x) { + + [init, model, Y, X, prior](const vec& x) { + vec extended = extend_hyper(init, model, x); - return log_posterior_hyper(p, extended, model, Y, X); + return log_posterior_hyper(extended, model, Y, X, prior); } ); diff --git a/src/sample_hyper.h b/src/sample_hyper.h index 399fbb1..cf93a50 100644 --- a/src/sample_hyper.h +++ b/src/sample_hyper.h @@ -25,12 +25,12 @@ Rcpp::List extend_dummy( arma::mat sample_hyper( const int& S, const int& start, - const int& p, const arma::vec& init, const arma::vec& model, const arma::mat& Y, const arma::mat& X, - const arma::mat& W + const arma::mat& W, + const Rcpp::List& prior ); #endif // _SAMPLE_HYPER_H_ From 2064a1a44b049193b99c81c45f72d0e68f6e770d Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Sat, 13 Jul 2024 10:56:46 +1000 Subject: [PATCH 13/22] online covariance algorithm --- R/RcppExports.R | 8 ---- R/specify_bsvarSIGN.R | 1 - inst/varia/nicetry.R | 26 +------------ src/RcppExports.cpp | 30 -------------- src/bsvars_sign.cpp | 6 +-- src/sample_hyper.cpp | 91 +------------------------------------------ src/sample_hyper.h | 15 ------- src/utils.cpp | 14 +++++-- 8 files changed, 16 insertions(+), 175 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 79b1fd5..105373a 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -13,14 +13,6 @@ 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) } -mn_prior <- function(p, lambda, psi) { - .Call(`_bsvarSIGNs_mn_prior`, p, lambda, psi) -} - -extend_dummy <- function(p, hyper, model, Y, X) { - .Call(`_bsvarSIGNs_extend_dummy`, p, hyper, model, Y, X) -} - log_dgamma <- function(x, k, theta) { .Call(`_bsvarSIGNs_log_dgamma`, x, k, theta) } diff --git a/R/specify_bsvarSIGN.R b/R/specify_bsvarSIGN.R index 76c6ac5..caf6c61 100644 --- a/R/specify_bsvarSIGN.R +++ b/R/specify_bsvarSIGN.R @@ -282,7 +282,6 @@ specify_prior_bsvarSIGN = R6::R6Class( # } # Xstar = cbind(Xstar, c(rep(0, N), 1), matrix(0, N + 1, d)) - self$Y = Y self$X = X self$Vp = Vp diff --git a/inst/varia/nicetry.R b/inst/varia/nicetry.R index b23a652..63a0560 100644 --- a/inst/varia/nicetry.R +++ b/inst/varia/nicetry.R @@ -4,13 +4,13 @@ name = sapply(1:7, \(i) data$ShortDescr[[i]][[1]]) data = data$y colnames(data) = name -# data = optimism +data = optimism spec = specify_bsvarSIGN$new(data, p = 4) start_time = Sys.time() -spec$prior$estimate_hyper(TRUE, TRUE, TRUE, TRUE, S = 20000, burn = 10000, start = 2000) +spec$prior$estimate_hyper(TRUE, TRUE, TRUE, TRUE) # post = estimate(spec, S = 1000) # irf = compute_impulse_responses(post, horizon = 40) # plot(irf, probability = 0.68) @@ -20,25 +20,3 @@ end_time - start_time plot.ts(t(spec$prior$hyper)) - - - -data(optimism) - -# optimism shock -# no contemporaneous effect on productivity -zero_irf = matrix(0, nrow = 5, ncol = 5) -zero_irf[1, 1] = 1 -# positive contemporaneous effect on stock prices -sign_irf = array(0, dim = c(5, 5, 1)) -sign_irf[2, 1, 1] = 1 - -specification = specify_bsvarSIGN$new(optimism, - p = 4, - sign_irf = sign_irf, - zero_irf = zero_irf) -specification$prior$estimate_hyper(FALSE, TRUE, TRUE, TRUE) - -posterior = estimate(specification, S = 10000) -irf = compute_impulse_responses(posterior, horizon = 40) -plot(irf, probability = 0.68) \ No newline at end of file diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index e5a54f3..cc78fa2 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -617,34 +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 -} -// 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) { @@ -958,8 +930,6 @@ static const R_CallMethodDef CallEntries[] = { {"_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_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}, diff --git a/src/bsvars_sign.cpp b/src/bsvars_sign.cpp index 677bdf1..00dff97 100644 --- a/src/bsvars_sign.cpp +++ b/src/bsvars_sign.cpp @@ -86,10 +86,8 @@ Rcpp::List bsvar_sign_cpp( lambda = hyper(2); psi = hyper.rows(3, N + 2); - prior_V = diagmat(join_vert( - lambda*lambda * kron(as(prior["Vp"]), 1 / psi), - as(prior["Vd"]) - )); + prior_V = diagmat(join_vert(lambda*lambda * kron(as(prior["Vp"]), 1 / psi), + as(prior["Vd"]))); prior_S = diagmat(psi); // update dummy observation prior diff --git a/src/sample_hyper.cpp b/src/sample_hyper.cpp index fdac00f..2b8784f 100644 --- a/src/sample_hyper.cpp +++ b/src/sample_hyper.cpp @@ -7,86 +7,6 @@ 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 - ); -} - - -// 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]] @@ -194,7 +114,8 @@ double log_ml( mat ehat = Y - X * Bhat; log_ml += - N * T / 2.0 * log(M_PI); - log_ml += log_mvgamma(N, (T + d) / 2.0) - log_mvgamma(N, d / 2.0); + 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); @@ -372,11 +293,3 @@ arma::mat sample_hyper( return hypers; } - - - - - - - - diff --git a/src/sample_hyper.h b/src/sample_hyper.h index cf93a50..ed0320e 100644 --- a/src/sample_hyper.h +++ b/src/sample_hyper.h @@ -7,21 +7,6 @@ using namespace Rcpp; using namespace arma; -Rcpp::List mn_prior( - const int& p, - const double& lambda, - const arma::vec& psi -); - - -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& start, diff --git a/src/utils.cpp b/src/utils.cpp index e0e9a3d..5e1454f 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -91,11 +91,11 @@ arma::mat metropolis( int n = x.n_elem; - double s = 1.0; + double s = 2.38 / sqrt(n); double d = log_target(x); double new_d, a; - vec new_x; + vec new_x, mu, delta; mat X(n, T); x = log(x); @@ -112,9 +112,15 @@ arma::mat metropolis( } X.col(t) = x; - if (t >= t0) { + + if (t == t0) { + mu = mean(X.cols(0, t), 1); Sigma = cov(X.cols(0, t).t()); - s += pow(t, -0.6) * (a - 0.234); + } else if (t > t0) { + delta = x - mu; + s = s + pow(t, -0.6) * (a - 0.234); + mu += delta / t; + Sigma += delta * delta.t() / t; } } From 76f9e9c65dbd799ef8193accb59c358a846d3751 Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Sat, 13 Jul 2024 15:01:31 +1000 Subject: [PATCH 14/22] better covariance scaling --- R/specify_bsvarSIGN.R | 27 +++---- inst/include/bsvarSIGNs_RcppExports.h | 6 +- inst/varia/nicetry.R | 2 +- src/RcppExports.cpp | 4 +- src/bsvars_sign.cpp | 110 +++++++++++++------------- src/sample_Q.cpp | 24 +++--- src/sample_Q.h | 2 +- src/sample_hyper.cpp | 48 ++++++----- src/utils.cpp | 12 +-- 9 files changed, 114 insertions(+), 121 deletions(-) diff --git a/R/specify_bsvarSIGN.R b/R/specify_bsvarSIGN.R index caf6c61..f46b0fc 100644 --- a/R/specify_bsvarSIGN.R +++ b/R/specify_bsvarSIGN.R @@ -355,8 +355,10 @@ specify_prior_bsvarSIGN = R6::R6Class( #' prior = specify_prior_bsvarSIGN$new(oil, p = 1, stationary = rep(TRUE, 3)) #' prior$estimate_hyper(S = 5) #' - estimate_hyper = function(mu = FALSE, delta = FALSE, lambda = TRUE, psi = FALSE, - S = 10000, start = 1000, burn = 5000){ + estimate_hyper = function( + S = 10000, burn = 5000, start_adaptive = 1000, + mu = FALSE, delta = FALSE, lambda = TRUE, psi = FALSE + ) { model = c(mu, delta, lambda, psi) @@ -368,15 +370,6 @@ specify_prior_bsvarSIGN = R6::R6Class( init = narrow_hyper(model, hyper) prior = self$get_prior() - if (!mu) { - prior$Ysoc = matrix(0, 0, ncol(self$Ysoc)) - prior$Xsoc = matrix(0, 0, ncol(self$Xsoc)) - } - if (!delta) { - prior$Ysur = matrix(0, 0, ncol(self$Ysur)) - prior$Xsur = matrix(0, 0, ncol(self$Xsur)) - } - result = stats::optim( init, \(x) -log_posterior_hyper(extend_hyper(hyper, model, matrix(x)), @@ -387,17 +380,17 @@ specify_prior_bsvarSIGN = R6::R6Class( hessian = TRUE ) - mode = extend_hyper(hyper, model, matrix(result$par)) - W = result$hessian + mode = extend_hyper(hyper, model, matrix(result$par)) + variance = result$hessian if (length(init) == 1){ - W = 1 / W + variance = 1 / variance } else { - e = eigen(W) - W = e$vectors %*% diag(as.vector(1 / abs(e$values))) %*% t(e$vectors) + e = eigen(variance) + variance = e$vectors %*% diag(as.vector(1 / abs(e$values))) %*% t(e$vectors) } - self$hyper = sample_hyper(S, start, mode, model, self$Y, self$X, W, prior) + self$hyper = sample_hyper(S, start_adaptive, mode, model, self$Y, self$X, variance, prior) self$hyper = self$hyper[, -(1:burn)] } # END estimate_hyper 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/varia/nicetry.R b/inst/varia/nicetry.R index 63a0560..bd63b1e 100644 --- a/inst/varia/nicetry.R +++ b/inst/varia/nicetry.R @@ -10,7 +10,7 @@ spec = specify_bsvarSIGN$new(data, p = 4) start_time = Sys.time() -spec$prior$estimate_hyper(TRUE, TRUE, TRUE, TRUE) +spec$prior$estimate_hyper(mu = FALSE, delta = FALSE, lambda = TRUE, psi = TRUE) # post = estimate(spec, S = 1000) # irf = compute_impulse_responses(post, horizon = 40) # plot(irf, probability = 0.68) diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index cc78fa2..aef846b 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -572,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; @@ -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&)"); diff --git a/src/bsvars_sign.cpp b/src/bsvars_sign.cpp index 00dff97..ff23f76 100644 --- a/src/bsvars_sign.cpp +++ b/src/bsvars_sign.cpp @@ -52,65 +52,70 @@ Rcpp::List bsvar_sign_cpp( } 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; - 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 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); - mat hypers = as(prior["hyper"]); + mat hypers = as(prior["hyper"]); - int s = 0, S_hyper = hypers.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; - - mat prior_B = as(prior["B"]), prior_V, prior_S; + int s = 0; + int S_hyper = hypers.n_cols - 1; int prior_nu = as(prior["nu"]); - field post; + int post_nu = prior_nu + T; + + double w, mu, delta, lambda; + + vec hyper, psi; - while (s < S) { + 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; + + for (s = 0; s < S; s++) { // Check for user interrupts if (s % 200 == 0) checkUserInterrupt(); - // update Minnesota prior 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 - mat Ystar(0, N), Xstar(0, K); - double mu = hyper(0); - Ystar = join_vert(Ystar, as(prior["Ysoc"]) / mu); - Xstar = join_vert(Xstar, as(prior["Xsoc"]) / mu); - - double delta = hyper(1); - Ystar = join_vert(Ystar, as(prior["Ysur"]) / delta); - Xstar = join_vert(Xstar, as(prior["Xsur"]) / delta); - - mat Yplus = join_vert(Ystar, Y); - mat Xplus = join_vert(Xstar, X); + 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 - post = niw_cpp(Yplus, Xplus, prior_B, prior_V, prior_S, prior_nu); - post_B = post(0); - post_V = post(1); - post_S = post(2); - post_nu = as_scalar(post(3)); + 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 (B, Sigma, Q) + // sample reduced-form parameters Sigma = iwishrnd(post_S, post_nu); chol_Sigma = chol(Sigma, "lower"); B = rmatnorm_cpp(post_B, post_V, Sigma); @@ -121,24 +126,21 @@ Rcpp::List bsvar_sign_cpp( prior, VB, sign_irf, sign_narrative, sign_B, Z, max_tries); - Q = as(result["Q"]); - w = as(result["w"]); - Epsilon = as(result["Epsilon"]); + Q = result(0); + shocks = result(1); + w = as_scalar(result(2)); - 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; + // Increment progress bar + if (any(prog_rep_points == s)) p.increment(); + + 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; - s++; - } } // END s loop return List::create( 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 2b8784f..aa98fdd 100644 --- a/src/sample_hyper.cpp +++ b/src/sample_hyper.cpp @@ -115,7 +115,7 @@ double log_ml( 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 += -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); @@ -141,31 +141,27 @@ double log_ml_dummy( const Rcpp::List& prior ) { - int N = Y.n_cols; - int K = X.n_cols; - - double lambda = hyper(2); - vec psi = hyper.rows(3, N + 2); - - 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"]); - - mat Ystar(0, N), Xstar(0, K); - double mu = hyper(0); - Ystar = join_vert(Ystar, as(prior["Ysoc"]) / mu); - Xstar = join_vert(Xstar, as(prior["Xsoc"]) / mu); - - double delta = hyper(1); - Ystar = join_vert(Ystar, as(prior["Ysur"]) / delta); - Xstar = join_vert(Xstar, as(prior["Xsur"]) / delta); - - mat Yplus = join_vert(Ystar, Y); - mat Xplus = join_vert(Xstar, X); + 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); diff --git a/src/utils.cpp b/src/utils.cpp index 5e1454f..ee62846 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -95,7 +95,7 @@ arma::mat metropolis( double d = log_target(x); double new_d, a; - vec new_x, mu, delta; + vec new_x, xbar, diff; mat X(n, T); x = log(x); @@ -114,13 +114,13 @@ arma::mat metropolis( if (t == t0) { - mu = mean(X.cols(0, t), 1); + xbar = mean(X.cols(0, t), 1); Sigma = cov(X.cols(0, t).t()); } else if (t > t0) { - delta = x - mu; - s = s + pow(t, -0.6) * (a - 0.234); - mu += delta / t; - Sigma += delta * delta.t() / t; + 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); } } From 4b337999dca6876fb442f75856c6607f847a51d6 Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Sat, 13 Jul 2024 19:07:48 +1000 Subject: [PATCH 15/22] openmp trial thread number within loop does not change --- src/Makevars | 4 ++++ src/Makevars.win | 4 ++++ src/bsvars_sign.cpp | 4 ++++ 3 files changed, 12 insertions(+) diff --git a/src/Makevars b/src/Makevars index 22c7566..075f362 100644 --- a/src/Makevars +++ b/src/Makevars @@ -17,3 +17,7 @@ PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) + +PKG_CFLAGS += -Xclang -fopenmp +PKG_FFLAGS += -Xclang -fopenmp +PKG_LIBS += -lomp \ No newline at end of file diff --git a/src/Makevars.win b/src/Makevars.win index 22c7566..075f362 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -17,3 +17,7 @@ PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) + +PKG_CFLAGS += -Xclang -fopenmp +PKG_FFLAGS += -Xclang -fopenmp +PKG_LIBS += -lomp \ No newline at end of file diff --git a/src/bsvars_sign.cpp b/src/bsvars_sign.cpp index ff23f76..fae7c7f 100644 --- a/src/bsvars_sign.cpp +++ b/src/bsvars_sign.cpp @@ -1,4 +1,5 @@ +#include #include #include "progress.hpp" #include "Rcpp/Rmath.h" @@ -86,7 +87,10 @@ Rcpp::List bsvar_sign_cpp( field result; + #pragma omp parallel for for (s = 0; s < S; s++) { + + cout << "Thread number " << omp_get_thread_num() << endl; // Check for user interrupts if (s % 200 == 0) checkUserInterrupt(); From ae2d5c32b9b3f618aab908ce4f5dbd530acaa0cf Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Sat, 13 Jul 2024 22:56:31 +1000 Subject: [PATCH 16/22] fix openmp but disables progress bar #13 --- man/specify_prior_bsvarSIGN.Rd | 16 ++++++++-------- src/Makevars | 4 +--- src/Makevars.win | 6 +----- src/bsvars_sign.cpp | 16 +++++----------- 4 files changed, 15 insertions(+), 27 deletions(-) diff --git a/man/specify_prior_bsvarSIGN.Rd b/man/specify_prior_bsvarSIGN.Rd index d9140e7..98ace19 100644 --- a/man/specify_prior_bsvarSIGN.Rd +++ b/man/specify_prior_bsvarSIGN.Rd @@ -177,19 +177,23 @@ prior$get_prior() # show the prior as list Estimates hyper-parameters with adaptive Metropolis algorithm. \subsection{Usage}{ \if{html}{\out{
}}\preformatted{specify_prior_bsvarSIGN$estimate_hyper( + S = 10000, + burn = 5000, + start_adaptive = 1000, mu = FALSE, delta = FALSE, lambda = TRUE, - psi = FALSE, - S = 10000, - start = S/10, - burn = 5000 + psi = FALSE )}\if{html}{\out{
}} } \subsection{Arguments}{ \if{html}{\out{
}} \describe{ +\item{\code{S}}{number of draws.} + +\item{\code{burn}}{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.} @@ -198,11 +202,7 @@ Estimates hyper-parameters with adaptive Metropolis algorithm. \item{\code{psi}}{whether to estimate the hyper-parameter of variances in the Minnesota prior.} -\item{\code{S}}{number of draws.} - \item{\code{start}}{starting point of the adaptive Metropolis algorithm.} - -\item{\code{burn}}{number of burn-in draws.} } \if{html}{\out{
}} } diff --git a/src/Makevars b/src/Makevars index 075f362..bd914ee 100644 --- a/src/Makevars +++ b/src/Makevars @@ -18,6 +18,4 @@ PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) -PKG_CFLAGS += -Xclang -fopenmp -PKG_FFLAGS += -Xclang -fopenmp -PKG_LIBS += -lomp \ No newline at end of file +PKG_CXXFLAGS += -Xclang -fopenmp \ No newline at end of file diff --git a/src/Makevars.win b/src/Makevars.win index 075f362..48f15f8 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -16,8 +16,4 @@ #CXX_STD = CXX11 PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) -PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) - -PKG_CFLAGS += -Xclang -fopenmp -PKG_FFLAGS += -Xclang -fopenmp -PKG_LIBS += -lomp \ No newline at end of file +PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) \ No newline at end of file diff --git a/src/bsvars_sign.cpp b/src/bsvars_sign.cpp index fae7c7f..b977cb3 100644 --- a/src/bsvars_sign.cpp +++ b/src/bsvars_sign.cpp @@ -1,10 +1,9 @@ -#include #include #include "progress.hpp" #include "Rcpp/Rmath.h" - #include +#include #include "sample_hyper.h" #include "sample_Q.h" @@ -87,13 +86,11 @@ Rcpp::List bsvar_sign_cpp( field result; - #pragma omp parallel for + #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 (s = 0; s < S; s++) { - cout << "Thread number " << omp_get_thread_num() << endl; - // Check for user interrupts - if (s % 200 == 0) checkUserInterrupt(); + // if (s % 200 == 0) checkUserInterrupt(); hyper = hypers.col(randi(distr_param(0, S_hyper))); mu = hyper(0); @@ -111,7 +108,7 @@ Rcpp::List bsvar_sign_cpp( 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); @@ -135,7 +132,7 @@ Rcpp::List bsvar_sign_cpp( w = as_scalar(result(2)); // Increment progress bar - if (any(prog_rep_points == s)) p.increment(); + // if (any(prog_rep_points == s)) p.increment(); posterior_w(s) = w; posterior_hyper.col(s) = hyper; @@ -144,7 +141,6 @@ Rcpp::List bsvar_sign_cpp( posterior_Sigma.slice(s) = Sigma; posterior_Theta0.slice(s) = chol_Sigma * Q; posterior_shocks.slice(s) = shocks; - } // END s loop return List::create( @@ -160,5 +156,3 @@ Rcpp::List bsvar_sign_cpp( ); } // END bsvar_sign_cpp - - From c2e95bfdbd267e0de66beec751348722c57d209a Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Sat, 13 Jul 2024 23:34:03 +1000 Subject: [PATCH 17/22] progress bar in parallel --- src/bsvars_sign.cpp | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/src/bsvars_sign.cpp b/src/bsvars_sign.cpp index b977cb3..23fd06f 100644 --- a/src/bsvars_sign.cpp +++ b/src/bsvars_sign.cpp @@ -38,7 +38,12 @@ Rcpp::List bsvar_sign_cpp( } // Progress bar setup - vec prog_rep_points = arma::round(arma::linspace(0, S, 50)); + double num_threads; + #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; @@ -89,9 +94,6 @@ Rcpp::List bsvar_sign_cpp( #pragma omp parallel for private(hyper, mu, delta, lambda, psi, prior_V, prior_S, Ystar, Xstar, Yplus, Xplus, result, post_B, post_V, post_S, Sigma, chol_Sigma, B, h_invp, Q, shocks, w) for (s = 0; s < S; s++) { - // Check for user interrupts - // if (s % 200 == 0) checkUserInterrupt(); - hyper = hypers.col(randi(distr_param(0, S_hyper))); mu = hyper(0); delta = hyper(1); @@ -131,9 +133,6 @@ Rcpp::List bsvar_sign_cpp( shocks = result(1); w = as_scalar(result(2)); - // Increment progress bar - // if (any(prog_rep_points == s)) p.increment(); - posterior_w(s) = w; posterior_hyper.col(s) = hyper; posterior_A.slice(s) = B.t(); @@ -141,6 +140,14 @@ Rcpp::List bsvar_sign_cpp( posterior_Sigma.slice(s) = Sigma; posterior_Theta0.slice(s) = chol_Sigma * Q; posterior_shocks.slice(s) = shocks; + + 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(); + } } // END s loop return List::create( From dd40b4cad45f9421eb6171fdaa7f3006de19a241 Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Sun, 14 Jul 2024 00:29:13 +1000 Subject: [PATCH 18/22] force valid draw with while loop --- R/estimate.BSVARSIGN.R | 3 +- src/Makevars | 5 ++- src/bsvars_sign.cpp | 81 ++++++++++++++++++++++-------------------- 3 files changed, 48 insertions(+), 41 deletions(-) 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/src/Makevars b/src/Makevars index bd914ee..cca1220 100644 --- a/src/Makevars +++ b/src/Makevars @@ -18,4 +18,7 @@ PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) -PKG_CXXFLAGS += -Xclang -fopenmp \ No newline at end of file +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/bsvars_sign.cpp b/src/bsvars_sign.cpp index 23fd06f..3f35861 100644 --- a/src/bsvars_sign.cpp +++ b/src/bsvars_sign.cpp @@ -94,44 +94,47 @@ Rcpp::List bsvar_sign_cpp( #pragma omp parallel for private(hyper, mu, delta, lambda, psi, prior_V, prior_S, Ystar, Xstar, Yplus, Xplus, result, post_B, post_V, post_S, Sigma, chol_Sigma, B, h_invp, Q, shocks, w) for (s = 0; s < S; s++) { - hyper = hypers.col(randi(distr_param(0, S_hyper))); - mu = hyper(0); - delta = hyper(1); - lambda = hyper(2); - psi = hyper.rows(3, N + 2); - - // update Minnesota prior - prior_V = diagmat(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)); + 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)); + } posterior_w(s) = w; posterior_hyper.col(s) = hyper; @@ -144,7 +147,7 @@ Rcpp::List bsvar_sign_cpp( 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(); } From 68558f53425bd190825459291aeb421cffaa3668 Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Sun, 14 Jul 2024 15:47:49 +1000 Subject: [PATCH 19/22] mcmc progress bar --- R/specify_bsvarSIGN.R | 23 +++++++++++++---------- README.Rmd | 8 ++++---- README.md | 9 ++++----- inst/varia/nicetry.R | 4 +++- man/specify_prior_bsvarSIGN.Rd | 21 +++++++++++---------- src/bsvars_sign.cpp | 14 +++++++------- src/utils.cpp | 22 +++++++++++++++------- 7 files changed, 57 insertions(+), 44 deletions(-) diff --git a/R/specify_bsvarSIGN.R b/R/specify_bsvarSIGN.R index f46b0fc..3584d01 100644 --- a/R/specify_bsvarSIGN.R +++ b/R/specify_bsvarSIGN.R @@ -341,13 +341,16 @@ specify_prior_bsvarSIGN = R6::R6Class( #' @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 shrinkage in the Minnesota prior. - #' @param psi whether to estimate the hyper-parameter of variances in the Minnesota prior. - #' @param S number of draws. - #' @param start starting point of the adaptive Metropolis algorithm. - #' @param burn number of burn-in draws. + #' @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 @@ -356,7 +359,7 @@ specify_prior_bsvarSIGN = R6::R6Class( #' prior$estimate_hyper(S = 5) #' estimate_hyper = function( - S = 10000, burn = 5000, start_adaptive = 1000, + S = 10000, burn_in = S / 2, mu = FALSE, delta = FALSE, lambda = TRUE, psi = FALSE ) { @@ -390,8 +393,8 @@ specify_prior_bsvarSIGN = R6::R6Class( variance = e$vectors %*% diag(as.vector(1 / abs(e$values))) %*% t(e$vectors) } - self$hyper = sample_hyper(S, start_adaptive, mode, model, self$Y, self$X, variance, prior) - self$hyper = self$hyper[, -(1:burn)] + 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 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/varia/nicetry.R b/inst/varia/nicetry.R index bd63b1e..6b428e9 100644 --- a/inst/varia/nicetry.R +++ b/inst/varia/nicetry.R @@ -10,7 +10,8 @@ spec = specify_bsvarSIGN$new(data, p = 4) start_time = Sys.time() -spec$prior$estimate_hyper(mu = FALSE, delta = FALSE, lambda = TRUE, psi = TRUE) +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) @@ -20,3 +21,4 @@ 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 98ace19..5bb1d2e 100644 --- a/man/specify_prior_bsvarSIGN.Rd +++ b/man/specify_prior_bsvarSIGN.Rd @@ -178,8 +178,7 @@ Estimates hyper-parameters with adaptive Metropolis algorithm. \subsection{Usage}{ \if{html}{\out{
}}\preformatted{specify_prior_bsvarSIGN$estimate_hyper( S = 10000, - burn = 5000, - start_adaptive = 1000, + burn_in = S/2, mu = FALSE, delta = FALSE, lambda = TRUE, @@ -190,19 +189,21 @@ Estimates hyper-parameters with adaptive Metropolis algorithm. \subsection{Arguments}{ \if{html}{\out{
}} \describe{ -\item{\code{S}}{number of draws.} +\item{\code{S}}{number of MCMC draws.} -\item{\code{burn}}{number of burn-in 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{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{delta}}{whether to estimate the hyper-parameter in the +single-unit-root dummy prior.} -\item{\code{lambda}}{whether to estimate the hyper-parameter of shrinkage in the Minnesota 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 variances in the Minnesota prior.} - -\item{\code{start}}{starting point of the adaptive Metropolis algorithm.} +\item{\code{psi}}{whether to estimate the hyper-parameter of the +variances in the Minnesota prior.} } \if{html}{\out{
}} } diff --git a/src/bsvars_sign.cpp b/src/bsvars_sign.cpp index 3f35861..ef84863 100644 --- a/src/bsvars_sign.cpp +++ b/src/bsvars_sign.cpp @@ -46,12 +46,13 @@ Rcpp::List bsvar_sign_cpp( 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 << " bsvarSIGNs: Bayesian Structural VAR with zero, |" << endl; + Rcout << " sign and narrative restriction |" << endl; Rcout << "**************************************************|" << endl; - Rcout << " Gibbs sampler for the SVAR model |" << 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; } @@ -71,7 +72,6 @@ Rcpp::List bsvar_sign_cpp( mat hypers = as(prior["hyper"]); - int s = 0; int S_hyper = hypers.n_cols - 1; int prior_nu = as(prior["nu"]); int post_nu = prior_nu + T; @@ -92,7 +92,7 @@ Rcpp::List bsvar_sign_cpp( 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 (s = 0; s < S; s++) { + for (int s = 0; s < S; s++) { w = 0; while (w == 0) { diff --git a/src/utils.cpp b/src/utils.cpp index ee62846..4760fcc 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -1,8 +1,9 @@ #include #include - #include +#include "Rcpp/Rmath.h" +#include "progress.hpp" using namespace arma; @@ -91,20 +92,23 @@ arma::mat metropolis( int n = x.n_elem; + x = log(x); double s = 2.38 / sqrt(n); - double d = log_target(x); + double d = log_target(x) + sum(x); double new_d, a; vec new_x, xbar, diff; mat X(n, T); - x = log(x); 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)); - a = std::min(1.0, exp(new_d - d + sum(x - new_x))); + 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; @@ -113,15 +117,19 @@ arma::mat metropolis( X.col(t) = x; - if (t == t0) { + if (t == 1) { xbar = mean(X.cols(0, t), 1); Sigma = cov(X.cols(0, t).t()); - } else if (t > t0) { + } else 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); From eb3f587c041bf8fa9164f27257150f039b839f56 Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Mon, 15 Jul 2024 14:26:18 +1000 Subject: [PATCH 20/22] disable openmp for now --- R/specify_bsvarSIGN.R | 8 +++---- R/utils.R | 2 -- inst/tinytest/test_estimate.R | 10 ++++---- man/specify_prior_bsvarSIGN.Rd | 12 +++++----- src/Makevars | 8 +++---- src/bsvars_sign.cpp | 42 +++++++++++++++++++--------------- 6 files changed, 42 insertions(+), 40 deletions(-) diff --git a/R/specify_bsvarSIGN.R b/R/specify_bsvarSIGN.R index 3584d01..cd04204 100644 --- a/R/specify_bsvarSIGN.R +++ b/R/specify_bsvarSIGN.R @@ -131,7 +131,7 @@ igamma_shape = function(mode, variance) { #' @examples #' # a prior for 3-variable example with one lag #' data(oil) -#' prior = specify_prior_bsvarSIGN$new(oil, p = 1, stationary = rep(TRUE, 3)) +#' prior = specify_prior_bsvarSIGN$new(oil, p = 1) #' prior$B # show autoregressive prior mean #' #' @export @@ -221,10 +221,10 @@ 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, p = 1, stationary = rep(TRUE, 3)) + #' prior = specify_prior_bsvarSIGN$new(oil, p = 1) #' prior$B # show autoregressive prior mean #' - initialize = function(data, p, exogenous = NULL, stationary = rep(FALSE, ncol(data))) { + 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) @@ -355,7 +355,7 @@ specify_prior_bsvarSIGN = R6::R6Class( #' @examples #' # a prior for 3-variable example with four lags #' data(oil) - #' prior = specify_prior_bsvarSIGN$new(oil, p = 1, stationary = rep(TRUE, 3)) + #' prior = specify_prior_bsvarSIGN$new(oil, p = 1) #' prior$estimate_hyper(S = 5) #' estimate_hyper = function( 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/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/man/specify_prior_bsvarSIGN.Rd b/man/specify_prior_bsvarSIGN.Rd index 5bb1d2e..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, p = 1, stationary = rep(TRUE, 3)) +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, p = 1, stationary = rep(TRUE, 3)) +prior = specify_prior_bsvarSIGN$new(oil, p = 1) prior$B # show autoregressive prior mean @@ -38,7 +38,7 @@ prior$get_prior() # show the prior as list # a prior for 3-variable example with four lags data(oil) -prior = specify_prior_bsvarSIGN$new(oil, p = 1, stationary = rep(TRUE, 3)) +prior = specify_prior_bsvarSIGN$new(oil, p = 1) prior$estimate_hyper(S = 5) } @@ -114,7 +114,7 @@ Create a new prior specification PriorBSVAR. data, p, exogenous = NULL, - stationary = rep(FALSE, ncol(data)) + stationary = rep(FALSE, dim(data)[2]) )}\if{html}{\out{
}} } @@ -140,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, p = 1, stationary = rep(TRUE, 3)) +prior = specify_prior_bsvarSIGN$new(oil, p = 1) prior$B # show autoregressive prior mean } @@ -211,7 +211,7 @@ variances in the Minnesota prior.} \if{html}{\out{
}} \preformatted{# a prior for 3-variable example with four lags data(oil) -prior = specify_prior_bsvarSIGN$new(oil, p = 1, stationary = rep(TRUE, 3)) +prior = specify_prior_bsvarSIGN$new(oil, p = 1) prior$estimate_hyper(S = 5) } diff --git a/src/Makevars b/src/Makevars index cca1220..62a255e 100644 --- a/src/Makevars +++ b/src/Makevars @@ -18,7 +18,7 @@ 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 +#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/bsvars_sign.cpp b/src/bsvars_sign.cpp index ef84863..ead6630 100644 --- a/src/bsvars_sign.cpp +++ b/src/bsvars_sign.cpp @@ -3,7 +3,7 @@ #include "progress.hpp" #include "Rcpp/Rmath.h" #include -#include +// #include #include "sample_hyper.h" #include "sample_Q.h" @@ -38,11 +38,11 @@ Rcpp::List bsvar_sign_cpp( } // Progress bar setup - double num_threads; - #pragma omp parallel - { - num_threads = omp_get_num_threads(); - } + 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; @@ -91,7 +91,7 @@ Rcpp::List bsvar_sign_cpp( 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) + // #pragma omp parallel for private(hyper, mu, delta, lambda, psi, prior_V, prior_S, Ystar, Xstar, Yplus, Xplus, result, post_B, post_V, post_S, Sigma, chol_Sigma, B, h_invp, Q, shocks, w) for (int s = 0; s < S; s++) { w = 0; @@ -126,11 +126,8 @@ Rcpp::List bsvar_sign_cpp( 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); + 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)); @@ -144,13 +141,20 @@ Rcpp::List bsvar_sign_cpp( posterior_Theta0.slice(s) = chol_Sigma * Q; posterior_shocks.slice(s) = shocks; - 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(); - } + // Check for user interrupts + if (s % 200 == 0) checkUserInterrupt(); + + // 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(); + // } + } // END s loop return List::create( From 9631e39ed76eca2441203f548d889b1795f6444e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tomasz=20Wo=C5=BAniak?= Date: Mon, 15 Jul 2024 15:20:14 +1000 Subject: [PATCH 21/22] fixed the bug wth initialisation of specify_prior #13 --- DESCRIPTION | 2 +- R/specify_bsvarSIGN.R | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index b8d120a..1c5edaa 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -22,5 +22,5 @@ Depends: bsvars Encoding: UTF-8 LazyData: true -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.2 Suggests: tinytest diff --git a/R/specify_bsvarSIGN.R b/R/specify_bsvarSIGN.R index cd04204..933f4b6 100644 --- a/R/specify_bsvarSIGN.R +++ b/R/specify_bsvarSIGN.R @@ -228,9 +228,9 @@ specify_prior_bsvarSIGN = R6::R6Class( stopifnot("Argument p must be a positive integer number." = p > 0 & p %% 1 == 0) - data = bsvars::specify_data_matrices$new(data, p, exogenous) - Y = t(data$Y) - X = t(data$X) + 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)) From a05aeff716952d89c2cf6da0c32c62b4ec518075 Mon Sep 17 00:00:00 2001 From: adamwang15 Date: Mon, 15 Jul 2024 15:28:09 +1000 Subject: [PATCH 22/22] import ramcmc --- DESCRIPTION | 3 +- src/mcmc.cpp | 70 ++++++++++++++++++++++++++++++++++++++++++++ src/mcmc.h | 23 +++++++++++++++ src/sample_hyper.cpp | 1 + src/utils.cpp | 65 ---------------------------------------- 5 files changed, 96 insertions(+), 66 deletions(-) create mode 100644 src/mcmc.cpp create mode 100644 src/mcmc.h 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/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_hyper.cpp b/src/sample_hyper.cpp index aa98fdd..9e8ee58 100644 --- a/src/sample_hyper.cpp +++ b/src/sample_hyper.cpp @@ -2,6 +2,7 @@ #include #include "utils.h" +#include "mcmc.h" using namespace Rcpp; using namespace arma; diff --git a/src/utils.cpp b/src/utils.cpp index 4760fcc..aaaad02 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -2,8 +2,6 @@ #include #include #include -#include "Rcpp/Rmath.h" -#include "progress.hpp" using namespace arma; @@ -79,66 +77,3 @@ arma::mat Df( return result; } - -// 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); - double s = 2.38 / sqrt(n); - double d = log_target(x) + sum(x); - - double new_d, a; - vec new_x, xbar, 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 == 1) { - xbar = mean(X.cols(0, t), 1); - Sigma = cov(X.cols(0, t).t()); - } else 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; -} - - - - - - -