Add functions for posterior simulation of covariance block

DESCRIPTION

@@ -1,11 +1,11 @@
 Package: bvartools
 Title: Bayesian Inference of Vector Autoregressive and Error Correction Models
-Version: 0.2.3
+Version: 0.2.3.9000
 Date: 2023-08-23
 Authors@R: person(c("Franz", "X."), "Mohr", email = "franz.x.mohr@outlook.com", role = c("aut","cre"), comment = c(ORCiD = "0009-0003-8890-7781"))
 Description: Assists in the set-up of algorithms for Bayesian inference of vector autoregressive (VAR) and error correction (VEC) models. Functions for posterior simulation, forecasting, impulse response analysis and forecast error variance decomposition are largely based on the introductory texts of Chan, Koop, Poirier and Tobias (2019, ISBN: 9781108437493), Koop and Korobilis (2010) <doi:10.1561/0800000013> and Luetkepohl (2006, ISBN: 9783540262398).
 License: GPL (>= 2)
-Depends: R (>= 3.4.0), coda
+Depends: R (>= 3.4.0), coda, Matrix
 Imports: grDevices, graphics, methods, parallel, Rcpp (>= 0.12.14), stats
 LinkingTo: Rcpp, RcppArmadillo
 Encoding: UTF-8
@@ -51,6 +51,8 @@ Collate:
     'plot.bvarprd.R'
     'plot.bvec.R'
     'plot.dfm.R'
+    'post_normal_covar_const.R'
+    'post_normal_covar_tvp.R'
     'predict.bvar.R'
     'summary.bvar.R'
     'print.summary.bvar.R'

NAMESPACE

@@ -42,9 +42,14 @@ export(gen_vec)
 export(inclusion_prior)
 export(irf)
 export(minnesota_prior)
+export(post_normal_covar_const)
+export(post_normal_covar_tvp)
 export(ssvs_prior)
 exportPattern("^[[:alpha:]]+")
 import(methods)
+importFrom(Matrix,chol)
+importFrom(Matrix,solve)
+importFrom(Matrix,t)
 importFrom(Rcpp,sourceCpp)
 importFrom(coda,thin)
 useDynLib(bvartools, .registration = TRUE)

NEWS.md

@@ -1,3 +1,8 @@
+# bvartools 0.2.4
+
+* Added `post_normal_covar_const` for posterior simulation of constant, lower triangular covariance matrices.
+* Added `post_normal_covar_tvp` for posterior simulation of time varying, lower triangular covariance matrices.
+
 # bvartools 0.2.3
 
 * Fixed alias issue resulting from use of `roxygen2`.

R/RcppExports.R

@@ -547,6 +547,10 @@ post_normal_sur <- function(y, z, sigma_i, a_prior, v_i_prior, svd = FALSE) {
     .Call(`_bvartools_post_normal_sur`, y, z, sigma_i, a_prior, v_i_prior, svd)
 }
 
+.prep_covar_data <- function(y, k, tt, tvp) {
+    .Call(`_bvartools_prep_covar_data`, y, k, tt, tvp)
+}
+
 #' Stochastic Volatility
 #' 
 #' Produces a draw of log-volatilities.

R/post_normal_covar_const.R

