Modify the functions for the Villani Steady-State prior to accept a m…

…atrix instead of just a double. This way it is possible to set an individual prior variance per variable, allowing the user to express how sure/unsure he is about the steady state.
jonlachmann · Feb 22, 2022 · e54e00b · e54e00b
1 parent 63565cd
commit e54e00b
Show file tree

Hide file tree

Showing 6 changed files with 14 additions and 14 deletions.
diff --git a/man/bvars.Rd b/man/bvars.Rd
@@ -39,15 +39,15 @@ psi_prior <- matrix(c(3,6,5),nrow=1)
 
 HP1 <- 0.5
 HP2 <- 0.5
-XiPsi <- 1
+psi_prior_var <- diag(3)
 gamma = 4
 
 bvar_obj <- new(bvars)
 
 #
 
 bvar_obj$build(bvar_data,TRUE,4)
-bvar_obj$prior(coef_prior,HP1,HP2,psi_prior,XiPsi,gamma)
+bvar_obj$prior(coef_prior,HP1,HP2,psi_prior,psi_prior_var,gamma)
 bvar_obj$gibbs(10000,5000)
 
 IRF(bvar_obj,20,var_names=colnames(USMacroData),save=FALSE)

diff --git a/src/bmlib/include/var/bvars.hpp b/src/bmlib/include/var/bvars.hpp
@@ -100,7 +100,7 @@ class bvars
         void reset_draws();
 
         void prior(const arma::vec& coef_prior, const double HP_1, const double HP_2, 
-                   const arma::mat& Psi_prior, const double Xi_psi, const int gamma, 
+                   const arma::mat& Psi_prior, const arma::mat& Psi_prior_var, const int gamma,
                    const bool full_cov_prior = true);
 
         void gibbs(const uint_t n_draws, const uint_t n_burnin);

