From 4144d63d35b7e2560771f28bb8d372c4d0914a4e Mon Sep 17 00:00:00 2001 From: Franz Mohr Date: Thu, 19 Oct 2023 21:11:19 +0200 Subject: [PATCH] Add post_gamma_state_variance --- R/RcppExports.R | 32 ++++++++++------ man/post_gamma_state_variance.Rd | 31 +++++++++------ src/RcppExports.cpp | 13 ++++--- src/post_gamma_state_variance.cpp | 63 ++++++++++++++++++------------- 4 files changed, 84 insertions(+), 55 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 690e24d..5fa2cc0 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -430,15 +430,23 @@ post_coint_kls_sur <- function(y, beta, w, sigma_i, v_i, p_tau_i, g_i, x = NULL, #' Posterior Draws of Error Variances #' -#' Produces a draw of error variances from a gamma posterior density. +#' Produces a draw of the constant diagonal error variance matrix of the +#' state equation of a state space model using an inverse gamma posterior density. #' -#' @param phi a \eqn{K \times T} matrix of time varying parameter draws. -#' @param phi_init a \eqn{K \times 1} vector of initial states. +#' @param a a \eqn{K \times T} matrix of time varying parameter draws. +#' @param a_init a \eqn{K \times 1} vector of initial states. #' @param shape_prior a \eqn{K \times 1} vector of prior shape parameters. #' @param rate_prior a \eqn{K \times 1} vector of prior rate parameters. +#' @param inverse logical. If \code{TRUE}, the function returns the precision matrix, +#' i.e. the inverse of the variance matrix. Defaults to \code{FALSE}. #' -#' @details The function produces a posterior draw of the variaces vector \eqn{a} for the model -#' Follow description in Chan eta al. +#' @details For the state space model with state equation +#' \deqn{a_t = a_{t-1} + v} +#' and measurement equation +#' \deqn{y_t = Z_{t} a_t + u_t} +#' with \eqn{v_t \sim N(0, \Sigma_{v})} and \eqn{u_t \sim N(0, \Sigma_{u,t})} +#' the function produces a draw of the constant diagonal error variances matrix of the +#' state equation \eqn{\Simga_v}. #' #' @references #' Chan, J., Koop, G., Poirier, D. J., & Tobias J. L. (2019). \emph{Bayesian econometric methods} @@ -451,25 +459,25 @@ post_coint_kls_sur <- function(y, beta, w, sigma_i, v_i, p_tau_i, g_i, x = NULL, #' set.seed(1234) # Set RNG seed #' #' # Generate artificial data according to a random walk -#' phi <- matrix(rnorm(k), k, tt + 1) +#' a <- matrix(rnorm(k), k, tt + 1) #' for (i in 2:(tt + 1)) { -#' phi[, i] <- phi[, i - 1] + rnorm(k, 0, sqrt(1 / 100)) +#' a[, i] <- a[, i - 1] + rnorm(k, 0, sqrt(1 / 100)) #' } #' -#' phi_init <- matrix(phi[, 1]) # Define inital state -#' phi <- phi[, -1] # Drop initial state from main sample +#' a_init <- matrix(a[, 1]) # Define inital state +#' a <- a[, -1] # Drop initial state from main sample #' #' # Define priors #' shape_prior <- matrix(1, k) #' rate_prior <- matrix(.0001, k) #' #' # Obtain posterior draw -#' post_gamma_state_variance(phi, phi_init, shape_prior, rate_prior) +#' post_gamma_state_variance(a, a_init, shape_prior, rate_prior) #' #' @return A matrix. #' -post_gamma_state_variance <- function(phi, phi_init, shape_prior, rate_prior) { - .Call(`_bvartools_post_gamma_state_variance`, phi, phi_init, shape_prior, rate_prior) +post_gamma_state_variance <- function(a, a_init, shape_prior, rate_prior, inverse = FALSE) { + .Call(`_bvartools_post_gamma_state_variance`, a, a_init, shape_prior, rate_prior, inverse) } #' Posterior Draw from a Normal Distribution diff --git a/man/post_gamma_state_variance.Rd b/man/post_gamma_state_variance.Rd index 4bd385d..aa95682 100644 --- a/man/post_gamma_state_variance.Rd +++ b/man/post_gamma_state_variance.Rd @@ -4,26 +4,35 @@ \alias{post_gamma_state_variance} \title{Posterior Draws of Error Variances} \usage{ -post_gamma_state_variance(phi, phi_init, shape_prior, rate_prior) +post_gamma_state_variance(a, a_init, shape_prior, rate_prior, inverse = FALSE) } \arguments{ -\item{phi}{a \eqn{K \times T} matrix of time varying parameter draws.} +\item{a}{a \eqn{K \times T} matrix of time varying parameter draws.} -\item{phi_init}{a \eqn{K \times 1} vector of initial states.} +\item{a_init}{a \eqn{K \times 1} vector of initial states.} \item{shape_prior}{a \eqn{K \times 1} vector of prior shape parameters.} \item{rate_prior}{a \eqn{K \times 1} vector of prior rate parameters.} + +\item{inverse}{logical. If \code{TRUE}, the function returns the precision matrix, +i.e. the inverse of the variance matrix. Defaults to \code{FALSE}.} } \value{ A matrix. } \description{ -Produces a draw of error variances from a gamma posterior density. +Produces a draw of the constant diagonal error variance matrix of the +state equation of a state space model using an inverse gamma posterior density. } \details{ -The function produces a posterior draw of the variaces vector \eqn{a} for the model -Follow description in Chan eta al. +For the state space model with state equation +\deqn{a_t = a_{t-1} + v} +and measurement equation +\deqn{y_t = Z_{t} a_t + u_t} +with \eqn{v_t \sim N(0, \Sigma_{v})} and \eqn{u_t \sim N(0, \Sigma_{u,t})} +the function produces a draw of the constant diagonal error variances matrix of the +state equation \eqn{\Simga_v}. } \examples{ k <- 10 # Number of artificial coefficients @@ -32,20 +41,20 @@ tt <- 1000 # Number of observations set.seed(1234) # Set RNG seed # Generate artificial data according to a random walk -phi <- matrix(rnorm(k), k, tt + 1) +a <- matrix(rnorm(k), k, tt + 1) for (i in 2:(tt + 1)) { - phi[, i] <- phi[, i - 1] + rnorm(k, 0, sqrt(1 / 100)) + a[, i] <- a[, i - 1] + rnorm(k, 0, sqrt(1 / 100)) } -phi_init <- matrix(phi[, 1]) # Define inital state -phi <- phi[, -1] # Drop initial state from main sample +a_init <- matrix(a[, 1]) # Define inital state +a <- a[, -1] # Drop initial state from main sample # Define priors shape_prior <- matrix(1, k) rate_prior <- matrix(.0001, k) # Obtain posterior draw -post_gamma_state_variance(phi, phi_init, shape_prior, rate_prior) +post_gamma_state_variance(a, a_init, shape_prior, rate_prior) } \references{ diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index e147392..65a1cdf 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -213,16 +213,17 @@ BEGIN_RCPP END_RCPP } // post_gamma_state_variance -arma::mat post_gamma_state_variance(arma::mat phi, arma::vec phi_init, arma::vec shape_prior, arma::vec rate_prior); -RcppExport SEXP _bvartools_post_gamma_state_variance(SEXP phiSEXP, SEXP phi_initSEXP, SEXP shape_priorSEXP, SEXP rate_priorSEXP) { +arma::mat post_gamma_state_variance(arma::mat a, arma::vec a_init, arma::vec shape_prior, arma::vec rate_prior, bool inverse); +RcppExport SEXP _bvartools_post_gamma_state_variance(SEXP aSEXP, SEXP a_initSEXP, SEXP shape_priorSEXP, SEXP rate_priorSEXP, SEXP inverseSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< arma::mat >::type phi(phiSEXP); - Rcpp::traits::input_parameter< arma::vec >::type phi_init(phi_initSEXP); + Rcpp::traits::input_parameter< arma::mat >::type a(aSEXP); + Rcpp::traits::input_parameter< arma::vec >::type a_init(a_initSEXP); Rcpp::traits::input_parameter< arma::vec >::type shape_prior(shape_priorSEXP); Rcpp::traits::input_parameter< arma::vec >::type rate_prior(rate_priorSEXP); - rcpp_result_gen = Rcpp::wrap(post_gamma_state_variance(phi, phi_init, shape_prior, rate_prior)); + Rcpp::traits::input_parameter< bool >::type inverse(inverseSEXP); + rcpp_result_gen = Rcpp::wrap(post_gamma_state_variance(a, a_init, shape_prior, rate_prior, inverse)); return rcpp_result_gen; END_RCPP } @@ -450,7 +451,7 @@ static const R_CallMethodDef CallEntries[] = { {"_bvartools_loglik_normal", (DL_FUNC) &_bvartools_loglik_normal, 2}, {"_bvartools_post_coint_kls", (DL_FUNC) &_bvartools_post_coint_kls, 10}, {"_bvartools_post_coint_kls_sur", (DL_FUNC) &_bvartools_post_coint_kls_sur, 11}, - {"_bvartools_post_gamma_state_variance", (DL_FUNC) &_bvartools_post_gamma_state_variance, 4}, + {"_bvartools_post_gamma_state_variance", (DL_FUNC) &_bvartools_post_gamma_state_variance, 5}, {"_bvartools_post_normal", (DL_FUNC) &_bvartools_post_normal, 5}, {"_bvartools_post_normal_sur", (DL_FUNC) &_bvartools_post_normal_sur, 6}, {"_bvartools_prep_covar_data", (DL_FUNC) &_bvartools_prep_covar_data, 4}, diff --git a/src/post_gamma_state_variance.cpp b/src/post_gamma_state_variance.cpp index 5162d8a..7e420a7 100644 --- a/src/post_gamma_state_variance.cpp +++ b/src/post_gamma_state_variance.cpp @@ -3,15 +3,23 @@ //' Posterior Draws of Error Variances //' -//' Produces a draw of error variances from a gamma posterior density. +//' Produces a draw of the constant diagonal error variance matrix of the +//' state equation of a state space model using an inverse gamma posterior density. //' -//' @param phi a \eqn{K \times T} matrix of time varying parameter draws. -//' @param phi_init a \eqn{K \times 1} vector of initial states. +//' @param a a \eqn{K \times T} matrix of time varying parameter draws. +//' @param a_init a \eqn{K \times 1} vector of initial states. //' @param shape_prior a \eqn{K \times 1} vector of prior shape parameters. //' @param rate_prior a \eqn{K \times 1} vector of prior rate parameters. +//' @param inverse logical. If \code{TRUE}, the function returns the precision matrix, +//' i.e. the inverse of the variance matrix. Defaults to \code{FALSE}. //' -//' @details The function produces a posterior draw of the variaces vector \eqn{a} for the model -//' Follow description in Chan eta al. +//' @details For the state space model with state equation +//' \deqn{a_t = a_{t-1} + v} +//' and measurement equation +//' \deqn{y_t = Z_{t} a_t + u_t} +//' with \eqn{v_t \sim N(0, \Sigma_{v})} and \eqn{u_t \sim N(0, \Sigma_{u,t})} +//' the function produces a draw of the constant diagonal error variances matrix of the +//' state equation \eqn{\Simga_v}. //' //' @references //' Chan, J., Koop, G., Poirier, D. J., & Tobias J. L. (2019). \emph{Bayesian econometric methods} @@ -24,40 +32,43 @@ //' set.seed(1234) # Set RNG seed //' //' # Generate artificial data according to a random walk -//' phi <- matrix(rnorm(k), k, tt + 1) +//' a <- matrix(rnorm(k), k, tt + 1) //' for (i in 2:(tt + 1)) { -//' phi[, i] <- phi[, i - 1] + rnorm(k, 0, sqrt(1 / 100)) +//' a[, i] <- a[, i - 1] + rnorm(k, 0, sqrt(1 / 100)) //' } //' -//' phi_init <- matrix(phi[, 1]) # Define inital state -//' phi <- phi[, -1] # Drop initial state from main sample +//' a_init <- matrix(a[, 1]) # Define inital state +//' a <- a[, -1] # Drop initial state from main sample //' //' # Define priors //' shape_prior <- matrix(1, k) //' rate_prior <- matrix(.0001, k) //' //' # Obtain posterior draw -//' post_gamma_state_variance(phi, phi_init, shape_prior, rate_prior) +//' post_gamma_state_variance(a, a_init, shape_prior, rate_prior) //' //' @return A matrix. //' // [[Rcpp::export]] -arma::mat post_gamma_state_variance(arma::mat phi, arma::vec phi_init, arma::vec shape_prior, arma::vec rate_prior) { +arma::mat post_gamma_state_variance(arma::mat a, arma::vec a_init, arma::vec shape_prior, arma::vec rate_prior, bool inverse = false) { - int k = phi.n_rows; - int tt = phi.n_cols; - arma::mat phi_lag = phi; - phi_lag.col(0) = phi_init; - phi_lag.cols(1, tt - 1) = phi.cols(0, tt - 2); - arma::mat phi_v = arma::trans(phi - phi_lag); - arma::vec psi_sigma_v_post_scale = 1 / (rate_prior + arma::vectorise(arma::sum(arma::pow(phi_v, 2))) * 0.5); - arma::mat psi_sigma_i = arma::zeros(k, k); + int k = a.n_rows; + int tt = a.n_cols; + arma::mat a_lag = a; + a_lag.col(0) = a_init; + a_lag.cols(1, tt - 1) = a.cols(0, tt - 2); + arma::mat a_v = arma::trans(a - a_lag); + arma::vec psi_sigma_v_post_scale = 1 / (rate_prior + arma::vectorise(arma::sum(arma::pow(a_v, 2))) * 0.5); + arma::mat psi_sigma = arma::zeros(k, k); arma::vec shape_post = shape_prior + tt * 0.5; for (int i = 0; i < k; i++) { - psi_sigma_i(i, i) = arma::randg(arma::distr_param(shape_post(i), psi_sigma_v_post_scale(i))); + psi_sigma(i, i) = arma::randg(arma::distr_param(shape_post(i), psi_sigma_v_post_scale(i))); + if (inverse) { + psi_sigma(i, i) = 1 / psi_sigma(i, i); + } } - return psi_sigma_i; + return psi_sigma; } /*** R @@ -67,19 +78,19 @@ tt <- 1000 # Number of observations set.seed(1234) # Set RNG seed # Generate artificial data according to a random walk -phi <- matrix(rnorm(k), k, tt + 1) +a <- matrix(rnorm(k), k, tt + 1) for (i in 2:(tt + 1)) { - phi[, i] <- phi[, i - 1] + rnorm(k, 0, sqrt(1 / 100)) + a[, i] <- a[, i - 1] + rnorm(k, 0, sqrt(1 / 100)) } -phi_init <- matrix(phi[, 1]) # Define inital state -phi <- phi[, -1] # Drop initial state from main sample +a_init <- matrix(a[, 1]) # Define inital state +a <- a[, -1] # Drop initial state from main sample # Define priors shape_prior <- matrix(1, k) rate_prior <- matrix(.0001, k) # Obtain posterior draw -post_gamma_state_variance(phi, phi_init, shape_prior, rate_prior) +post_gamma_state_variance(a, a_init, shape_prior, rate_prior, inverse = FALSE) */ \ No newline at end of file