@@ -0,0 +1,73 @@
+#' Posterior Simulation of Error Covariance Coefficients
+#' 
+#' Produces posterior draws of constant error covariance coefficients.
+#' 
+#' @param y a \eqn{K \times T} matrix of data with \eqn{K} as the number of
+#' endogenous variables and \eqn{T} the number of observations.
+#' @param u_omega_i matrix of error variances of the measurement equation.
+#' Either a \eqn{K \times K} matrix for constant variances or
+#' a \eqn{KT \times KT} matrix for time varying variances.
+#' @param prior_mean vector of prior means. In case of TVP, this vector is used
+#' as initial condition.
+#' @param prior_covariance_i inverse prior covariance matrix. In case of TVP, this matrix
+#' is used as initial condition.
+#' 
+#' @details For the multivariate model \eqn{A_0 y_t = u_t} with \eqn{u_t \sim N(0, \Omega_t)}
+#' the function produces a draw of the lower triangular part of \eqn{A_0} similar as in
+#' Primiceri (2005), i.e., using \deqn{y_t = Z_t \psi + u_t,}
+#' where 
+#' \deqn{Z_{t} = \begin{bmatrix} 0 & \dotsm & \dotsm & 0 \\ -y_{1, t} & 0 & \dotsm & 0 \\ 0 & -y_{[1,2], t} & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \dotsm & 0 & -y_{[1,...,K-1], t} \end{bmatrix}}
+#' and \eqn{y_{[1,...,K-1], t}} denotes the first to \eqn{(K-1)}th elements of the vector \eqn{y_t}.
+#' 
+#' @references Primiceri, G. E. (2005). Time varying structural vector autoregressions and monetary policy.
+#' \emph{The Review of Economic Studies, 72}(3), 821--852. \doi{10.1111/j.1467-937X.2005.00353.x}
+#' 
+#' @returns A matrix.
+#' 
+#' @examples
+#' # Load example data
+#' data("e1")
+#' y <- log(t(e1))
+#' 
+#' # Generate artificial draws of other matrices
+#' u_omega_i <- diag(1, 3)
+#' prior_mean <- matrix(0, 3)
+#' prior_covariance_i <- diag(0, 3)
+#' 
+#' # Obtain posterior draw
+#' post_normal_covar_const(y, u_omega_i, prior_mean, prior_covariance_i)
+#'
+#' @export
+post_normal_covar_const <- function(y, u_omega_i, prior_mean, prior_covariance_i) {
+
+  k <- nrow(y)
+  if (k == 1L) {
+    stop("Argument 'y' must contain at least two variables.")
+  }
+  n_covar <- k * (k - 1) / 2
+  tt <- ncol(y)
+  y <- matrix(y)
+
+  # Generate z for lower triangular design
+  # Use C++ to speed up the for loop
+  z <- .prep_covar_data(y, k, tt, FALSE)
+
+  # Get positions of values, on which variables in z are regressed
+  pos_used <- rep(k * 0:(tt - 1), each = k - 1) + 2:k
+
+  # Trim endogenous variables
+  y <- matrix(y[pos_used, ])
+
+  # Adapt error variance matrix
+  if (NCOL(u_omega_i) == k & NCOL(u_omega_i) == k) {
+    u_omega_i <- kronecker(Diagonal(tt, 1), u_omega_i)
+  }
+  # Trim error variance matrix
+  u_omega_i <- u_omega_i[pos_used, pos_used]
+
+  # Draw coefficients
+  v_post <- prior_covariance_i + crossprod(z, u_omega_i) %*% z
+  mu_post <- solve(v_post, prior_covariance_i %*% prior_mean + crossprod(z, u_omega_i) %*% y)
+
+  return(matrix(mu_post + solve(chol(v_post), matrix(rnorm(n_covar)))))
+}

R/post_normal_covar_tvp.R

