From e8194e04d950d76a75b95166bb8c4d998c21bb44 Mon Sep 17 00:00:00 2001 From: Franz Mohr Date: Mon, 11 Dec 2023 20:23:42 +0100 Subject: [PATCH] Multi column input for stochvol functions --- NEWS.md | 2 + src/stochvol_ksc1998.cpp | 87 +++++++++++++++++++++++----------------- 2 files changed, 53 insertions(+), 36 deletions(-) diff --git a/NEWS.md b/NEWS.md index f1031e2..d1ac63d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,7 @@ # bvartools 0.2.4 +* `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. * Added `post_normal_covar_tvp` for posterior simulation of time varying, lower triangular covariance matrices. * Added `post_normal_covar_const` for posterior simulation of constant, lower triangular covariance matrices. diff --git a/src/stochvol_ksc1998.cpp b/src/stochvol_ksc1998.cpp index 8b031f8..94db48d 100644 --- a/src/stochvol_ksc1998.cpp +++ b/src/stochvol_ksc1998.cpp @@ -3,23 +3,27 @@ // [[Rcpp::interfaces(r, cpp)]] //' Stochastic Volatility -//' +//' //' 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 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'. //' -//' @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'. -//' -//' @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,} //' with \eqn{u_t \sim N(0, \sigma^2)}. //' -//' The implementation follows the algorithm of Kim, Shephard and Chip (1998) and performs the +//' The implementation is based on the algorithm of Kim, Shephard and Chip (1998) and performs the //' following steps: //' \enumerate{ //' \item Perform the transformation \eqn{y_t^* = ln(y_t^2 + constant)}. @@ -28,7 +32,8 @@ //' \item Obtain a draw of log-volatilities. //' } //' -//' The implementation follows the code provided on the website to the textbook by Chan, Koop, Poirier, and Tobias (2019). +//' The implementation follows the code provided on the website to the textbook +//' by Chan, Koop, Poirier, and Tobias (2019). //' //' @return A vector of log-volatility draws. //' @@ -52,19 +57,12 @@ //' with ARCH models. \emph{Review of Economic Studies 65}(3), 361--393. \doi{10.1111/1467-937X.00050} //' // [[Rcpp::export]] -arma::vec stochvol_ksc1998(arma::vec y, arma::vec h, double sigma, double h_init, double constant) { +arma::mat stochvol_ksc1998(arma::mat y, arma::mat h, arma::vec sigma, arma::vec h_init, arma::vec constant) { 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 (tt != h.n_elem) { - Rcpp::stop("Arguments 'y' and 'h' do not have the same length."); - } - // Components of the mixture model arma::rowvec p_i(7), mu(7), sigma2(7); p_i(0) = 0.00730; p_i(1) = 0.10556; p_i(2) = 0.00002; p_i(3) = 0.04395; p_i(4) = 0.34001; p_i(5) = 0.24566; p_i(6) = 0.25750; @@ -72,30 +70,47 @@ arma::vec stochvol_ksc1998(arma::vec y, arma::vec h, double sigma, double h_init mu = mu - 1.2704; // Already subtract the constant sigma2(0) = 5.79596; sigma2(1) = 2.61369; sigma2(2) = 5.17950; sigma2(3) = 0.16735; sigma2(4) = 0.64009; sigma2(5) = 0.34023; sigma2(6) = 1.26261; - // Sample s - arma::mat q = arma::repmat(p_i, tt, 1) % arma::normpdf(arma::repmat(y, 1, 7), arma::repmat(h, 1, 7) + arma::repmat(mu, tt, 1), arma::repmat(sqrt(sigma2), tt, 1)); - q = q / arma::repmat(sum(q, 1), 1, 7); - arma::umat s = 7 - sum(arma::repmat(arma::randu(tt), 1, 7) < cumsum(q, 1), 1); + 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); hh = hh.t() * hh; - arma::mat sigh_hh = 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))); - return post_h_mu + arma::solve(arma::chol(arma::mat(post_h_v)), arma::randn(tt)); + if (tt != 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)); + + // Sample s + q = arma::repmat(p_i, tt, 1) % arma::normpdf(arma::repmat(y.col(i), 1, 7), arma::repmat(h.col(i), 1, 7) + arma::repmat(mu, tt, 1), arma::repmat(sqrt(sigma2), tt, 1)); + q = q / arma::repmat(sum(q, 1), 1, 7); + s = 7 - sum(arma::repmat(arma::randu(tt), 1, 7) < 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 h; } /*** R data("us_macrodata") -y <- us_macrodata[, 1] -h_init <- log(var(y)) -h <- rep(h_init, length(y)) -stochvol_ksc1998(y - mean(y), h, .05, h_init, 0.0001) - -***/ \ No newline at end of file +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_ksc1998(y, h, sigma_h, h_init, const) +***/