From 2dc13d681b17402e9427f3fa2c9d4ad25d87c52f Mon Sep 17 00:00:00 2001 From: Franz Mohr Date: Sat, 6 Jan 2024 10:46:33 +0100 Subject: [PATCH] Adapt C++ functions to new SV functionality --- NEWS.md | 2 +- R/RcppExports.R | 25 ++++++--- inst/include/bvartools_RcppExports.h | 10 ++-- man/stochvol_ksc1998.Rd | 2 +- man/stochvol_ocsn2007.Rd | 23 +++++--- src/RcppExports.cpp | 18 +++--- src/bvaralg.cpp | 68 ++++++++++++++++------- src/bvartvpalg.cpp | 44 +++++++++------ src/bvecalg.cpp | 42 ++++++++------ src/bvectvpalg.cpp | 31 +++++++---- src/stoch_vol.cpp | 2 +- src/stochvol_ksc1998.cpp | 14 +++-- src/stochvol_ocsn2007.cpp | 83 +++++++++++++++++----------- 13 files changed, 230 insertions(+), 134 deletions(-) diff --git a/NEWS.md b/NEWS.md index 3ee3c1e..d29f825 100644 --- a/NEWS.md +++ b/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. diff --git a/R/RcppExports.R b/R/RcppExports.R index 7bfcf5b..1823716 100644 --- a/R/RcppExports.R +++ b/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 #' diff --git a/inst/include/bvartools_RcppExports.h b/inst/include/bvartools_RcppExports.h index 21cbb67..0022c1f 100644 --- a/inst/include/bvartools_RcppExports.h +++ b/inst/include/bvartools_RcppExports.h @@ -46,11 +46,11 @@ namespace bvartools { return Rcpp::as(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(rcpp_result_gen).c_str()); - return Rcpp::as(rcpp_result_gen); + return Rcpp::as(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(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; diff --git a/man/stochvol_ksc1998.Rd b/man/stochvol_ksc1998.Rd index 725ef2a..cd463f9 100644 --- a/man/stochvol_ksc1998.Rd +++ b/man/stochvol_ksc1998.Rd @@ -9,7 +9,7 @@ stochvol_ksc1998(y, h, sigma, h_init, constant) \arguments{ \item{y}{a \eqn{T \times K} matrix containing the time series.} -\item{h}{a \eqn{T \times K} vector of log-volatilities.} +\item{h}{a \eqn{T \times K} vector of the current draw of log-volatilities.} \item{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}.} diff --git a/man/stochvol_ocsn2007.Rd b/man/stochvol_ocsn2007.Rd index 3a75993..5bf984a 100644 --- a/man/stochvol_ocsn2007.Rd +++ b/man/stochvol_ocsn2007.Rd @@ -7,16 +7,19 @@ stochvol_ocsn2007(y, h, sigma, h_init, constant) } \arguments{ -\item{y}{a \eqn{T \times 1} vector containing the time series.} +\item{y}{a \eqn{T \times K} matrix containing the time series.} -\item{h}{a \eqn{T \times 1} vector of log-volatilities.} +\item{h}{a \eqn{T \times K} vector of the current draw of log-volatilities.} -\item{sigma}{a numeric of the variance of the log-volatilites.} +\item{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}.} -\item{h_init}{a numeric of the initial state of log-volatilities.} +\item{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}.} -\item{constant}{a numeric of the constant that should be added to \eqn{y^2} -before taking the natural logarithm. See 'Details'.} +\item{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'.} } \value{ A vector of log-volatility draws. @@ -25,7 +28,8 @@ A vector of log-volatility draws. Produces a draw of log-volatilities based on Omori, Chib, Shephard and Nakajima (2007). } \details{ -The function produces a posterior draw of the log-volatility \eqn{h} for the model +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,} @@ -39,6 +43,9 @@ following steps: approximating the log-\eqn{\chi_1^2} distribution. \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). } \examples{ data("us_macrodata") @@ -49,7 +56,7 @@ h_init <- matrix(log(var(y))) 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{ diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index c69c0ae..6ab702c 100644 --- a/src/RcppExports.cpp +++ b/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 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(); } diff --git a/src/bvaralg.cpp b/src/bvaralg.cpp index 01df564..4224a54 100644 --- a/src/bvaralg.cpp +++ b/src/bvaralg.cpp @@ -287,7 +287,7 @@ Rcpp::List bvaralg(Rcpp::List object) { if (covar && psi_varsel) { draws_lambda_a0 = arma::zeros(n_psi, iter); } - + // Start Gibbs sampler for (int draw = 0; draw < draws; draw++) { @@ -295,7 +295,8 @@ Rcpp::List bvaralg(Rcpp::List object) { 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(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); -} \ No newline at end of file +} + +/*** 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) + +*/ diff --git a/src/bvartvpalg.cpp b/src/bvartvpalg.cpp index 8841ffe..30e3ef5 100644 --- a/src/bvartvpalg.cpp +++ b/src/bvartvpalg.cpp @@ -406,17 +406,13 @@ Rcpp::List bvartvpalg(Rcpp::List object) { u = arma::reshape(Psi * u_vec, k, tt); } + //////////////////////////////////////////////////////////////////////////// + // Draw measurement 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_sigma_u_i = arma::trans(Psi) * diag_omega_i * Psi; - } else { - diag_sigma_u_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(); h_lag.rows(1, tt - 1) = h.rows(0, tt - 2); @@ -425,35 +421,51 @@ Rcpp::List bvartvpalg(Rcpp::List object) { for (int i = 0; i < k; i++) { sigma_h(i) = 1 / arma::randg(arma::distr_param(sigma_post_shape(i), sigma_post_scale(i))); } - + // Draw h_init sigma_h_i = arma::diagmat(1 / sigma_h); h_init_post_v = sigma_prior_vi + sigma_h_i; h_init_post_mu = arma::solve(h_init_post_v, sigma_prior_vi * sigma_prior_mu + sigma_h_i * h.row(0).t()); h_init = h_init_post_mu + arma::solve(arma::chol(h_init_post_v), arma::randn(k)); - + } else { - + if (use_gamma) { - sse = u * u.t(); for (int i = 0; i < k; i++) { omega_i(i, i) = arma::randg(arma::distr_param(sigma_post_shape(i), 1 / arma::as_scalar(sigma_prior_rate(i) + sse(i, i) * 0.5))); } - diag_omega_i = arma::kron(diag_tt, omega_i); + } + + } + + //////////////////////////////////////////////////////////////////////////// + // Combine Psi and Omega resp. draw from Wishart + + if (sv) { + diag_omega_i.diag() = 1 / exp(arma::vectorise(h.t())); + if (covar) { + diag_sigma_u_i = arma::trans(Psi) * diag_omega_i * Psi; + } else { + diag_sigma_u_i = diag_omega_i; + } + } else { + if (use_gamma) { if (covar) { diag_sigma_u_i = arma::trans(Psi) * diag_omega_i * Psi; } else { sigma_u_i = omega_i; diag_sigma_u_i = diag_omega_i; } - } else { sigma_u_i = arma::wishrnd(arma::solve(sigma_prior_scale + u * u.t(), diag_k), sigma_post_df); diag_sigma_u_i = arma::kron(diag_tt, sigma_u_i); } } + //////////////////////////////////////////////////////////////////////////// + // Draw state variances and initial states + if (use_a) { // Draw sigma_v_i a_lag.subvec(0, n_tot - 1) = a_init; diff --git a/src/bvecalg.cpp b/src/bvecalg.cpp index 3339552..035b228 100644 --- a/src/bvecalg.cpp +++ b/src/bvecalg.cpp @@ -589,18 +589,13 @@ Rcpp::List bvecalg(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(); h_lag.rows(1, tt - 1) = h.rows(0, tt - 2); @@ -609,28 +604,43 @@ Rcpp::List bvecalg(Rcpp::List object) { for (int i = 0; i < k; i++) { sigma_h(i) = 1 / arma::randg(arma::distr_param(sigma_post_shape(i), sigma_post_scale(i))); } - + // Draw h_init sigma_h_i = arma::diagmat(1 / sigma_h); h_init_post_v = sigma_prior_vi + sigma_h_i; h_init_post_mu = arma::solve(h_init_post_v, sigma_prior_vi * sigma_prior_mu + sigma_h_i * h.row(0).t()); h_init = h_init_post_mu + arma::solve(arma::chol(h_init_post_v), arma::randn(k)); - + } else { - + if (use_gamma) { - sse = u * u.t(); for (int i = 0; i < k; i++) { omega_i(i, i) = arma::randg(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(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; + } + } else { + if (use_gamma) { if (covar) { diag_omega_i = arma::kron(diag_tt, omega_i); sigma_i = arma::trans(Psi) * omega_i * Psi; } else { sigma_i = omega_i; } - } else { if (use_rr) { sigma_i = arma::wishrnd(arma::solve(coint_v_i * alpha * (beta.t() * p_tau_i * beta) * alpha.t() + u * u.t(), diag_k), sigma_post_df); @@ -638,9 +648,7 @@ Rcpp::List bvecalg(Rcpp::List object) { sigma_i = arma::wishrnd(arma::solve(sigma_prior_scale + u * u.t(), diag_k), sigma_post_df); } } - diag_sigma_i = arma::kron(diag_tt, sigma_i); - } // Update g_i diff --git a/src/bvectvpalg.cpp b/src/bvectvpalg.cpp index eeb9b4e..97fb170 100644 --- a/src/bvectvpalg.cpp +++ b/src/bvectvpalg.cpp @@ -625,16 +625,12 @@ Rcpp::List bvectvpalg(Rcpp::List object) { u = arma::reshape(Psi * arma::vectorise(u), k, tt); } + //////////////////////////////////////////////////////////////////////////// + // Draw measurement 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_sigma_u_i = arma::trans(Psi) * diag_omega_i * Psi; - } else { - diag_sigma_u_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(); @@ -654,19 +650,32 @@ Rcpp::List bvectvpalg(Rcpp::List object) { } else { if (use_gamma) { - sse = u * u.t(); for (int i = 0; i < k; i++) { omega_i(i, i) = arma::randg(arma::distr_param(sigma_post_shape(i), 1 / arma::as_scalar(sigma_prior_rate(i) + sse(i, i) * 0.5))); } - diag_omega_i = arma::kron(diag_tt, omega_i); + } + + } + + //////////////////////////////////////////////////////////////////////////// + // Combine Psi and Omega resp. draw from Wishart + + if (sv) { + diag_omega_i.diag() = 1 / exp(arma::vectorise(h.t())); + if (covar) { + diag_sigma_u_i = arma::trans(Psi) * diag_omega_i * Psi; + } else { + diag_sigma_u_i = diag_omega_i; + } + } else { + if (use_gamma) { if (covar) { diag_sigma_u_i = arma::trans(Psi) * diag_omega_i * Psi; } else { sigma_u_i = omega_i; diag_sigma_u_i = diag_omega_i; } - } else { sigma_u_i = arma::wishrnd(arma::solve(sigma_prior_scale + u * u.t(), diag_k), sigma_post_df); diag_sigma_u_i = arma::kron(diag_tt, sigma_u_i); diff --git a/src/stoch_vol.cpp b/src/stoch_vol.cpp index 9f2203f..f4eacc5 100644 --- a/src/stoch_vol.cpp +++ b/src/stoch_vol.cpp @@ -38,7 +38,7 @@ //' with ARCH models. \emph{Review of Economic Studies 65}(3), 361--393. \doi{10.1111/1467-937X.00050} //' // [[Rcpp::export]] -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) { return bvartools::stochvol_ksc1998(y, h, sigma, h_init, constant); } diff --git a/src/stochvol_ksc1998.cpp b/src/stochvol_ksc1998.cpp index 232d00b..1231081 100644 --- a/src/stochvol_ksc1998.cpp +++ b/src/stochvol_ksc1998.cpp @@ -7,7 +7,7 @@ //' 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, @@ -59,9 +59,16 @@ // [[Rcpp::export]] arma::mat stochvol_ksc1998(arma::mat y, arma::mat h, arma::vec sigma, arma::vec h_init, arma::vec constant) { + // Checks if (y.has_nan()) { Rcpp::stop("Argument 'y' contains NAs."); } + if (y.n_rows != h.n_rows) { + Rcpp::stop("Arguments 'y' and 'h' do not have the same number of rows."); + } + if (y.n_cols != h.n_cols) { + Rcpp::stop("Arguments 'y' and 'h' do not have the same number of columns."); + } // Components of the mixture model arma::rowvec p_i(7), mu(7), sigma2(7); @@ -79,10 +86,6 @@ arma::mat stochvol_ksc1998(arma::mat y, arma::mat h, arma::vec sigma, arma::vec hh.diag(-1) = -arma::ones(tt - 1); hh = hh.t() * hh; - if (y.n_rows != h.n_rows) { - Rcpp::stop("Arguments 'y' and 'h' do not have the same length."); - } - for (int i = 0; i < k; i++) { // Prepare series y.col(i) = log(arma::pow(y.col(i), 2) + constant(i)); @@ -113,4 +116,5 @@ h <- t(matrix(h_init, 3, nrow(y))) sigma_h <- rep(.05, 3) const <- rep(.0001, 3) stochvol_ksc1998(y, h, sigma_h, h_init, const) + ***/ diff --git a/src/stochvol_ocsn2007.cpp b/src/stochvol_ocsn2007.cpp index b7f801e..05968ad 100644 --- a/src/stochvol_ocsn2007.cpp +++ b/src/stochvol_ocsn2007.cpp @@ -6,14 +6,18 @@ //' //' 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,} @@ -28,6 +32,9 @@ //' \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 @@ -39,7 +46,7 @@ //' 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 //' @@ -50,17 +57,17 @@ //' \emph{Journal of Econometrics 140}(2), 425--449. \doi{10.1016/j.jeconom.2006.07.008} //' // [[Rcpp::export]] -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) { + // Checks if (y.has_nan()) { Rcpp::stop("Argument 'y' contains NAs."); } - - // Prepare series - y = log(arma::pow(y, 2) + constant); - arma::uword tt = y.n_elem; - if (y.n_elem != h.n_elem) { - Rcpp::stop("Arguments 'y' and 'h' do not have the same length."); + if (y.n_rows != h.n_rows) { + Rcpp::stop("Arguments 'y' and 'h' do not have the same number of rows."); + } + if (y.n_cols != h.n_cols) { + Rcpp::stop("Arguments 'y' and 'h' do not have the same number of columns."); } // Components of the mixture model @@ -76,30 +83,44 @@ arma::mat stochvol_ocsn2007(arma::vec y, arma::vec h, double sigma, double h_ini p_i(8) = 0.01575; mu(8) = -8.68384; sigma2(8) = 4.16591; p_i(9) = 0.00115; mu(9) = -14.65000; sigma2(9) = 7.33342; - // Choose, which component should be used per observation - arma::mat q = arma::repmat(p_i, tt, 1) % arma::normpdf(arma::repmat(y, 1, 10), arma::repmat(h, 1, 10) + arma::repmat(mu, tt, 1), arma::repmat(sqrt(sigma2), tt, 1)); - q = q / arma::repmat(sum(q, 1), 1, 10); - arma::umat s = 10 - sum(arma::repmat(arma::randu(tt), 1, 10) < cumsum(q, 1), 1); - // umat is a matrix index + int k = y.n_cols; + int tt = y.n_rows; + arma::mat q, sigh_hh, sigs, post_h_v, post_h_mu; + arma::umat s; - // Sample log-volatility arma::mat hh = arma::eye(tt, tt); hh.diag(-1) = -arma::ones(tt - 1); - arma::mat sigh_hh = hh.t() * hh / sigma; - arma::mat sigs = arma::eye(tt, tt); - sigs.diag() = 1 / sigma2.elem(s); - arma::mat post_h_v = sigh_hh + sigs; - arma::mat post_h_mu = arma::solve(post_h_v, sigh_hh * arma::ones(tt) * h_init + sigs * (y - mu.elem(s))); + hh = hh.t() * hh; + + for (int i = 0; i < k; i++) { + // Prepare series + y.col(i) = log(arma::pow(y.col(i), 2) + constant(i)); + + // Sample s + q = arma::repmat(p_i, tt, 1) % arma::normpdf(arma::repmat(y.col(i), 1, 10), arma::repmat(h.col(i), 1, 10) + arma::repmat(mu, tt, 1), arma::repmat(sqrt(sigma2), tt, 1)); + q = q / arma::repmat(sum(q, 1), 1, 10); + s = 10 - sum(arma::repmat(arma::randu(tt), 1, 10) < cumsum(q, 1), 1); + + // Sample log-volatility + sigh_hh = hh / sigma(i); + sigs = arma::eye(tt, tt); + sigs.diag() = 1 / sigma2.elem(s); + post_h_v = sigh_hh + sigs; + post_h_mu = arma::solve(post_h_v, sigh_hh * arma::ones(tt) * h_init(i) + sigs * (y.col(i) - mu.elem(s))); + h.col(i) = post_h_mu + arma::solve(arma::chol(post_h_v), arma::randn(tt)); + } - return post_h_mu + arma::solve(arma::chol(arma::mat(post_h_v)), arma::randn(tt)); + return h; } /*** R data("us_macrodata") -aud <- us_macrodata[, 1] -h_init <- matrix(log(var(aud))) -h <- matrix(rep(h_init, length(aud))) -stochvol_ocsn2007(aud - mean(aud), h, matrix(.05), h_init, matrix(0.0001)) +y <- us_macrodata +h_init <- log(diag(var(y))) +h <- t(matrix(h_init, 3, nrow(y))) +sigma_h <- rep(.05, 3) +const <- rep(.0001, 3) +stochvol_ocsn2007(y, h, sigma_h, h_init, const) ***/