@@ -0,0 +1,91 @@
+#' Posterior Simulation of Error Covariance Coefficients
+#' 
+#' Produces posterior draws of time varying error covariance coefficients.
+#' 
+#' @param y a \eqn{K \times T} matrix of data with \eqn{K} as the number of
+#' endogenous variables and \eqn{T} the number of observations.
+#' @param u_omega_i matrix of error variances of the measurement equation.
+#' Either a \eqn{K \times K} matrix for constant variances or
+#' a \eqn{KT \times KT} matrix for time varying variances.
+#' @param v_sigma_i matrix of error variances of the state equation.
+#' Either an \eqn{M \times M} matrix for constant variances or
+#' an \eqn{MT \times MT} matrix for time varying variances, where \eqn{M} is the
+#' number of estimated variables.
+#' @param psi_init a vector of inital values of the state equation.
+#' 
+#' @details For the multivariate model \eqn{A_{0,t} y_t = u_t} with \eqn{u_t \sim N(0, \Omega_t)}
+#' the function produces a draw of the lower triangular part of \eqn{A_{0,t}} similar as in
+#' Primiceri (2005), i.e., using \deqn{y_t = Z_t \psi_t + u_t,}
+#' where 
+#' \deqn{Z_{t} = \begin{bmatrix} 0 & \dotsm & \dotsm & 0 \\ -y_{1, t} & 0 & \dotsm & 0 \\ 0 & -y_{[1,2], t} & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \dotsm & 0 & -y_{[1,...,K-1], t} \end{bmatrix}}
+#' and \eqn{y_{[1,...,K-1], t}} denotes the first to \eqn{(K-1)}th elements of the vector \eqn{y_t}.
+#' 
+#' The algorithm of Chan and Jeliazkov (2009) is used to obtain time varying coefficients.
+#' 
+#' @references
+#' 
+#' Chan, J., & Jeliazkov, I. (2009). Efficient simulation and integrated likelihood estimation in state space
+#' models. \emph{International Journal of Mathematical Modelling and Numerical Optimisation, 1}(1/2), 101–120.
+#' \doi{10.1504/IJMMNO.2009.030090}
+#' 
+#' Primiceri, G. E. (2005). Time varying structural vector autoregressions and monetary policy.
+#' \emph{The Review of Economic Studies, 72}(3), 821--852. \doi{10.1111/j.1467-937X.2005.00353.x}
+#' 
+#' @returns A matrix.
+#' 
+#' @examples
+#' # Load example data
+#' data("e1")
+#' y <- log(t(e1))
+#' 
+#' # Generate artificial draws of other matrices
+#' u_omega_i <- diag(1, 3)
+#' v_sigma_i <- diag(1000, 3)
+#' psi_init <- matrix(0, 3)
+#' 
+#' # Obtain posterior draw
+#' post_normal_covar_tvp(y, u_omega_i, v_sigma_i, psi_init)
+#'
+#' @export
+post_normal_covar_tvp <- function(y, u_omega_i, v_sigma_i, psi_init) {
+
+  k <- nrow(y)
+  if (k == 1L) {
+    stop("Argument 'y' must contain at least two variables.")
+  }
+  n_covar <- k * (k - 1) / 2
+  tt <- ncol(y)
+  y <- matrix(y)
+
+  # Generate z for lower triangular design
+  z <- .prep_covar_data(y, k, tt, TRUE)
+
+  # Get positions of values, on which variables in z are regressed
+  pos_used <- rep(k * 0:(tt - 1), each = k - 1) + 2:k
+
+  # Trim endogenous variables
+  y <- matrix(y[pos_used, ])
+
+  # Trim 
+  if (NCOL(u_omega_i) == k & NCOL(u_omega_i) == k) {
+    u_omega_i <- kronecker(Diagonal(tt, 1), u_omega_i)
+  }
+  # Trim error variance matrix
+  u_omega_i <- u_omega_i[pos_used, pos_used]
+
+  # Draw coefficients
+  hh <- Matrix(0, n_covar * tt, n_covar * tt)
+  diag(hh) <- 1
+  diag(hh[-(1:n_covar), -(n_covar * (tt - 1) + 1:n_covar)]) <- -1
+
+  if (NCOL(v_sigma_i) == n_covar & NCOL(v_sigma_i) == n_covar) {
+    v_sigma_i <- kronecker(Diagonal(tt, 1), v_sigma_i)
+  }
+
+  hh_temp <- t(hh) %*% v_sigma_i %*% hh
+  x_temp <- t(z) %*% u_omega_i
+  v_post <- hh_temp + x_temp %*% z
+  mu_post <- solve(v_post, hh_temp %*% kronecker(matrix(1, tt), psi_init) + x_temp %*% y)
+
+  return(matrix(mu_post + solve(chol(v_post), matrix(rnorm(n_covar * tt)))))
+}

R/zzz.R

@@ -1,3 +1,9 @@
+# Load needed functions from different packages for fast access
+
+#' @importFrom Matrix chol
+#' @importFrom Matrix solve
+#' @importFrom Matrix t
+
 # Unload the DLL when the package is unloaded
 .onUnload <- function (libpath) {
   library.dynam.unload("bvartools", libpath)

src/RcppExports.cpp

@@ -243,6 +243,20 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// prep_covar_data
+arma::sp_mat prep_covar_data(arma::vec y, int k, int tt, bool tvp);
+RcppExport SEXP _bvartools_prep_covar_data(SEXP ySEXP, SEXP kSEXP, SEXP ttSEXP, SEXP tvpSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< arma::vec >::type y(ySEXP);
+    Rcpp::traits::input_parameter< int >::type k(kSEXP);
+    Rcpp::traits::input_parameter< int >::type tt(ttSEXP);
+    Rcpp::traits::input_parameter< bool >::type tvp(tvpSEXP);
+    rcpp_result_gen = Rcpp::wrap(prep_covar_data(y, k, tt, tvp));
+    return rcpp_result_gen;
+END_RCPP
+}
 // stoch_vol
 arma::vec stoch_vol(arma::vec y, arma::vec h, double sigma, double h_init, double constant);
 static SEXP _bvartools_stoch_vol_try(SEXP ySEXP, SEXP hSEXP, SEXP sigmaSEXP, SEXP h_initSEXP, SEXP constantSEXP) {
@@ -424,6 +438,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_bvartools_post_coint_kls_sur", (DL_FUNC) &_bvartools_post_coint_kls_sur, 11},
     {"_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},
     {"_bvartools_stoch_vol", (DL_FUNC) &_bvartools_stoch_vol, 5},
     {"_bvartools_stochvol_ksc1998", (DL_FUNC) &_bvartools_stochvol_ksc1998, 5},
     {"_bvartools_stochvol_ocsn2007", (DL_FUNC) &_bvartools_stochvol_ocsn2007, 5},