diff --git a/src/bmlib/src/var/bvars.cpp b/src/bmlib/src/var/bvars.cpp
@@ -151,7 +151,7 @@ bm::bvars::minn_pr_var()
 
 void
 bm::bvars::prior(const arma::vec& coef_prior, const double HP_1, const double HP_2, 
-                 const arma::mat& Psi_prior, const double Xi_psi, const int gamma,
+                 const arma::mat& Psi_prior, const arma::mat& Psi_prior_var, const int gamma,
                  const bool full_cov_prior)
 {
     arma::mat Z = arma::join_rows(X,d);
@@ -175,7 +175,7 @@ bm::bvars::prior(const arma::vec& coef_prior, const double HP_1, const double HP
     //
 
     psi_pr_mean = arma::vectorise(Psi_prior.t());
-    psi_pr_var = Xi_psi*arma::eye(q*M,q*M);
+    psi_pr_var = Psi_prior_var;
 
     //
 

diff --git a/src/modules/bvars_R.cpp b/src/modules/bvars_R.cpp
@@ -30,7 +30,7 @@ RCPP_MODULE(bvars_module)
     void (bvars_R::*build_1)(const arma::mat&, bool, int) = &bvars_R::build_R;
     void (bvars_R::*build_2)(const arma::mat&, const arma::mat&, bool, int) = &bvars_R::build_R;
 
-    void (bvars_R::*prior_3)(const arma::vec&, double, double, const arma::mat&, double, int, bool) = &bvars_R::prior_R;
+    void (bvars_R::*prior_3)(const arma::vec&, double, double, const arma::mat&, const arma::mat&, int, bool) = &bvars_R::prior_R;
 
     SEXP (bvars_R::*forecast_1)(int, bool) = &bvars_R::forecast_R;
     SEXP (bvars_R::*forecast_2)(const arma::mat&, int, bool) = &bvars_R::forecast_R;
@@ -136,11 +136,11 @@ void bvars_R::reset_draws_R()
 }
 
 void bvars_R::prior_R(const arma::vec& coef_prior, double HP_1, double HP_2, 
-                      const arma::mat& Psi_prior, double Xi_psi, int gamma,
+                      const arma::mat& Psi_prior, const arma::mat& Psi_prior_var, int gamma,
                       bool full_cov_prior)
 {
     try {
-        this->prior(coef_prior,HP_1,HP_2,Psi_prior,Xi_psi,gamma,full_cov_prior);
+        this->prior(coef_prior,HP_1,HP_2,Psi_prior,Psi_prior_var,gamma,full_cov_prior);
     } catch( std::exception &ex ) {
         forward_exception_to_r( ex );
     } catch(...) {

diff --git a/src/modules/vars_R.hpp b/src/modules/vars_R.hpp
@@ -49,7 +49,7 @@ class bvars_R : public bm::bvars
         void reset_draws_R();
 
         void prior_R(const arma::vec& coef_prior, double HP_1, double HP_2, 
-                     const arma::mat& Psi_prior, double Xi_Psi, int gamma,
+                     const arma::mat& Psi_prior, const arma::mat& Psi_prior_var, int gamma,
                      bool full_cov_prior);
 
         void gibbs_R(int n_draws, int n_burnin);

diff --git a/tests/var/bvars.R b/tests/var/bvars.R
@@ -17,7 +17,7 @@ psi_prior <- matrix(c(3,6,5),nrow=1)
 
 HP1 <- 0.5
 HP2 <- 0.5
-XiPsi <- 1
+psi_prior_var <- diag(3)
 gamma = 4
 
 bvar_obj <- new(bvars)
@@ -28,7 +28,7 @@ bvar_obj <- new(bvars)
 # p = 1
 
 bvar_obj$build(bvar_data,TRUE,1)
-bvar_obj$prior(coef_prior,HP1,HP2,psi_prior,XiPsi,gamma,FALSE)
+bvar_obj$prior(coef_prior,HP1,HP2,psi_prior,psi_prior_var,gamma,FALSE)
 bvar_obj$gibbs(10000,5000)
 
 IRF(bvar_obj,20,var_names=colnames(USMacroData),save=FALSE)
@@ -40,7 +40,7 @@ FEVD(bvar_obj,20,var_names=colnames(bvar_data),save=FALSE)
 
 bvar_obj$reset_draws()
 bvar_obj$build(bvar_data,TRUE,2)
-bvar_obj$prior(coef_prior,HP1,HP2,psi_prior,XiPsi,gamma,FALSE)
+bvar_obj$prior(coef_prior,HP1,HP2,psi_prior,psi_prior_var,gamma,FALSE)
 bvar_obj$gibbs(10000,5000)
 
 IRF(bvar_obj,20,var_names=colnames(USMacroData),save=FALSE)
@@ -51,7 +51,7 @@ forecast(bvar_obj,shocks=TRUE,var_names=colnames(USMacroData),back_data=10,save=
 
 bvar_obj$reset_draws()
 bvar_obj$build(bvar_data,TRUE,3)
-bvar_obj$prior(coef_prior,HP1,HP2,psi_prior,XiPsi,gamma,FALSE)
+bvar_obj$prior(coef_prior,HP1,HP2,psi_prior,psi_prior_var,gamma,FALSE)
 bvar_obj$gibbs(10000,5000)
 
 IRF(bvar_obj,20,var_names=colnames(USMacroData),save=FALSE)
@@ -62,7 +62,7 @@ forecast(bvar_obj,shocks=TRUE,var_names=colnames(USMacroData),back_data=10,save=
 
 bvar_obj$reset_draws()
 bvar_obj$build(bvar_data,TRUE,4)
-bvar_obj$prior(coef_prior,HP1,HP2,psi_prior,XiPsi,gamma,FALSE)
+bvar_obj$prior(coef_prior,HP1,HP2,psi_prior,psi_prior_var,gamma,FALSE)
 bvar_obj$gibbs(10000,5000)
 
 IRF(bvar_obj,20,var_names=colnames(USMacroData),save=FALSE)