hyperparameters in prior

bsvars · Jun 17, 2024 · 8fc1b68 · 8fc1b68
1 parent 9f313be
commit 8fc1b68
Show file tree

Hide file tree

Showing 11 changed files with 238 additions and 305 deletions.
diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -13,12 +13,16 @@ niw_cpp <- function(Y, X, prior) {
     .Call(`_bsvarSIGNs_niw_cpp`, Y, X, prior)
 }
 
-update_prior <- function(p, hyper, prior) {
-    .Call(`_bsvarSIGNs_update_prior`, p, hyper, prior)
+mn_prior <- function(p, lambda, psi) {
+    .Call(`_bsvarSIGNs_mn_prior`, p, lambda, psi)
 }
 
-extend_dummy <- function(p, hyper, Y, X) {
-    .Call(`_bsvarSIGNs_extend_dummy`, p, hyper, Y, X)
+update_prior <- function(p, hyper, model, prior) {
+    .Call(`_bsvarSIGNs_update_prior`, p, hyper, model, prior)
+}
+
+extend_dummy <- function(p, hyper, model, Y, X) {
+    .Call(`_bsvarSIGNs_extend_dummy`, p, hyper, model, Y, X)
 }
 
 log_dgamma <- function(x, k, theta) {
@@ -29,8 +33,8 @@ log_dinvgamma <- function(x, alpha, beta) {
     .Call(`_bsvarSIGNs_log_dinvgamma`, x, alpha, beta)
 }
 
-log_prior_hyper <- function(hyper, prior) {
-    .Call(`_bsvarSIGNs_log_prior_hyper`, hyper, prior)
+log_prior_hyper <- function(hyper, model, prior) {
+    .Call(`_bsvarSIGNs_log_prior_hyper`, hyper, model, prior)
 }
 
 log_mvgamma <- function(n, x) {
@@ -41,16 +45,12 @@ log_ml <- function(p, b, Omega, Psi, d, inv_Omega, Y, X) {
     .Call(`_bsvarSIGNs_log_ml`, p, b, Omega, Psi, d, inv_Omega, Y, X)
 }
 
-log_ml_dummy <- function(p, hyper, Y, X, prior) {
-    .Call(`_bsvarSIGNs_log_ml_dummy`, p, hyper, Y, X, prior)
-}
-
-log_posterior_hyper <- function(p, hyper, Y, X, prior) {
-    .Call(`_bsvarSIGNs_log_posterior_hyper`, p, hyper, Y, X, prior)
+log_ml_dummy <- function(p, hyper, model, Y, X, prior) {
+    .Call(`_bsvarSIGNs_log_ml_dummy`, p, hyper, model, Y, X, prior)
 }
 
-sample_hyper <- function(S, p, c, Y, X, prior) {
-    .Call(`_bsvarSIGNs_sample_hyper`, S, p, c, Y, X, prior)
+log_posterior_hyper <- function(p, hyper, model, Y, X, prior) {
+    .Call(`_bsvarSIGNs_log_posterior_hyper`, p, hyper, model, Y, X, prior)
 }
 
 # Register entry points for exported C++ functions

diff --git a/R/estimate.BSVARSIGN.R b/R/estimate.BSVARSIGN.R
@@ -108,7 +108,7 @@ estimate.BSVARSIGN <- function(specification, S, thin = 1, show_progress = TRUE)
                               data_matrices$Y, data_matrices$X, identification$VB, 
                               identification$sign_irf, identification$sign_narrative, 
                               identification$sign_B, identification$zero_irf,
-                              prior, starting_values, thin, show_progress, max_tries)
+                              prior, starting_values, show_progress, thin, max_tries)
 
   specification$starting_values$set_starting_values(qqq$last_draw)
   output              = specify_posterior_bsvarSIGN$new(specification, qqq$posterior)

diff --git a/R/sample_hyper.R b/R/sample_hyper.R
diff --git a/R/specify_bsvarSIGN.R b/R/specify_bsvarSIGN.R
@@ -124,21 +124,12 @@ specify_prior_bsvarSIGN = R6::R6Class(
   "PriorBSVARSIGN",
 
   public = list(
+    #' @field hyper a \code{(N+3)xS} matrix of posterior draws of hyperparameters.
+    hyper      = matrix(),
 
-    #' @field B an \code{KxN} matrix, the mean of the normal prior distribution for the parameter matrix \eqn{B}. 
-    B          = matrix(),
-
-    #' @field V a \code{KxK} covariance matrix of the normal prior distribution for each of 
-    #' the column of the parameter matrix \eqn{B}. This covariance matrix is equation invariant.
-    V          = matrix(),
-
-    #' @field S an \code{NxN} scale matrix of the inverse-Wishart prior distribution 
-    #' for the covariance matrix \eqn{\Sigma}. This scale matrix is equation invariant.
-    S          = matrix(),
-
-    #' @field nu a positive real number greater of equal than \code{N}, a degree of freedom parameter
-    #' of the inverse-Wishart prior distribution for the covariance matrix \eqn{\Sigma}.
-    nu         = NA,
+    #' @field model a \code{4x1} vector of Boolean values indicating prior specifications,
+    #' model[1] = dummy soc, model[2] = dummy sur, model[3] = mn lambda, model[4] = mn psi.
+    model      = c(),
 
     #' @description
     #' Create a new prior specification PriorBSVAR.
@@ -162,11 +153,17 @@ specify_prior_bsvarSIGN = R6::R6Class(
       stopifnot("Argument d must be a non-negative integer number." = d >= 0 & d %% 1 == 0)
       stopifnot("Argument stationary must be a logical vector of length N." = length(stationary) == N & is.logical(stationary))
 
-      prior   = niw_prior(data, p, !stationary)
-      self$B  = prior$B
-      self$V  = prior$V
-      self$S  = prior$S
-      self$nu = prior$nu
+      Y      = data
+      T      = nrow(Y)
+      lambda = 0.2
+      psi    = sapply(1:N, \(i) summary(lm(Y[2:T, i] ~ Y[1:(T - 1), i]))$sigma^2)
+
+      hyper              = matrix(NA, N + 3, 1)
+      hyper[3, ]         = lambda
+      hyper[4:(N + 3), ] = psi
+
+      self$model         = c(FALSE, FALSE, FALSE, FALSE)
+      self$hyper         = hyper
 
     }, # END initialize
 
@@ -180,10 +177,8 @@ specify_prior_bsvarSIGN = R6::R6Class(
     #' 
     get_prior = function(){
       list(
-        B  = self$B,
-        V  = self$V,
-        S  = self$S,
-        nu = self$nu
+        model = self$model,
+        hyper = self$hyper
       )
     } # END get_prior
 

diff --git a/R/utils.R b/R/utils.R
@@ -13,7 +13,7 @@ importance_sampling <- function(posterior) {
 
   posterior$posterior$A      = posterior$posterior$A[, , indices]
   posterior$posterior$B      = posterior$posterior$B[, , indices]
-  posterior$posterior$hyper  = posterior$posterior$hyper[, , indices]
+  posterior$posterior$hyper  = posterior$posterior$hyper[, indices]
   posterior$posterior$Q      = posterior$posterior$Q[, , indices]
   posterior$posterior$Sigma  = posterior$posterior$Sigma[, , indices]
   posterior$posterior$Theta0 = posterior$posterior$Theta0[, , indices]

diff --git a/inst/include/bsvarSIGNs_RcppExports.h b/inst/include/bsvarSIGNs_RcppExports.h
@@ -88,17 +88,17 @@ namespace bsvarSIGNs {
         return Rcpp::as<arma::field<arma::cube> >(rcpp_result_gen);
     }
 
-    inline Rcpp::List bsvar_sign_cpp(const int& S, const int& lags, const arma::mat& Y, const arma::mat& X, const arma::field<arma::mat>& VB, const arma::cube& sign_irf, const arma::mat& sign_narrative, const arma::mat& sign_B, const arma::field<arma::mat>& Z, const Rcpp::List& prior, const Rcpp::List& starting_values, const int thin = 100, const bool show_progress = true, const int& max_tries = 10000) {
+    inline Rcpp::List bsvar_sign_cpp(const int& S, const int& lags, const arma::mat& Y, const arma::mat& X, const arma::field<arma::mat>& VB, const arma::cube& sign_irf, const arma::mat& sign_narrative, const arma::mat& sign_B, const arma::field<arma::mat>& Z, const Rcpp::List& prior, const Rcpp::List& starting_values, const bool show_progress = true, const int thin = 100, const int& max_tries = 10000) {
         typedef SEXP(*Ptr_bsvar_sign_cpp)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP);
         static Ptr_bsvar_sign_cpp p_bsvar_sign_cpp = NULL;
         if (p_bsvar_sign_cpp == NULL) {
-            validateSignature("Rcpp::List(*bsvar_sign_cpp)(const int&,const int&,const arma::mat&,const arma::mat&,const arma::field<arma::mat>&,const arma::cube&,const arma::mat&,const arma::mat&,const arma::field<arma::mat>&,const Rcpp::List&,const Rcpp::List&,const int,const bool,const int&)");
+            validateSignature("Rcpp::List(*bsvar_sign_cpp)(const int&,const int&,const arma::mat&,const arma::mat&,const arma::field<arma::mat>&,const arma::cube&,const arma::mat&,const arma::mat&,const arma::field<arma::mat>&,const Rcpp::List&,const Rcpp::List&,const bool,const int,const int&)");
             p_bsvar_sign_cpp = (Ptr_bsvar_sign_cpp)R_GetCCallable("bsvarSIGNs", "_bsvarSIGNs_bsvar_sign_cpp");
         }
         RObject rcpp_result_gen;
         {
             RNGScope RCPP_rngScope_gen;
-            rcpp_result_gen = p_bsvar_sign_cpp(Shield<SEXP>(Rcpp::wrap(S)), Shield<SEXP>(Rcpp::wrap(lags)), Shield<SEXP>(Rcpp::wrap(Y)), Shield<SEXP>(Rcpp::wrap(X)), Shield<SEXP>(Rcpp::wrap(VB)), Shield<SEXP>(Rcpp::wrap(sign_irf)), Shield<SEXP>(Rcpp::wrap(sign_narrative)), Shield<SEXP>(Rcpp::wrap(sign_B)), Shield<SEXP>(Rcpp::wrap(Z)), Shield<SEXP>(Rcpp::wrap(prior)), Shield<SEXP>(Rcpp::wrap(starting_values)), Shield<SEXP>(Rcpp::wrap(thin)), Shield<SEXP>(Rcpp::wrap(show_progress)), Shield<SEXP>(Rcpp::wrap(max_tries)));
+            rcpp_result_gen = p_bsvar_sign_cpp(Shield<SEXP>(Rcpp::wrap(S)), Shield<SEXP>(Rcpp::wrap(lags)), Shield<SEXP>(Rcpp::wrap(Y)), Shield<SEXP>(Rcpp::wrap(X)), Shield<SEXP>(Rcpp::wrap(VB)), Shield<SEXP>(Rcpp::wrap(sign_irf)), Shield<SEXP>(Rcpp::wrap(sign_narrative)), Shield<SEXP>(Rcpp::wrap(sign_B)), Shield<SEXP>(Rcpp::wrap(Z)), Shield<SEXP>(Rcpp::wrap(prior)), Shield<SEXP>(Rcpp::wrap(starting_values)), Shield<SEXP>(Rcpp::wrap(show_progress)), Shield<SEXP>(Rcpp::wrap(thin)), Shield<SEXP>(Rcpp::wrap(max_tries)));
         }
         if (rcpp_result_gen.inherits("interrupted-error"))
             throw Rcpp::internal::InterruptedException();