diff --git a/DESCRIPTION b/DESCRIPTION index 474bf89..aeb783c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -20,4 +20,4 @@ Suggests: tinytest LinkingTo: Rcpp, RcppProgress, RcppArmadillo, RcppTN BugReports: https://github.com/donotdespair/bsvars/issues -RoxygenNote: 7.2.1 +RoxygenNote: 7.2.2 diff --git a/NEWS.md b/NEWS.md index 38f830f..ef9ef38 100644 --- a/NEWS.md +++ b/NEWS.md @@ -15,7 +15,7 @@ 6. Include citation info for the package [#12](https://github.com/donotdespair/bsvars/issues/12) 7. Corrected sampler for AR parameter of the SV equations [#19](https://github.com/donotdespair/bsvars/issues/19) 8. Added samplers from joint predictive densities [#15](https://github.com/donotdespair/bsvars/issues/15) - +9. A new centred Stochastic Volatility heteroskedastic process is implemented [#22](https://github.com/donotdespair/bsvars/issues/22) # bsvars 1.0.0 diff --git a/R/bsvars-package.R b/R/bsvars-package.R index 14d32a8..abf0d07 100644 --- a/R/bsvars-package.R +++ b/R/bsvars-package.R @@ -56,7 +56,8 @@ #' \itemize{ #' \item homoskedastic model with unit variances #' \item heteroskedastic model with stationary Markov switching in the variances -#' \item heteroskedastic model with Stochastic Volatility process for variances +#' \item heteroskedastic model with non-centred Stochastic Volatility process for variances +#' \item heteroskedastic model with centred Stochastic Volatility process for variances #' \item non-normal model with a finite mixture of normal components and component-specific variances #' \item heteroskedastic model with sparse Markov switching in the variances where the number of heteroskedastic components is estimated #' \item non-normal model with a sparse mixture of normal components and component-specific variances where the number of heteroskedastic components is estimated diff --git a/R/estimate.BSVARSV.R b/R/estimate.BSVARSV.R index 5e6e182..110fa26 100644 --- a/R/estimate.BSVARSV.R +++ b/R/estimate.BSVARSV.R @@ -2,13 +2,17 @@ #' @title Bayesian estimation of a Structural Vector Autoregression with #' Stochastic Volatility heteroskedasticity via Gibbs sampler #' -#' @description Estimates the SVAR with Stochastic Volatility (SV) heteroskedasticity proposed by Lütkepohl, Shang, Uzeda, and Woźniak (2022). +#' @description Estimates the SVAR with Stochastic Volatility (SV) heteroskedasticity +#' proposed by Lütkepohl, Shang, Uzeda, and Woźniak (2022). #' Implements the Gibbs sampler proposed by Waggoner & Zha (2003) -#' for the structural matrix \eqn{B} and the equation-by-equation sampler by Chan, Koop, & Yu (2021) -#' for the autoregressive slope parameters \eqn{A}. Additionally, the parameter matrices \eqn{A} and \eqn{B} -#' follow a Minnesota prior and generalised-normal prior distributions respectively with the matrix-specific -#' overall shrinkage parameters estimated thanks to a 3-level hierarchical prior distribution. The SV model -#' is estimated in a non-centred parameterisation using a range of techniques including: +#' for the structural matrix \eqn{B} and the equation-by-equation sampler +#' by Chan, Koop, & Yu (2021) +#' for the autoregressive slope parameters \eqn{A}. Additionally, +#' the parameter matrices \eqn{A} and \eqn{B} +#' follow a Minnesota prior and generalised-normal prior distributions +#' respectively with the matrix-specific +#' overall shrinkage parameters estimated thanks to a 3-level hierarchical prior distribution. +#' The SV model is estimated using a range of techniques including: #' simulation smoother, auxiliary mixture, ancillarity-sufficiency interweaving strategy, #' and generalised inverse Gaussian distribution summarised by Kastner & Frühwirth-Schnatter (2014). #' See section \bold{Details} for the model equations. @@ -23,15 +27,30 @@ #' \deqn{BE = U} #' where \eqn{U} is an \code{NxT} matrix of structural form error terms, and #' \eqn{B} is an \code{NxN} matrix of contemporaneous relationships. -#' #' Finally, the structural shocks, \eqn{U}, are temporally and contemporaneously independent and jointly normally distributed with zero mean. -#' The conditional variance of the \code{n}th shock at time \code{t} is given by: +#' +#' Two alternative specifications of the conditional variance of the \code{n}th shock at time \code{t} +#' can be estimated: non-centred Stochastic Volatility by Lütkepohl, Shang, Uzeda, and Woźniak (2022) +#' or centred Stochastic Volatility by Chan, Koop, & Yu (2021). +#' +#' The non-centred Stochastic Volatility by Lütkepohl, Shang, Uzeda, and Woźniak (2022) +#' is selected by setting argument \code{centred_sv} of function \code{specify_bsvar_sv$new()} to value \code{FALSE}. +#' It has the conditional variances given by: #' \deqn{Var_{t-1}[u_{n.t}] = exp(w_n h_{n.t})} #' where \eqn{w_n} is the estimated conditional standard deviation of the log-conditional variance #' and the log-volatility process \eqn{h_{n.t}} follows an autoregressive process: #' \deqn{h_{n.t} = g_n h_{n.t-1} + v_{n.t}} #' where \eqn{h_{n.0}=0}, \eqn{g_n} is an autoregressive parameter and \eqn{v_{n.t}} is a standard normal error term. #' +#' The centred Stochastic Volatility by Chan, Koop, & Yu (2021) +#' is selected by setting argument \code{centred_sv} of function \code{specify_bsvar_sv$new()} to value \code{TRUE}. +#' Its conditional variances are given by: +#' \deqn{Var_{t-1}[u_{n.t}] = exp(h_{n.t})} +#' where the log-conditional variances \eqn{h_{n.t}} follow an autoregressive process: +#' \deqn{h_{n.t} = g_n h_{n.t-1} + v_{n.t}} +#' where \eqn{h_{n.0}=0}, \eqn{g_n} is an autoregressive parameter and \eqn{v_{n.t}} is a zero-mean normal error term +#' with variance \eqn{s_{v.n}^2}. +#' #' @param specification an object of class BSVARSV generated using the \code{specify_bsvar_sv$new()} function. #' @param S a positive integer, the number of posterior draws to be generated #' @param thin a positive integer, specifying the frequency of MCMC output thinning @@ -112,9 +131,10 @@ estimate.BSVARSV <- function(specification, S, thin = 10, show_progress = TRUE) starting_values = specification$starting_values$get_starting_values() VB = specification$identification$get_identification() data_matrices = specification$data_matrices$get_data_matrices() + centred_sv = specification$centred_sv # estimation - qqq = .Call(`_bsvars_bsvar_sv_cpp`, S, data_matrices$Y, data_matrices$X, prior, VB, starting_values, thin, TRUE, show_progress) + qqq = .Call(`_bsvars_bsvar_sv_cpp`, S, data_matrices$Y, data_matrices$X, prior, VB, starting_values, thin, centred_sv, show_progress) specification$starting_values$set_starting_values(qqq$last_draw) output = specify_posterior_bsvar_sv$new(specification, qqq$posterior) @@ -171,9 +191,10 @@ estimate.PosteriorBSVARSV <- function(specification, S, thin = 10, show_progress starting_values = specification$last_draw$starting_values$get_starting_values() VB = specification$last_draw$identification$get_identification() data_matrices = specification$last_draw$data_matrices$get_data_matrices() + centred_sv = specification$last_draw$centred_sv # estimation - qqq = .Call(`_bsvars_bsvar_sv_cpp`, S, data_matrices$Y, data_matrices$X, prior, VB, starting_values, thin, TRUE, show_progress) + qqq = .Call(`_bsvars_bsvar_sv_cpp`, S, data_matrices$Y, data_matrices$X, prior, VB, starting_values, thin, centred_sv, show_progress) specification$last_draw$starting_values$set_starting_values(qqq$last_draw) output = specify_posterior_bsvar_sv$new(specification$last_draw, qqq$posterior) diff --git a/R/forecast.R b/R/forecast.R index a0f7792..f566941 100644 --- a/R/forecast.R +++ b/R/forecast.R @@ -227,8 +227,9 @@ forecast.PosteriorBSVARSV = function(posterior, horizon) { T = ncol(posterior$last_draw$data_matrices$X) X_T = posterior$last_draw$data_matrices$X[,T] posterior_h_T = posterior$posterior$h[,T,] + centred_sv = posterior$last_draw$centred_sv - fore = .Call(`_bsvars_forecast_bsvar_sv`, posterior_B, posterior_A, posterior_h_T, posterior_rho, posterior_omega, X_T, horizon) + fore = .Call(`_bsvars_forecast_bsvar_sv`, posterior_B, posterior_A, posterior_h_T, posterior_rho, posterior_omega, X_T, horizon, centred_sv) class(fore) = "Forecasts" return(fore) diff --git a/R/specify_bsvar_sv.R b/R/specify_bsvar_sv.R index f3b64d6..d4c2aa0 100644 --- a/R/specify_bsvar_sv.R +++ b/R/specify_bsvar_sv.R @@ -125,6 +125,9 @@ specify_starting_values_bsvar_sv = R6::R6Class( #' @field omega an \code{N}-vector with values of SV process conditional standard deviations. omega = numeric(), + #' @field sigma2v an \code{N}-vector with values of SV process conditional variances. + sigma2v = numeric(), + #' @field S an \code{NxT} integer matrix with the auxiliary mixture component indicators. S = matrix(), @@ -150,6 +153,7 @@ specify_starting_values_bsvar_sv = R6::R6Class( self$h = matrix(rnorm(N * T, sd = .01), N, T) self$rho = rep(.5, N) self$omega = rep(.1, N) + self$sigma2v = rep(.1^2, N) self$S = matrix(1, N, T) self$sigma2_omega = rep(1, N) self$s_ = rep(0.05, N) @@ -171,6 +175,7 @@ specify_starting_values_bsvar_sv = R6::R6Class( h = self$h, rho = self$rho, omega = self$omega, + sigma2v = self$sigma2v, S = self$S, sigma2_omega = self$sigma2_omega, s_ = self$s_ @@ -198,6 +203,7 @@ specify_starting_values_bsvar_sv = R6::R6Class( self$h = last_draw$h self$rho = last_draw$rho self$omega = last_draw$omega + self$sigma2v = last_draw$sigma2v self$S = last_draw$S self$sigma2_omega = last_draw$sigma2_omega self$s_ = last_draw$s_ @@ -242,17 +248,22 @@ specify_bsvar_sv = R6::R6Class( #' @field starting_values an object StartingValuesBSVARSV with the starting values. starting_values = list(), + #' @field centred_sv a logical value - if true a centred parameterisation of the Stochastic Volatility process is estimated. Otherwise, its non-centred parameterisation is estimated. See Lütkepohl, Shang, Uzeda, Woźniak (2022) for more info. + centred_sv = logical(), + #' @description #' Create a new specification of the BSVAR model with Stochastic Volatility heteroskedasticity, BSVARSV. #' @param data a \code{(T+p)xN} matrix with time series data. #' @param p a positive integer providing model's autoregressive lag order. #' @param B a logical \code{NxN} matrix containing value \code{TRUE} for the elements of the structural matrix \eqn{B} to be estimated and value \code{FALSE} for exclusion restrictions to be set to zero. + #' @param centred_sv a logical value. If \code{FALSE} a non-centred Stochastic Volatility processes for conditional variances are estimated. Otherwise, a centred process is estimated. #' @param stationary an \code{N} logical vector - its element set to \code{FALSE} sets the prior mean for the autoregressive parameters of the \code{N}th equation to the white noise process, otherwise to random walk. #' @return A new complete specification for the bsvar model with Stochastic Volatility heteroskedasticity, BSVARSV. initialize = function( data, p = 1L, B, + centred_sv = FALSE, stationary = rep(FALSE, ncol(data)) ) { stopifnot("Argument p has to be a positive integer." = ((p %% 1) == 0 & p > 0)) @@ -273,6 +284,7 @@ specify_bsvar_sv = R6::R6Class( self$identification = specify_identification_bsvars$new(N, B) self$prior = specify_prior_bsvar_sv$new(N, p, stationary) self$starting_values = specify_starting_values_bsvar_sv$new(N, self$p, T) + self$centred_sv = centred_sv }, # END initialize #' @description diff --git a/inst/include/bsvars_RcppExports.h b/inst/include/bsvars_RcppExports.h index c1c42ea..6e13a7c 100644 --- a/inst/include/bsvars_RcppExports.h +++ b/inst/include/bsvars_RcppExports.h @@ -193,7 +193,7 @@ namespace bsvars { return Rcpp::as(rcpp_result_gen); } - inline Rcpp::List bsvar_sv_cpp(const int& S, const arma::mat& Y, const arma::mat& X, const Rcpp::List& prior, const arma::field& VB, const Rcpp::List& starting_values, const int thin = 100, const bool sample_s_ = true, const bool show_progress = true) { + inline Rcpp::List bsvar_sv_cpp(const int& S, const arma::mat& Y, const arma::mat& X, const Rcpp::List& prior, const arma::field& VB, const Rcpp::List& starting_values, const int thin = 100, const bool centred_sv = false, const bool show_progress = true) { typedef SEXP(*Ptr_bsvar_sv_cpp)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP); static Ptr_bsvar_sv_cpp p_bsvar_sv_cpp = NULL; if (p_bsvar_sv_cpp == NULL) { @@ -203,7 +203,7 @@ namespace bsvars { RObject rcpp_result_gen; { RNGScope RCPP_rngScope_gen; - rcpp_result_gen = p_bsvar_sv_cpp(Shield(Rcpp::wrap(S)), Shield(Rcpp::wrap(Y)), Shield(Rcpp::wrap(X)), Shield(Rcpp::wrap(prior)), Shield(Rcpp::wrap(VB)), Shield(Rcpp::wrap(starting_values)), Shield(Rcpp::wrap(thin)), Shield(Rcpp::wrap(sample_s_)), Shield(Rcpp::wrap(show_progress))); + rcpp_result_gen = p_bsvar_sv_cpp(Shield(Rcpp::wrap(S)), Shield(Rcpp::wrap(Y)), Shield(Rcpp::wrap(X)), Shield(Rcpp::wrap(prior)), Shield(Rcpp::wrap(VB)), Shield(Rcpp::wrap(starting_values)), Shield(Rcpp::wrap(thin)), Shield(Rcpp::wrap(centred_sv)), Shield(Rcpp::wrap(show_progress))); } if (rcpp_result_gen.inherits("interrupted-error")) throw Rcpp::internal::InterruptedException(); @@ -256,17 +256,17 @@ namespace bsvars { return Rcpp::as(rcpp_result_gen); } - inline Rcpp::List forecast_bsvar_sv(arma::cube& posterior_B, arma::cube& posterior_A, arma::vec& posterior_h_T, arma::mat& posterior_rho, arma::mat& posterior_omega, arma::vec& X_T, const int horizon) { - typedef SEXP(*Ptr_forecast_bsvar_sv)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP); + inline Rcpp::List forecast_bsvar_sv(arma::cube& posterior_B, arma::cube& posterior_A, arma::vec& posterior_h_T, arma::mat& posterior_rho, arma::mat& posterior_omega, arma::vec& X_T, const int horizon, const bool centred_sv = FALSE) { + typedef SEXP(*Ptr_forecast_bsvar_sv)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP); static Ptr_forecast_bsvar_sv p_forecast_bsvar_sv = NULL; if (p_forecast_bsvar_sv == NULL) { - validateSignature("Rcpp::List(*forecast_bsvar_sv)(arma::cube&,arma::cube&,arma::vec&,arma::mat&,arma::mat&,arma::vec&,const int)"); + validateSignature("Rcpp::List(*forecast_bsvar_sv)(arma::cube&,arma::cube&,arma::vec&,arma::mat&,arma::mat&,arma::vec&,const int,const bool)"); p_forecast_bsvar_sv = (Ptr_forecast_bsvar_sv)R_GetCCallable("bsvars", "_bsvars_forecast_bsvar_sv"); } RObject rcpp_result_gen; { RNGScope RCPP_rngScope_gen; - rcpp_result_gen = p_forecast_bsvar_sv(Shield(Rcpp::wrap(posterior_B)), Shield(Rcpp::wrap(posterior_A)), Shield(Rcpp::wrap(posterior_h_T)), Shield(Rcpp::wrap(posterior_rho)), Shield(Rcpp::wrap(posterior_omega)), Shield(Rcpp::wrap(X_T)), Shield(Rcpp::wrap(horizon))); + rcpp_result_gen = p_forecast_bsvar_sv(Shield(Rcpp::wrap(posterior_B)), Shield(Rcpp::wrap(posterior_A)), Shield(Rcpp::wrap(posterior_h_T)), Shield(Rcpp::wrap(posterior_rho)), Shield(Rcpp::wrap(posterior_omega)), Shield(Rcpp::wrap(X_T)), Shield(Rcpp::wrap(horizon)), Shield(Rcpp::wrap(centred_sv))); } if (rcpp_result_gen.inherits("interrupted-error")) throw Rcpp::internal::InterruptedException(); @@ -793,17 +793,38 @@ namespace bsvars { return Rcpp::as(rcpp_result_gen); } - inline Rcpp::List svar_nc1(arma::rowvec& aux_h_n, double& aux_rho_n, double& aux_omega_n, double& aux_sigma2_omega_n, double& aux_s_n, arma::urowvec& aux_S_n, const arma::rowvec& u, const Rcpp::List& prior, bool sample_s_ = true) { - typedef SEXP(*Ptr_svar_nc1)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP); + inline Rcpp::List svar_nc1(arma::rowvec& aux_h_n, double& aux_rho_n, double& aux_omega_n, double& aux_sigma2v_n, double& aux_sigma2_omega_n, double& aux_s_n, arma::urowvec& aux_S_n, const arma::rowvec& u, const Rcpp::List& prior, bool sample_s_ = true) { + typedef SEXP(*Ptr_svar_nc1)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP); static Ptr_svar_nc1 p_svar_nc1 = NULL; if (p_svar_nc1 == NULL) { - validateSignature("Rcpp::List(*svar_nc1)(arma::rowvec&,double&,double&,double&,double&,arma::urowvec&,const arma::rowvec&,const Rcpp::List&,bool)"); + validateSignature("Rcpp::List(*svar_nc1)(arma::rowvec&,double&,double&,double&,double&,double&,arma::urowvec&,const arma::rowvec&,const Rcpp::List&,bool)"); p_svar_nc1 = (Ptr_svar_nc1)R_GetCCallable("bsvars", "_bsvars_svar_nc1"); } RObject rcpp_result_gen; { RNGScope RCPP_rngScope_gen; - rcpp_result_gen = p_svar_nc1(Shield(Rcpp::wrap(aux_h_n)), Shield(Rcpp::wrap(aux_rho_n)), Shield(Rcpp::wrap(aux_omega_n)), Shield(Rcpp::wrap(aux_sigma2_omega_n)), Shield(Rcpp::wrap(aux_s_n)), Shield(Rcpp::wrap(aux_S_n)), Shield(Rcpp::wrap(u)), Shield(Rcpp::wrap(prior)), Shield(Rcpp::wrap(sample_s_))); + rcpp_result_gen = p_svar_nc1(Shield(Rcpp::wrap(aux_h_n)), Shield(Rcpp::wrap(aux_rho_n)), Shield(Rcpp::wrap(aux_omega_n)), Shield(Rcpp::wrap(aux_sigma2v_n)), Shield(Rcpp::wrap(aux_sigma2_omega_n)), Shield(Rcpp::wrap(aux_s_n)), Shield(Rcpp::wrap(aux_S_n)), Shield(Rcpp::wrap(u)), Shield(Rcpp::wrap(prior)), Shield(Rcpp::wrap(sample_s_))); + } + if (rcpp_result_gen.inherits("interrupted-error")) + throw Rcpp::internal::InterruptedException(); + if (Rcpp::internal::isLongjumpSentinel(rcpp_result_gen)) + 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); + } + + inline Rcpp::List svar_ce1(arma::rowvec& aux_h_n, double& aux_rho_n, double& aux_omega_n, double& aux_sigma2v_n, double& aux_sigma2_omega_n, double& aux_s_n, arma::urowvec& aux_S_n, const arma::rowvec& u, const Rcpp::List& prior, bool sample_s_ = true) { + typedef SEXP(*Ptr_svar_ce1)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP); + static Ptr_svar_ce1 p_svar_ce1 = NULL; + if (p_svar_ce1 == NULL) { + validateSignature("Rcpp::List(*svar_ce1)(arma::rowvec&,double&,double&,double&,double&,double&,arma::urowvec&,const arma::rowvec&,const Rcpp::List&,bool)"); + p_svar_ce1 = (Ptr_svar_ce1)R_GetCCallable("bsvars", "_bsvars_svar_ce1"); + } + RObject rcpp_result_gen; + { + RNGScope RCPP_rngScope_gen; + rcpp_result_gen = p_svar_ce1(Shield(Rcpp::wrap(aux_h_n)), Shield(Rcpp::wrap(aux_rho_n)), Shield(Rcpp::wrap(aux_omega_n)), Shield(Rcpp::wrap(aux_sigma2v_n)), Shield(Rcpp::wrap(aux_sigma2_omega_n)), Shield(Rcpp::wrap(aux_s_n)), Shield(Rcpp::wrap(aux_S_n)), Shield(Rcpp::wrap(u)), Shield(Rcpp::wrap(prior)), Shield(Rcpp::wrap(sample_s_))); } if (rcpp_result_gen.inherits("interrupted-error")) throw Rcpp::internal::InterruptedException(); diff --git a/inst/tinytest/test_estimate_bsvar_sv.R b/inst/tinytest/test_estimate_bsvar_sv.R index acc6509..379a377 100644 --- a/inst/tinytest/test_estimate_bsvar_sv.R +++ b/inst/tinytest/test_estimate_bsvar_sv.R @@ -35,3 +35,22 @@ expect_identical( run_no3$last_draw$starting_values$B[1,1], info = "estimate_bsvar_sv: the last_draw(s) of a normal and pipe run to be identical." ) + + +set.seed(1) +suppressMessages( + specification_no1 <- specify_bsvar_sv$new(us_fiscal_lsuw, centred_sv = TRUE) +) +run_no1 <- estimate(specification_no1, 3, 1, show_progress = FALSE) + +set.seed(1) +suppressMessages( + specification_no2 <- specify_bsvar_sv$new(us_fiscal_lsuw, centred_sv = TRUE) +) +run_no2 <- estimate(specification_no2, 3, 1, show_progress = FALSE) + +expect_identical( + run_no1$last_draw$starting_values$B[1,1], + run_no2$last_draw$starting_values$B[1,1], + info = "estimate_bsvar_sv centred: the last_draw(s) of two runs to be identical." +) diff --git a/inst/tinytest/test_forecast.R b/inst/tinytest/test_forecast.R index 236fb43..3a15f28 100644 --- a/inst/tinytest/test_forecast.R +++ b/inst/tinytest/test_forecast.R @@ -188,3 +188,31 @@ expect_error( info = "forecast: sv: specify horizon." ) + +# for bsvar_sv centred +set.seed(1) +suppressMessages( + specification_no1 <- specify_bsvar_sv$new(us_fiscal_lsuw, centred_sv = TRUE) +) +run_no1 <- estimate(specification_no1, 3, 1, show_progress = FALSE) +ff <- forecast(run_no1, horizon = 2) + +set.seed(1) +suppressMessages( + ff2 <- us_fiscal_lsuw |> + specify_bsvar_sv$new(centred_sv = TRUE) |> + estimate(S = 3, thin = 1, show_progress = FALSE) |> + forecast(horizon = 2) +) + + +expect_identical( + ff$forecasts[1,1,1], ff2$forecasts[1,1,1], + info = "forecast: sv centred: forecast identical for normal and pipe workflow." +) + +expect_identical( + ff$forecasts_sigma[1,1,1], ff2$forecasts_sigma[1,1,1], + info = "forecast: sv centred: sigma forecast identical for normal and pipe workflow." +) + diff --git a/man/bsvars-package.Rd b/man/bsvars-package.Rd index c63d14c..3622c02 100644 --- a/man/bsvars-package.Rd +++ b/man/bsvars-package.Rd @@ -41,7 +41,8 @@ variances. The different models include: \itemize{ \item homoskedastic model with unit variances \item heteroskedastic model with stationary Markov switching in the variances - \item heteroskedastic model with Stochastic Volatility process for variances + \item heteroskedastic model with non-centred Stochastic Volatility process for variances + \item heteroskedastic model with centred Stochastic Volatility process for variances \item non-normal model with a finite mixture of normal components and component-specific variances \item heteroskedastic model with sparse Markov switching in the variances where the number of heteroskedastic components is estimated \item non-normal model with a sparse mixture of normal components and component-specific variances where the number of heteroskedastic components is estimated diff --git a/man/estimate.BSVARSV.Rd b/man/estimate.BSVARSV.Rd index 713164a..edfe1c3 100644 --- a/man/estimate.BSVARSV.Rd +++ b/man/estimate.BSVARSV.Rd @@ -35,13 +35,17 @@ An object of class PosteriorBSVARSV containing the Bayesian estimation output an \code{last_draw} an object of class BSVARSV with the last draw of the current MCMC run as the starting value to be passed to the continuation of the MCMC estimation using \code{estimate()}. } \description{ -Estimates the SVAR with Stochastic Volatility (SV) heteroskedasticity proposed by Lütkepohl, Shang, Uzeda, and Woźniak (2022). +Estimates the SVAR with Stochastic Volatility (SV) heteroskedasticity +proposed by Lütkepohl, Shang, Uzeda, and Woźniak (2022). Implements the Gibbs sampler proposed by Waggoner & Zha (2003) -for the structural matrix \eqn{B} and the equation-by-equation sampler by Chan, Koop, & Yu (2021) -for the autoregressive slope parameters \eqn{A}. Additionally, the parameter matrices \eqn{A} and \eqn{B} -follow a Minnesota prior and generalised-normal prior distributions respectively with the matrix-specific -overall shrinkage parameters estimated thanks to a 3-level hierarchical prior distribution. The SV model -is estimated in a non-centred parameterisation using a range of techniques including: +for the structural matrix \eqn{B} and the equation-by-equation sampler +by Chan, Koop, & Yu (2021) +for the autoregressive slope parameters \eqn{A}. Additionally, +the parameter matrices \eqn{A} and \eqn{B} +follow a Minnesota prior and generalised-normal prior distributions +respectively with the matrix-specific +overall shrinkage parameters estimated thanks to a 3-level hierarchical prior distribution. +The SV model is estimated using a range of techniques including: simulation smoother, auxiliary mixture, ancillarity-sufficiency interweaving strategy, and generalised inverse Gaussian distribution summarised by Kastner & Frühwirth-Schnatter (2014). See section \bold{Details} for the model equations. @@ -56,14 +60,29 @@ The structural equation is given by \deqn{BE = U} where \eqn{U} is an \code{NxT} matrix of structural form error terms, and \eqn{B} is an \code{NxN} matrix of contemporaneous relationships. - Finally, the structural shocks, \eqn{U}, are temporally and contemporaneously independent and jointly normally distributed with zero mean. -The conditional variance of the \code{n}th shock at time \code{t} is given by: + +Two alternative specifications of the conditional variance of the \code{n}th shock at time \code{t} +can be estimated: non-centred Stochastic Volatility by Lütkepohl, Shang, Uzeda, and Woźniak (2022) +or centred Stochastic Volatility by Chan, Koop, & Yu (2021). + +The non-centred Stochastic Volatility by Lütkepohl, Shang, Uzeda, and Woźniak (2022) +is selected by setting argument \code{centred_sv} of function \code{specify_bsvar_sv$new()} to value \code{FALSE}. +It has the conditional variances given by: \deqn{Var_{t-1}[u_{n.t}] = exp(w_n h_{n.t})} where \eqn{w_n} is the estimated conditional standard deviation of the log-conditional variance and the log-volatility process \eqn{h_{n.t}} follows an autoregressive process: \deqn{h_{n.t} = g_n h_{n.t-1} + v_{n.t}} where \eqn{h_{n.0}=0}, \eqn{g_n} is an autoregressive parameter and \eqn{v_{n.t}} is a standard normal error term. + +The centred Stochastic Volatility by Chan, Koop, & Yu (2021) +is selected by setting argument \code{centred_sv} of function \code{specify_bsvar_sv$new()} to value \code{TRUE}. +Its conditional variances are given by: +\deqn{Var_{t-1}[u_{n.t}] = exp(h_{n.t})} +where the log-conditional variances \eqn{h_{n.t}} follow an autoregressive process: +\deqn{h_{n.t} = g_n h_{n.t-1} + v_{n.t}} +where \eqn{h_{n.0}=0}, \eqn{g_n} is an autoregressive parameter and \eqn{v_{n.t}} is a zero-mean normal error term +with variance \eqn{s_{v.n}^2}. } \examples{ # simple workflow diff --git a/man/estimate.PosteriorBSVARSV.Rd b/man/estimate.PosteriorBSVARSV.Rd index df302c1..acb9db5 100644 --- a/man/estimate.PosteriorBSVARSV.Rd +++ b/man/estimate.PosteriorBSVARSV.Rd @@ -36,13 +36,17 @@ An object of class PosteriorBSVARSV containing the Bayesian estimation output an \code{last_draw} an object of class BSVARSV with the last draw of the current MCMC run as the starting value to be passed to the continuation of the MCMC estimation using \code{estimate()}. } \description{ -Estimates the SVAR with Stochastic Volatility (SV) heteroskedasticity proposed by Lütkepohl, Shang, Uzeda, and Woźniak (2022). +Estimates the SVAR with Stochastic Volatility (SV) heteroskedasticity +proposed by Lütkepohl, Shang, Uzeda, and Woźniak (2022). Implements the Gibbs sampler proposed by Waggoner & Zha (2003) -for the structural matrix \eqn{B} and the equation-by-equation sampler by Chan, Koop, & Yu (2021) -for the autoregressive slope parameters \eqn{A}. Additionally, the parameter matrices \eqn{A} and \eqn{B} -follow a Minnesota prior and generalised-normal prior distributions respectively with the matrix-specific -overall shrinkage parameters estimated thanks to a 3-level hierarchical prior distribution. The SV model -is estimated in a non-centred parameterisation using a range of techniques including: +for the structural matrix \eqn{B} and the equation-by-equation sampler +by Chan, Koop, & Yu (2021) +for the autoregressive slope parameters \eqn{A}. Additionally, +the parameter matrices \eqn{A} and \eqn{B} +follow a Minnesota prior and generalised-normal prior distributions +respectively with the matrix-specific +overall shrinkage parameters estimated thanks to a 3-level hierarchical prior distribution. +The SV model is estimated using a range of techniques including: simulation smoother, auxiliary mixture, ancillarity-sufficiency interweaving strategy, and generalised inverse Gaussian distribution summarised by Kastner & Frühwirth-Schnatter (2014). See section \bold{Details} for the model equations. @@ -57,14 +61,29 @@ The structural equation is given by \deqn{BE = U} where \eqn{U} is an \code{NxT} matrix of structural form error terms, and \eqn{B} is an \code{NxN} matrix of contemporaneous relationships. - Finally, the structural shocks, \eqn{U}, are temporally and contemporaneously independent and jointly normally distributed with zero mean. -The conditional variance of the \code{n}th shock at time \code{t} is given by: + +Two alternative specifications of the conditional variance of the \code{n}th shock at time \code{t} +can be estimated: non-centred Stochastic Volatility by Lütkepohl, Shang, Uzeda, and Woźniak (2022) +or centred Stochastic Volatility by Chan, Koop, & Yu (2021). + +The non-centred Stochastic Volatility by Lütkepohl, Shang, Uzeda, and Woźniak (2022) +is selected by setting argument \code{centred_sv} of function \code{specify_bsvar_sv$new()} to value \code{FALSE}. +It has the conditional variances given by: \deqn{Var_{t-1}[u_{n.t}] = exp(w_n h_{n.t})} where \eqn{w_n} is the estimated conditional standard deviation of the log-conditional variance and the log-volatility process \eqn{h_{n.t}} follows an autoregressive process: \deqn{h_{n.t} = g_n h_{n.t-1} + v_{n.t}} where \eqn{h_{n.0}=0}, \eqn{g_n} is an autoregressive parameter and \eqn{v_{n.t}} is a standard normal error term. + +The centred Stochastic Volatility by Chan, Koop, & Yu (2021) +is selected by setting argument \code{centred_sv} of function \code{specify_bsvar_sv$new()} to value \code{TRUE}. +Its conditional variances are given by: +\deqn{Var_{t-1}[u_{n.t}] = exp(h_{n.t})} +where the log-conditional variances \eqn{h_{n.t}} follow an autoregressive process: +\deqn{h_{n.t} = g_n h_{n.t-1} + v_{n.t}} +where \eqn{h_{n.0}=0}, \eqn{g_n} is an autoregressive parameter and \eqn{v_{n.t}} is a zero-mean normal error term +with variance \eqn{s_{v.n}^2}. } \examples{ # simple workflow diff --git a/man/specify_bsvar_sv.Rd b/man/specify_bsvar_sv.Rd index 3fcb6bd..10b718e 100644 --- a/man/specify_bsvar_sv.Rd +++ b/man/specify_bsvar_sv.Rd @@ -77,6 +77,8 @@ spec$get_starting_values() \item{\code{data_matrices}}{an object DataMatricesBSVAR with the data matrices.} \item{\code{starting_values}}{an object StartingValuesBSVARSV with the starting values.} + +\item{\code{centred_sv}}{a logical value - if true a centred parameterisation of the Stochastic Volatility process is estimated. Otherwise, its non-centred parameterisation is estimated. See Lütkepohl, Shang, Uzeda, Woźniak (2022) for more info.} } \if{html}{\out{}} } @@ -97,7 +99,13 @@ spec$get_starting_values() \subsection{Method \code{new()}}{ Create a new specification of the BSVAR model with Stochastic Volatility heteroskedasticity, BSVARSV. \subsection{Usage}{ -\if{html}{\out{
}}\preformatted{specify_bsvar_sv$new(data, p = 1L, B, stationary = rep(FALSE, ncol(data)))}\if{html}{\out{
}} +\if{html}{\out{
}}\preformatted{specify_bsvar_sv$new( + data, + p = 1L, + B, + centred_sv = FALSE, + stationary = rep(FALSE, ncol(data)) +)}\if{html}{\out{
}} } \subsection{Arguments}{ @@ -109,6 +117,8 @@ Create a new specification of the BSVAR model with Stochastic Volatility heteros \item{\code{B}}{a logical \code{NxN} matrix containing value \code{TRUE} for the elements of the structural matrix \eqn{B} to be estimated and value \code{FALSE} for exclusion restrictions to be set to zero.} +\item{\code{centred_sv}}{a logical value. If \code{FALSE} a non-centred Stochastic Volatility processes for conditional variances are estimated. Otherwise, a centred process is estimated.} + \item{\code{stationary}}{an \code{N} logical vector - its element set to \code{FALSE} sets the prior mean for the autoregressive parameters of the \code{N}th equation to the white noise process, otherwise to random walk.} } \if{html}{\out{}} diff --git a/man/specify_starting_values_bsvar_sv.Rd b/man/specify_starting_values_bsvar_sv.Rd index dc30623..6b9c216 100644 --- a/man/specify_starting_values_bsvar_sv.Rd +++ b/man/specify_starting_values_bsvar_sv.Rd @@ -51,6 +51,8 @@ sv$set_starting_values(sv_list) # providing to the class object \item{\code{omega}}{an \code{N}-vector with values of SV process conditional standard deviations.} +\item{\code{sigma2v}}{an \code{N}-vector with values of SV process conditional variances.} + \item{\code{S}}{an \code{NxT} integer matrix with the auxiliary mixture component indicators.} \item{\code{sigma2_omega}}{an \code{N}-vector with variances of the zero-mean normal prior for \eqn{\omega_n}.} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 047c582..d3f3468 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -316,8 +316,8 @@ RcppExport SEXP _bsvars_bsvar_msh_cpp(SEXP SSEXP, SEXP YSEXP, SEXP XSEXP, SEXP p return rcpp_result_gen; } // bsvar_sv_cpp -Rcpp::List bsvar_sv_cpp(const int& S, const arma::mat& Y, const arma::mat& X, const Rcpp::List& prior, const arma::field& VB, const Rcpp::List& starting_values, const int thin, const bool sample_s_, const bool show_progress); -static SEXP _bsvars_bsvar_sv_cpp_try(SEXP SSEXP, SEXP YSEXP, SEXP XSEXP, SEXP priorSEXP, SEXP VBSEXP, SEXP starting_valuesSEXP, SEXP thinSEXP, SEXP sample_s_SEXP, SEXP show_progressSEXP) { +Rcpp::List bsvar_sv_cpp(const int& S, const arma::mat& Y, const arma::mat& X, const Rcpp::List& prior, const arma::field& VB, const Rcpp::List& starting_values, const int thin, const bool centred_sv, const bool show_progress); +static SEXP _bsvars_bsvar_sv_cpp_try(SEXP SSEXP, SEXP YSEXP, SEXP XSEXP, SEXP priorSEXP, SEXP VBSEXP, SEXP starting_valuesSEXP, SEXP thinSEXP, SEXP centred_svSEXP, SEXP show_progressSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::traits::input_parameter< const int& >::type S(SSEXP); @@ -327,17 +327,17 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const arma::field& >::type VB(VBSEXP); Rcpp::traits::input_parameter< const Rcpp::List& >::type starting_values(starting_valuesSEXP); Rcpp::traits::input_parameter< const int >::type thin(thinSEXP); - Rcpp::traits::input_parameter< const bool >::type sample_s_(sample_s_SEXP); + Rcpp::traits::input_parameter< const bool >::type centred_sv(centred_svSEXP); Rcpp::traits::input_parameter< const bool >::type show_progress(show_progressSEXP); - rcpp_result_gen = Rcpp::wrap(bsvar_sv_cpp(S, Y, X, prior, VB, starting_values, thin, sample_s_, show_progress)); + rcpp_result_gen = Rcpp::wrap(bsvar_sv_cpp(S, Y, X, prior, VB, starting_values, thin, centred_sv, show_progress)); return rcpp_result_gen; END_RCPP_RETURN_ERROR } -RcppExport SEXP _bsvars_bsvar_sv_cpp(SEXP SSEXP, SEXP YSEXP, SEXP XSEXP, SEXP priorSEXP, SEXP VBSEXP, SEXP starting_valuesSEXP, SEXP thinSEXP, SEXP sample_s_SEXP, SEXP show_progressSEXP) { +RcppExport SEXP _bsvars_bsvar_sv_cpp(SEXP SSEXP, SEXP YSEXP, SEXP XSEXP, SEXP priorSEXP, SEXP VBSEXP, SEXP starting_valuesSEXP, SEXP thinSEXP, SEXP centred_svSEXP, SEXP show_progressSEXP) { SEXP rcpp_result_gen; { Rcpp::RNGScope rcpp_rngScope_gen; - rcpp_result_gen = PROTECT(_bsvars_bsvar_sv_cpp_try(SSEXP, YSEXP, XSEXP, priorSEXP, VBSEXP, starting_valuesSEXP, thinSEXP, sample_s_SEXP, show_progressSEXP)); + rcpp_result_gen = PROTECT(_bsvars_bsvar_sv_cpp_try(SSEXP, YSEXP, XSEXP, priorSEXP, VBSEXP, starting_valuesSEXP, thinSEXP, centred_svSEXP, show_progressSEXP)); } Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error"); if (rcpp_isInterrupt_gen) { @@ -435,8 +435,8 @@ RcppExport SEXP _bsvars_forecast_bsvar_msh(SEXP posterior_BSEXP, SEXP posterior_ return rcpp_result_gen; } // forecast_bsvar_sv -Rcpp::List forecast_bsvar_sv(arma::cube& posterior_B, arma::cube& posterior_A, arma::vec& posterior_h_T, arma::mat& posterior_rho, arma::mat& posterior_omega, arma::vec& X_T, const int horizon); -static SEXP _bsvars_forecast_bsvar_sv_try(SEXP posterior_BSEXP, SEXP posterior_ASEXP, SEXP posterior_h_TSEXP, SEXP posterior_rhoSEXP, SEXP posterior_omegaSEXP, SEXP X_TSEXP, SEXP horizonSEXP) { +Rcpp::List forecast_bsvar_sv(arma::cube& posterior_B, arma::cube& posterior_A, arma::vec& posterior_h_T, arma::mat& posterior_rho, arma::mat& posterior_omega, arma::vec& X_T, const int horizon, const bool centred_sv); +static SEXP _bsvars_forecast_bsvar_sv_try(SEXP posterior_BSEXP, SEXP posterior_ASEXP, SEXP posterior_h_TSEXP, SEXP posterior_rhoSEXP, SEXP posterior_omegaSEXP, SEXP X_TSEXP, SEXP horizonSEXP, SEXP centred_svSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::traits::input_parameter< arma::cube& >::type posterior_B(posterior_BSEXP); @@ -446,15 +446,16 @@ BEGIN_RCPP Rcpp::traits::input_parameter< arma::mat& >::type posterior_omega(posterior_omegaSEXP); Rcpp::traits::input_parameter< arma::vec& >::type X_T(X_TSEXP); Rcpp::traits::input_parameter< const int >::type horizon(horizonSEXP); - rcpp_result_gen = Rcpp::wrap(forecast_bsvar_sv(posterior_B, posterior_A, posterior_h_T, posterior_rho, posterior_omega, X_T, horizon)); + Rcpp::traits::input_parameter< const bool >::type centred_sv(centred_svSEXP); + rcpp_result_gen = Rcpp::wrap(forecast_bsvar_sv(posterior_B, posterior_A, posterior_h_T, posterior_rho, posterior_omega, X_T, horizon, centred_sv)); return rcpp_result_gen; END_RCPP_RETURN_ERROR } -RcppExport SEXP _bsvars_forecast_bsvar_sv(SEXP posterior_BSEXP, SEXP posterior_ASEXP, SEXP posterior_h_TSEXP, SEXP posterior_rhoSEXP, SEXP posterior_omegaSEXP, SEXP X_TSEXP, SEXP horizonSEXP) { +RcppExport SEXP _bsvars_forecast_bsvar_sv(SEXP posterior_BSEXP, SEXP posterior_ASEXP, SEXP posterior_h_TSEXP, SEXP posterior_rhoSEXP, SEXP posterior_omegaSEXP, SEXP X_TSEXP, SEXP horizonSEXP, SEXP centred_svSEXP) { SEXP rcpp_result_gen; { Rcpp::RNGScope rcpp_rngScope_gen; - rcpp_result_gen = PROTECT(_bsvars_forecast_bsvar_sv_try(posterior_BSEXP, posterior_ASEXP, posterior_h_TSEXP, posterior_rhoSEXP, posterior_omegaSEXP, X_TSEXP, horizonSEXP)); + rcpp_result_gen = PROTECT(_bsvars_forecast_bsvar_sv_try(posterior_BSEXP, posterior_ASEXP, posterior_h_TSEXP, posterior_rhoSEXP, posterior_omegaSEXP, X_TSEXP, horizonSEXP, centred_svSEXP)); } Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error"); if (rcpp_isInterrupt_gen) { @@ -1384,28 +1385,72 @@ RcppExport SEXP _bsvars_find_mixture_indicator_cdf(SEXP datanormSEXP) { return rcpp_result_gen; } // svar_nc1 -Rcpp::List svar_nc1(arma::rowvec& aux_h_n, double& aux_rho_n, double& aux_omega_n, double& aux_sigma2_omega_n, double& aux_s_n, arma::urowvec& aux_S_n, const arma::rowvec& u, const Rcpp::List& prior, bool sample_s_); -static SEXP _bsvars_svar_nc1_try(SEXP aux_h_nSEXP, SEXP aux_rho_nSEXP, SEXP aux_omega_nSEXP, SEXP aux_sigma2_omega_nSEXP, SEXP aux_s_nSEXP, SEXP aux_S_nSEXP, SEXP uSEXP, SEXP priorSEXP, SEXP sample_s_SEXP) { +Rcpp::List svar_nc1(arma::rowvec& aux_h_n, double& aux_rho_n, double& aux_omega_n, double& aux_sigma2v_n, double& aux_sigma2_omega_n, double& aux_s_n, arma::urowvec& aux_S_n, const arma::rowvec& u, const Rcpp::List& prior, bool sample_s_); +static SEXP _bsvars_svar_nc1_try(SEXP aux_h_nSEXP, SEXP aux_rho_nSEXP, SEXP aux_omega_nSEXP, SEXP aux_sigma2v_nSEXP, SEXP aux_sigma2_omega_nSEXP, SEXP aux_s_nSEXP, SEXP aux_S_nSEXP, SEXP uSEXP, SEXP priorSEXP, SEXP sample_s_SEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::traits::input_parameter< arma::rowvec& >::type aux_h_n(aux_h_nSEXP); + Rcpp::traits::input_parameter< double& >::type aux_rho_n(aux_rho_nSEXP); + Rcpp::traits::input_parameter< double& >::type aux_omega_n(aux_omega_nSEXP); + Rcpp::traits::input_parameter< double& >::type aux_sigma2v_n(aux_sigma2v_nSEXP); + Rcpp::traits::input_parameter< double& >::type aux_sigma2_omega_n(aux_sigma2_omega_nSEXP); + Rcpp::traits::input_parameter< double& >::type aux_s_n(aux_s_nSEXP); + Rcpp::traits::input_parameter< arma::urowvec& >::type aux_S_n(aux_S_nSEXP); + Rcpp::traits::input_parameter< const arma::rowvec& >::type u(uSEXP); + Rcpp::traits::input_parameter< const Rcpp::List& >::type prior(priorSEXP); + Rcpp::traits::input_parameter< bool >::type sample_s_(sample_s_SEXP); + rcpp_result_gen = Rcpp::wrap(svar_nc1(aux_h_n, aux_rho_n, aux_omega_n, aux_sigma2v_n, aux_sigma2_omega_n, aux_s_n, aux_S_n, u, prior, sample_s_)); + return rcpp_result_gen; +END_RCPP_RETURN_ERROR +} +RcppExport SEXP _bsvars_svar_nc1(SEXP aux_h_nSEXP, SEXP aux_rho_nSEXP, SEXP aux_omega_nSEXP, SEXP aux_sigma2v_nSEXP, SEXP aux_sigma2_omega_nSEXP, SEXP aux_s_nSEXP, SEXP aux_S_nSEXP, SEXP uSEXP, SEXP priorSEXP, SEXP sample_s_SEXP) { + SEXP rcpp_result_gen; + { + Rcpp::RNGScope rcpp_rngScope_gen; + rcpp_result_gen = PROTECT(_bsvars_svar_nc1_try(aux_h_nSEXP, aux_rho_nSEXP, aux_omega_nSEXP, aux_sigma2v_nSEXP, aux_sigma2_omega_nSEXP, aux_s_nSEXP, aux_S_nSEXP, uSEXP, priorSEXP, sample_s_SEXP)); + } + Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error"); + if (rcpp_isInterrupt_gen) { + UNPROTECT(1); + Rf_onintr(); + } + bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen); + if (rcpp_isLongjump_gen) { + Rcpp::internal::resumeJump(rcpp_result_gen); + } + Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error"); + if (rcpp_isError_gen) { + SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen); + UNPROTECT(1); + Rf_error(CHAR(rcpp_msgSEXP_gen)); + } + UNPROTECT(1); + return rcpp_result_gen; +} +// svar_ce1 +Rcpp::List svar_ce1(arma::rowvec& aux_h_n, double& aux_rho_n, double& aux_omega_n, double& aux_sigma2v_n, double& aux_sigma2_omega_n, double& aux_s_n, arma::urowvec& aux_S_n, const arma::rowvec& u, const Rcpp::List& prior, bool sample_s_); +static SEXP _bsvars_svar_ce1_try(SEXP aux_h_nSEXP, SEXP aux_rho_nSEXP, SEXP aux_omega_nSEXP, SEXP aux_sigma2v_nSEXP, SEXP aux_sigma2_omega_nSEXP, SEXP aux_s_nSEXP, SEXP aux_S_nSEXP, SEXP uSEXP, SEXP priorSEXP, SEXP sample_s_SEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::traits::input_parameter< arma::rowvec& >::type aux_h_n(aux_h_nSEXP); Rcpp::traits::input_parameter< double& >::type aux_rho_n(aux_rho_nSEXP); Rcpp::traits::input_parameter< double& >::type aux_omega_n(aux_omega_nSEXP); + Rcpp::traits::input_parameter< double& >::type aux_sigma2v_n(aux_sigma2v_nSEXP); Rcpp::traits::input_parameter< double& >::type aux_sigma2_omega_n(aux_sigma2_omega_nSEXP); Rcpp::traits::input_parameter< double& >::type aux_s_n(aux_s_nSEXP); Rcpp::traits::input_parameter< arma::urowvec& >::type aux_S_n(aux_S_nSEXP); Rcpp::traits::input_parameter< const arma::rowvec& >::type u(uSEXP); Rcpp::traits::input_parameter< const Rcpp::List& >::type prior(priorSEXP); Rcpp::traits::input_parameter< bool >::type sample_s_(sample_s_SEXP); - rcpp_result_gen = Rcpp::wrap(svar_nc1(aux_h_n, aux_rho_n, aux_omega_n, aux_sigma2_omega_n, aux_s_n, aux_S_n, u, prior, sample_s_)); + rcpp_result_gen = Rcpp::wrap(svar_ce1(aux_h_n, aux_rho_n, aux_omega_n, aux_sigma2v_n, aux_sigma2_omega_n, aux_s_n, aux_S_n, u, prior, sample_s_)); return rcpp_result_gen; END_RCPP_RETURN_ERROR } -RcppExport SEXP _bsvars_svar_nc1(SEXP aux_h_nSEXP, SEXP aux_rho_nSEXP, SEXP aux_omega_nSEXP, SEXP aux_sigma2_omega_nSEXP, SEXP aux_s_nSEXP, SEXP aux_S_nSEXP, SEXP uSEXP, SEXP priorSEXP, SEXP sample_s_SEXP) { +RcppExport SEXP _bsvars_svar_ce1(SEXP aux_h_nSEXP, SEXP aux_rho_nSEXP, SEXP aux_omega_nSEXP, SEXP aux_sigma2v_nSEXP, SEXP aux_sigma2_omega_nSEXP, SEXP aux_s_nSEXP, SEXP aux_S_nSEXP, SEXP uSEXP, SEXP priorSEXP, SEXP sample_s_SEXP) { SEXP rcpp_result_gen; { Rcpp::RNGScope rcpp_rngScope_gen; - rcpp_result_gen = PROTECT(_bsvars_svar_nc1_try(aux_h_nSEXP, aux_rho_nSEXP, aux_omega_nSEXP, aux_sigma2_omega_nSEXP, aux_s_nSEXP, aux_S_nSEXP, uSEXP, priorSEXP, sample_s_SEXP)); + rcpp_result_gen = PROTECT(_bsvars_svar_ce1_try(aux_h_nSEXP, aux_rho_nSEXP, aux_omega_nSEXP, aux_sigma2v_nSEXP, aux_sigma2_omega_nSEXP, aux_s_nSEXP, aux_S_nSEXP, uSEXP, priorSEXP, sample_s_SEXP)); } Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error"); if (rcpp_isInterrupt_gen) { @@ -1509,7 +1554,7 @@ static int _bsvars_RcppExport_validate(const char* sig) { signatures.insert("Rcpp::List(*bsvar_sv_cpp)(const int&,const arma::mat&,const arma::mat&,const Rcpp::List&,const arma::field&,const Rcpp::List&,const int,const bool,const bool)"); signatures.insert("Rcpp::List(*forecast_bsvar)(arma::cube&,arma::cube&,arma::vec&,const int)"); signatures.insert("Rcpp::List(*forecast_bsvar_msh)(arma::cube&,arma::cube&,arma::cube&,arma::cube&,arma::vec&,arma::vec&,const int)"); - signatures.insert("Rcpp::List(*forecast_bsvar_sv)(arma::cube&,arma::cube&,arma::vec&,arma::mat&,arma::mat&,arma::vec&,const int)"); + signatures.insert("Rcpp::List(*forecast_bsvar_sv)(arma::cube&,arma::cube&,arma::vec&,arma::mat&,arma::mat&,arma::vec&,const int,const bool)"); signatures.insert("arma::vec(*Ergodic_PR_TR)(const arma::mat&)"); signatures.insert("arma::mat(*count_regime_transitions)(const arma::mat&)"); signatures.insert("arma::rowvec(*rDirichlet1)(const arma::rowvec&)"); @@ -1535,7 +1580,8 @@ static int _bsvars_RcppExport_validate(const char* sig) { signatures.insert("arma::vec(*precision_sampler_ar1)(const arma::vec&,const double&,const arma::vec&)"); signatures.insert("arma::uvec(*inverse_transform_sampling)(const arma::vec&,const int)"); signatures.insert("arma::vec(*find_mixture_indicator_cdf)(const arma::vec&)"); - signatures.insert("Rcpp::List(*svar_nc1)(arma::rowvec&,double&,double&,double&,double&,arma::urowvec&,const arma::rowvec&,const Rcpp::List&,bool)"); + signatures.insert("Rcpp::List(*svar_nc1)(arma::rowvec&,double&,double&,double&,double&,double&,arma::urowvec&,const arma::rowvec&,const Rcpp::List&,bool)"); + signatures.insert("Rcpp::List(*svar_ce1)(arma::rowvec&,double&,double&,double&,double&,double&,arma::urowvec&,const arma::rowvec&,const Rcpp::List&,bool)"); signatures.insert("arma::mat(*orthogonal_complement_matrix_TW)(const arma::mat&)"); signatures.insert("arma::vec(*log_mean)(arma::mat)"); } @@ -1582,6 +1628,7 @@ RcppExport SEXP _bsvars_RcppExport_registerCCallable() { R_RegisterCCallable("bsvars", "_bsvars_inverse_transform_sampling", (DL_FUNC)_bsvars_inverse_transform_sampling_try); R_RegisterCCallable("bsvars", "_bsvars_find_mixture_indicator_cdf", (DL_FUNC)_bsvars_find_mixture_indicator_cdf_try); R_RegisterCCallable("bsvars", "_bsvars_svar_nc1", (DL_FUNC)_bsvars_svar_nc1_try); + R_RegisterCCallable("bsvars", "_bsvars_svar_ce1", (DL_FUNC)_bsvars_svar_ce1_try); R_RegisterCCallable("bsvars", "_bsvars_orthogonal_complement_matrix_TW", (DL_FUNC)_bsvars_orthogonal_complement_matrix_TW_try); R_RegisterCCallable("bsvars", "_bsvars_log_mean", (DL_FUNC)_bsvars_log_mean_try); R_RegisterCCallable("bsvars", "_bsvars_RcppExport_validate", (DL_FUNC)_bsvars_RcppExport_validate); @@ -1600,7 +1647,7 @@ static const R_CallMethodDef CallEntries[] = { {"_bsvars_bsvar_sv_cpp", (DL_FUNC) &_bsvars_bsvar_sv_cpp, 9}, {"_bsvars_forecast_bsvar", (DL_FUNC) &_bsvars_forecast_bsvar, 4}, {"_bsvars_forecast_bsvar_msh", (DL_FUNC) &_bsvars_forecast_bsvar_msh, 7}, - {"_bsvars_forecast_bsvar_sv", (DL_FUNC) &_bsvars_forecast_bsvar_sv, 7}, + {"_bsvars_forecast_bsvar_sv", (DL_FUNC) &_bsvars_forecast_bsvar_sv, 8}, {"_bsvars_Ergodic_PR_TR", (DL_FUNC) &_bsvars_Ergodic_PR_TR, 1}, {"_bsvars_count_regime_transitions", (DL_FUNC) &_bsvars_count_regime_transitions, 1}, {"_bsvars_rDirichlet1", (DL_FUNC) &_bsvars_rDirichlet1, 1}, @@ -1626,7 +1673,8 @@ static const R_CallMethodDef CallEntries[] = { {"_bsvars_precision_sampler_ar1", (DL_FUNC) &_bsvars_precision_sampler_ar1, 3}, {"_bsvars_inverse_transform_sampling", (DL_FUNC) &_bsvars_inverse_transform_sampling, 2}, {"_bsvars_find_mixture_indicator_cdf", (DL_FUNC) &_bsvars_find_mixture_indicator_cdf, 1}, - {"_bsvars_svar_nc1", (DL_FUNC) &_bsvars_svar_nc1, 9}, + {"_bsvars_svar_nc1", (DL_FUNC) &_bsvars_svar_nc1, 10}, + {"_bsvars_svar_ce1", (DL_FUNC) &_bsvars_svar_ce1, 10}, {"_bsvars_orthogonal_complement_matrix_TW", (DL_FUNC) &_bsvars_orthogonal_complement_matrix_TW, 1}, {"_bsvars_log_mean", (DL_FUNC) &_bsvars_log_mean, 1}, {"_bsvars_RcppExport_registerCCallable", (DL_FUNC) &_bsvars_RcppExport_registerCCallable, 0}, diff --git a/src/bsvar_sv.cpp b/src/bsvar_sv.cpp index abbc7c6..aaddaca 100644 --- a/src/bsvar_sv.cpp +++ b/src/bsvar_sv.cpp @@ -21,16 +21,25 @@ Rcpp::List bsvar_sv_cpp ( const arma::field& VB, // restrictions on B0 const Rcpp::List& starting_values, const int thin = 100, // introduce thinning - const bool sample_s_ = true, + const bool centred_sv = false, const bool show_progress = true ) { // Progress bar setup vec prog_rep_points = arma::round(arma::linspace(0, S, 50)); + + std::string name_model = ""; + if ( centred_sv ) { + name_model = " Centred"; + } else { + name_model = "Non-centred"; + } + if (show_progress) { Rcout << "**************************************************|" << endl; Rcout << "bsvars: Bayesian Structural Vector Autoregressions|" << endl; Rcout << "**************************************************|" << endl; Rcout << " Gibbs sampler for the SVAR-SV model |" << endl; + Rcout << " " << name_model << " SV model is estimated |" << endl; Rcout << "**************************************************|" << endl; Rcout << " Progress of the MCMC simulation for " << S << " draws" << endl; Rcout << " Every " << thin << "th draw is saved via MCMC thinning" << endl; @@ -49,13 +58,20 @@ Rcpp::List bsvar_sv_cpp ( mat aux_h = as(starting_values["h"]); vec aux_rho = as(starting_values["rho"]); vec aux_omega = as(starting_values["omega"]); + vec aux_sigma2v = as(starting_values["sigma2v"]); umat aux_S = as(starting_values["S"]); vec aux_sigma2_omega = as(starting_values["sigma2_omega"]); vec aux_s_ = as(starting_values["s_"]); mat aux_sigma(N, T); - for (int n=0; n(sv_n["aux_h_n"]); aux_rho(n) = as(sv_n["aux_rho_n"]); aux_omega(n) = as(sv_n["aux_omega_n"]); + aux_sigma2v(n) = as(sv_n["aux_sigma2v_n"]); aux_S.row(n) = as(sv_n["aux_S_n"]); aux_sigma2_omega(n) = as(sv_n["aux_sigma2_omega_n"]); aux_s_(n) = as(sv_n["aux_s_n"]); - aux_sigma.row(n) = exp(0.5 * aux_omega(n) * aux_h.row(n)); + if ( centred_sv ) { + aux_sigma.row(n) = exp(0.5 * aux_h.row(n)); + } else { + aux_sigma.row(n) = exp(0.5 * aux_omega(n) * aux_h.row(n)); + } } if (s % thin == 0) { @@ -120,6 +148,7 @@ Rcpp::List bsvar_sv_cpp ( posterior_h.slice(ss) = aux_h; posterior_rho.col(ss) = aux_rho; posterior_omega.col(ss) = aux_omega; + posterior_sigma2v.col(ss) = aux_sigma2v; posterior_S.slice(ss) = aux_S; posterior_sigma2_omega.col(ss) = aux_sigma2_omega; posterior_s_.col(ss) = aux_s_; @@ -136,6 +165,7 @@ Rcpp::List bsvar_sv_cpp ( _["h"] = aux_h, _["rho"] = aux_rho, _["omega"] = aux_omega, + _["sigma2v"] = aux_sigma2v, _["S"] = aux_S, _["sigma2_omega"] = aux_sigma2_omega, _["s_"] = aux_s_, @@ -148,6 +178,7 @@ Rcpp::List bsvar_sv_cpp ( _["h"] = posterior_h, _["rho"] = posterior_rho, _["omega"] = posterior_omega, + _["sigma2v"] = posterior_sigma2v, _["S"] = posterior_S, _["sigma2_omega"] = posterior_sigma2_omega, _["s_"] = posterior_s_, diff --git a/src/bsvar_sv.h b/src/bsvar_sv.h index 7c86a8b..fbbc899 100644 --- a/src/bsvar_sv.h +++ b/src/bsvar_sv.h @@ -12,7 +12,7 @@ Rcpp::List bsvar_sv_cpp ( const arma::field& VB, // restrictions on B0 const Rcpp::List& starting_values, const int thin = 100, // introduce thinning - const bool sample_s_ = true, + const bool centred_sv = false, const bool show_progress = true ); diff --git a/src/forecast.cpp b/src/forecast.cpp index a824612..cafc2c6 100644 --- a/src/forecast.cpp +++ b/src/forecast.cpp @@ -108,7 +108,8 @@ Rcpp::List forecast_bsvar_sv ( arma::mat& posterior_rho, // NxS arma::mat& posterior_omega, // NxS arma::vec& X_T, // (K) - const int horizon + const int horizon, + const bool centred_sv = FALSE ) { const int N = posterior_B.n_rows; @@ -133,7 +134,11 @@ Rcpp::List forecast_bsvar_sv ( for (int n=0; n