Adapt C++ functions to new SV functionality

NEWS.md

@@ -1,6 +1,6 @@
 # bvartools 0.2.4
 
-* Dropped some error messages to fix a security concern with `Rcpp::stop`.
+* Using an updated version of `Rcpp` to address an issue with `Rcpp::stop`.
 * `stochvol_ocsn2007` can handle multi-column input.
 * `stochvol_ksc1998` can handle multi-column input.
 * Added `post_gamma_state_variance` for posterior simulation of constant error variances of the state equation.

R/RcppExports.R

@@ -646,7 +646,7 @@ stoch_vol <- function(y, h, sigma, h_init, constant) {
 #' Produces a draw of log-volatilities.
 #'
 #' @param y a \eqn{T \times K} matrix containing the time series.
-#' @param h a \eqn{T \times K} vector of log-volatilities.
+#' @param h a \eqn{T \times K} vector of the current draw of log-volatilities.
 #' @param sigma a \eqn{K \times 1} vector of variances of log-volatilities,
 #' where the \eqn{i}th element corresponds to the \eqn{i}th column in \code{y}.
 #' @param h_init a \eqn{K \times 1} vector of the initial states of log-volatilities,
@@ -703,14 +703,18 @@ stochvol_ksc1998 <- function(y, h, sigma, h_init, constant) {
 #' 
 #' Produces a draw of log-volatilities based on Omori, Chib, Shephard and Nakajima (2007).
 #' 
-#' @param y a \eqn{T \times 1} vector containing the time series.
-#' @param h a \eqn{T \times 1} vector of log-volatilities.
-#' @param sigma a numeric of the variance of the log-volatilites.
-#' @param h_init a numeric of the initial state of log-volatilities.
-#' @param constant a numeric of the constant that should be added to \eqn{y^2}
-#' before taking the natural logarithm. See 'Details'.
+#' @param y a \eqn{T \times K} matrix containing the time series.
+#' @param h a \eqn{T \times K} vector of the current draw of log-volatilities.
+#' @param sigma a \eqn{K \times 1} vector of variances of log-volatilities,
+#' where the \eqn{i}th element corresponds to the \eqn{i}th column in \code{y}.
+#' @param h_init a \eqn{K \times 1} vector of the initial states of log-volatilities,
+#' where the \eqn{i}th element corresponds to the \eqn{i}th column in \code{y}.
+#' @param constant a \eqn{K \times 1} vector of constants that should be added to \eqn{y^2}
+#' before taking the natural logarithm. The \eqn{i}th element corresponds to
+#' the \eqn{i}th column in \code{y}. See 'Details'.
 #' 
-#' @details The function produces a posterior draw of the log-volatility \eqn{h} for the model
+#' @details For each column in \code{y} the function produces a posterior
+#' draw of the log-volatility \eqn{h} for the model
 #' \deqn{y_{t} = e^{\frac{1}{2}h_t} \epsilon_{t},}
 #' where \eqn{\epsilon_t \sim N(0, 1)} and \eqn{h_t} is assumed to evolve according to a random walk
 #' \deqn{h_t = h_{t - 1} + u_t,}
@@ -725,6 +729,9 @@ stochvol_ksc1998 <- function(y, h, sigma, h_init, constant) {
 #'   \item Obtain a draw of log-volatilities.
 #' }
 #' 
+#' The implementation is an adaption of the code provided on the website to the textbook
+#' by Chan, Koop, Poirier, and Tobias (2019).
+#' 
 #' @return A vector of log-volatility draws.
 #' 
 #' @examples
@@ -736,7 +743,7 @@ stochvol_ksc1998 <- function(y, h, sigma, h_init, constant) {
 #' h <- matrix(rep(h_init, length(y)))
 #' 
 #' # Obtain draw
-#' stochvol_ocsn2007(y - mean(y), matrix(h), matrix(.05), h_init, matrix(0.0001))
+#' stochvol_ocsn2007(y - mean(y), h, matrix(.05), h_init, matrix(0.0001))
 #' 
 #' @references
 #' 

inst/include/bvartools_RcppExports.h

@@ -46,11 +46,11 @@ namespace bvartools {
         return Rcpp::as<arma::mat >(rcpp_result_gen);
     }
 
-    inline arma::vec stoch_vol(arma::mat y, arma::mat h, arma::vec sigma, arma::vec h_init, arma::vec constant) {
+    inline arma::mat stoch_vol(arma::mat y, arma::mat h, arma::vec sigma, arma::vec h_init, arma::vec constant) {
         typedef SEXP(*Ptr_stoch_vol)(SEXP,SEXP,SEXP,SEXP,SEXP);
         static Ptr_stoch_vol p_stoch_vol = NULL;
         if (p_stoch_vol == NULL) {
-            validateSignature("arma::vec(*stoch_vol)(arma::mat,arma::mat,arma::vec,arma::vec,arma::vec)");
+            validateSignature("arma::mat(*stoch_vol)(arma::mat,arma::mat,arma::vec,arma::vec,arma::vec)");
             p_stoch_vol = (Ptr_stoch_vol)R_GetCCallable("bvartools", "_bvartools_stoch_vol");
         }
         RObject rcpp_result_gen;
@@ -64,7 +64,7 @@ namespace bvartools {
             throw Rcpp::LongjumpException(rcpp_result_gen);
         if (rcpp_result_gen.inherits("try-error"))
             throw Rcpp::exception(Rcpp::as<std::string>(rcpp_result_gen).c_str());
-        return Rcpp::as<arma::vec >(rcpp_result_gen);
+        return Rcpp::as<arma::mat >(rcpp_result_gen);
     }
 
     inline arma::mat stochvol_ksc1998(arma::mat y, arma::mat h, arma::vec sigma, arma::vec h_init, arma::vec constant) {
@@ -88,11 +88,11 @@ namespace bvartools {
         return Rcpp::as<arma::mat >(rcpp_result_gen);
     }
 
-    inline arma::mat stochvol_ocsn2007(arma::vec y, arma::vec h, double sigma, double h_init, double constant) {
+    inline arma::mat stochvol_ocsn2007(arma::mat y, arma::mat h, arma::vec sigma, arma::vec h_init, arma::vec constant) {
         typedef SEXP(*Ptr_stochvol_ocsn2007)(SEXP,SEXP,SEXP,SEXP,SEXP);
         static Ptr_stochvol_ocsn2007 p_stochvol_ocsn2007 = NULL;
         if (p_stochvol_ocsn2007 == NULL) {
-            validateSignature("arma::mat(*stochvol_ocsn2007)(arma::vec,arma::vec,double,double,double)");
+            validateSignature("arma::mat(*stochvol_ocsn2007)(arma::mat,arma::mat,arma::vec,arma::vec,arma::vec)");
             p_stochvol_ocsn2007 = (Ptr_stochvol_ocsn2007)R_GetCCallable("bvartools", "_bvartools_stochvol_ocsn2007");
         }
         RObject rcpp_result_gen;

src/RcppExports.cpp

@@ -273,7 +273,7 @@ BEGIN_RCPP
 END_RCPP
 }
 // stoch_vol
-arma::vec stoch_vol(arma::mat y, arma::mat h, arma::vec sigma, arma::vec h_init, arma::vec constant);
+arma::mat stoch_vol(arma::mat y, arma::mat h, arma::vec sigma, arma::vec h_init, arma::vec constant);
 static SEXP _bvartools_stoch_vol_try(SEXP ySEXP, SEXP hSEXP, SEXP sigmaSEXP, SEXP h_initSEXP, SEXP constantSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
@@ -349,15 +349,15 @@ RcppExport SEXP _bvartools_stochvol_ksc1998(SEXP ySEXP, SEXP hSEXP, SEXP sigmaSE
     return rcpp_result_gen;
 }
 // stochvol_ocsn2007
-arma::mat stochvol_ocsn2007(arma::vec y, arma::vec h, double sigma, double h_init, double constant);
+arma::mat stochvol_ocsn2007(arma::mat y, arma::mat h, arma::vec sigma, arma::vec h_init, arma::vec constant);
 static SEXP _bvartools_stochvol_ocsn2007_try(SEXP ySEXP, SEXP hSEXP, SEXP sigmaSEXP, SEXP h_initSEXP, SEXP constantSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
-    Rcpp::traits::input_parameter< arma::vec >::type y(ySEXP);
-    Rcpp::traits::input_parameter< arma::vec >::type h(hSEXP);
-    Rcpp::traits::input_parameter< double >::type sigma(sigmaSEXP);
-    Rcpp::traits::input_parameter< double >::type h_init(h_initSEXP);
-    Rcpp::traits::input_parameter< double >::type constant(constantSEXP);
+    Rcpp::traits::input_parameter< arma::mat >::type y(ySEXP);
+    Rcpp::traits::input_parameter< arma::mat >::type h(hSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type sigma(sigmaSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type h_init(h_initSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type constant(constantSEXP);
     rcpp_result_gen = Rcpp::wrap(stochvol_ocsn2007(y, h, sigma, h_init, constant));
     return rcpp_result_gen;
 END_RCPP_RETURN_ERROR
@@ -421,9 +421,9 @@ static int _bvartools_RcppExport_validate(const char* sig) {
     static std::set<std::string> signatures;
     if (signatures.empty()) {
         signatures.insert("arma::mat(*kalman_dk)(arma::mat,arma::mat,arma::mat,arma::mat,arma::mat,arma::vec,arma::mat)");
-        signatures.insert("arma::vec(*stoch_vol)(arma::mat,arma::mat,arma::vec,arma::vec,arma::vec)");
+        signatures.insert("arma::mat(*stoch_vol)(arma::mat,arma::mat,arma::vec,arma::vec,arma::vec)");
         signatures.insert("arma::mat(*stochvol_ksc1998)(arma::mat,arma::mat,arma::vec,arma::vec,arma::vec)");
-        signatures.insert("arma::mat(*stochvol_ocsn2007)(arma::vec,arma::vec,double,double,double)");
+        signatures.insert("arma::mat(*stochvol_ocsn2007)(arma::mat,arma::mat,arma::vec,arma::vec,arma::vec)");
     }
     return signatures.find(sig) != signatures.end();
 }

src/bvaralg.cpp

@@ -287,15 +287,16 @@ Rcpp::List bvaralg(Rcpp::List object) {
   if (covar && psi_varsel) {
     draws_lambda_a0 = arma::zeros<arma::mat>(n_psi, iter);  
   }
-
+  
   // Start Gibbs sampler
   for (int draw = 0; draw < draws; draw++) {
 
     if (draw % 20 == 0) { // Check for user interruption every now and then
       Rcpp::checkUserInterrupt();
     }
 
-    // Draw coefficients ----
+    ////////////////////////////////////////////////////////////////////////////
+    // Draw coefficients
     if (use_a) {
 
       if (bvs) {
@@ -366,7 +367,9 @@ Rcpp::List bvaralg(Rcpp::List object) {
       u = y;
     }
 
-    // Covariances
+    ////////////////////////////////////////////////////////////////////////////
+    // Draw error covariances
+
     if (covar) {
 
       // Prepare data
@@ -451,17 +454,12 @@ Rcpp::List bvaralg(Rcpp::List object) {
       u = Psi * u;
     }
 
+    ////////////////////////////////////////////////////////////////////////////
+    // Draw error variances
+
     if (sv) {
 
-      // Draw variances  
-      h = bvartools::stoch_vol(u, h, sigma_h, h_init, h_constant);
-      diag_omega_i.diag() = 1 / exp(arma::vectorise(h.t()));
-      if (covar) {
-        diag_Psi = arma::kron(diag_tt, Psi);
-        diag_sigma_i = arma::trans(diag_Psi) * diag_omega_i * diag_Psi;
-      } else {
-        diag_sigma_i = diag_omega_i;
-      }
+      h = bvartools::stochvol_ksc1998(arma::trans(u), h, sigma_h, h_init, h_constant);
 
       // Draw sigma_h
       h_lag.row(0) = h_init.t();
@@ -481,24 +479,39 @@ Rcpp::List bvaralg(Rcpp::List object) {
     } else {
 
       if (use_gamma) {
-
         sse = u * u.t();
         for (int i = 0; i < k; i++) {
           omega_i(i, i) = arma::randg<double>(arma::distr_param(sigma_post_shape(i), 1 / arma::as_scalar(sigma_prior_rate(i) + sse(i, i) * 0.5)));
         }
+      }
+
+    }
+
+    ////////////////////////////////////////////////////////////////////////////
+    // Combine Psi and Omega resp. draw from Wishart
+
+    if (sv) {
+      diag_omega_i.diag() = 1 / exp(arma::vectorise(arma::trans(h)));
+      if (covar) {
+        diag_Psi = arma::kron(diag_tt, Psi);
+        diag_sigma_i = arma::trans(diag_Psi) * diag_omega_i * diag_Psi;
+      } else {
+        diag_sigma_i = diag_omega_i;
+      }
+    } else {
+      if (use_gamma) {
         if (covar) {
-          diag_omega_i = arma::kron(diag_tt, omega_i);
-          sigma_i = arma::trans(Psi) * omega_i * Psi;
+          diag_omega_i = arma::kron(diag_tt, omega_i); // Used if covar is estimated
+          sigma_i = arma::trans(Psi) * omega_i * Psi; // Update sigma
         } else {
-          sigma_i = omega_i;
+          sigma_i = omega_i; // Since no covar, sigma = omega
         }
-
       } else {
+        // Draw from Wishart
         sigma_i = arma::wishrnd(arma::solve(sigma_prior_scale + u * u.t(), diag_k), sigma_post_df);
       }
-
+      // Final sigma block diagonal matrix for block A
       diag_sigma_i = arma::kron(diag_tt, sigma_i);
-
     }
 
     // Store draws
@@ -613,4 +626,19 @@ Rcpp::List bvaralg(Rcpp::List object) {
 
   return result;
   // return Rcpp::List::create(Rcpp::Named("test") = sigma_i);
-}
+}
+
+/*** R
+
+data("us_macrodata")
+
+object <- gen_var(us_macrodata, p = 0, deterministic = "none",
+                  sv = TRUE,
+                  iterations = 20, burnin = 10)
+
+object <- add_priors(object,
+                     sigma = list(shape = 3, rate = .4, mu = 10, v_i = .01, sigma_h = .05, constant = .0001))
+
+.bvaralg(object)
+
+*/