diff --git a/DESCRIPTION b/DESCRIPTION index 890b60b..5a54312 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -5,7 +5,7 @@ Date: 2023-06-30 Authors@R: person("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 and inference of Bayesian global vector autoregressive (BGVAR) and error correction (BGVEC) models. License: GPL (>= 2) -Depends: R (>= 3.3.0), bvartools (>= 0.2.1) +Depends: R (>= 3.3.0), bvartools (>= 0.2.3) Imports: coda, grDevices, @@ -24,8 +24,10 @@ Suggests: knitr, rmarkdown Collate: 'RcppExports.R' 'add_priors.gvarsubmodels.R' + 'add_priors.gvecsubmodels.R' 'bgvars-package.R' 'bvarxpost.R' + 'bvecxpost.R' 'combine_submodels.R' 'country_fill_helper.R' 'create_models.R' @@ -37,6 +39,7 @@ Collate: 'data.R' 'diff_variables.R' 'draw_posterior.gvarsubmodels.R' + 'draw_posterior.gvecsubmodels.R' 'gen_varx.R' 'gen_vecx.R' 'get_regressor_names.R' @@ -50,16 +53,23 @@ Collate: 'plot.bgvarfevd.R' 'plot.bgvarirf.R' 'plot.bgvarprd.R' + 'plot.bgvecest.R' 'plot.ctryvarest.R' + 'plot.ctryvecest.R' 'plot.teststats.bgvarest.R' 'predict.bgvar.R' 'summary.bgvarest.R' 'print.summary.bgvarest.R' + 'summary.bgvecest.R' + 'print.summary.bgvecest.R' 'summary.ctryvarest.R' 'print.summary.ctryvarest.R' + 'summary.ctryvecest.R' + 'print.summary.ctryvecest.R' 'select_submodels.R' 'submodel_test_statistics.R' 'thin.bgvar.R' 'thin.bgvarest.R' + 'thin.bgvecest.R' 'tvpribbon.R' 'zzz.R' diff --git a/NAMESPACE b/NAMESPACE index 2d848f9..f84ee20 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,7 +1,9 @@ # Generated by roxygen2: do not edit by hand S3method(add_priors,gvarsubmodels) +S3method(add_priors,gvecsubmodels) S3method(draw_posterior,gvarsubmodels) +S3method(draw_posterior,gvecsubmodels) S3method(irf,bgvarest) S3method(irf,ctryvarest) S3method(plot,bgvarest) @@ -9,16 +11,24 @@ S3method(plot,bgvarestirf) S3method(plot,bgvarfevd) S3method(plot,bgvarirf) S3method(plot,bgvarprd) +S3method(plot,bgvecest) S3method(plot,ctryvarest) +S3method(plot,ctryvecest) S3method(plot,teststats.bgvarest) S3method(predict,bgvar) S3method(print,summary.bgvarest) +S3method(print,summary.bgvecest) S3method(print,summary.ctryvarest) +S3method(print,summary.ctryvecest) S3method(summary,bgvarest) +S3method(summary,bgvecest) S3method(summary,ctryvarest) +S3method(summary,ctryvecest) S3method(thin,bgvar) S3method(thin,bgvarest) +S3method(thin,bgvecest) export(bvarxpost) +export(bvecxpost) export(combine_submodels) export(create_models) export(create_regions) diff --git a/R/RcppExports.R b/R/RcppExports.R index cda8ebb..e4e5716 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -9,6 +9,14 @@ .Call(`_bgvars_bgvartvpalg`, object) } +.bgvecalg <- function(object) { + .Call(`_bgvars_bgvecalg`, object) +} + +.bgvectvpalg <- function(object) { + .Call(`_bgvars_bgvectvpalg`, object) +} + .draw_forecast <- function(i, k, a0, a, b_, c_, sigma, pred) { .Call(`_bgvars_draw_forecast`, i, k, a0, a, b_, c_, sigma, pred) } diff --git a/R/add_priors.gvecsubmodels.R b/R/add_priors.gvecsubmodels.R new file mode 100644 index 0000000..62c3d3b --- /dev/null +++ b/R/add_priors.gvecsubmodels.R @@ -0,0 +1,867 @@ +#' Add Priors to Country-Specific VECX Models of a GVAR Model +#' +#' Adds prior specifications to a list of country models, which was produced by +#' the function \code{\link{create_models}}. +#' +#' @param object a named list, usually, the output of a call to \code{\link{create_models}}. +#' @param coef a named list of prior specifications for the coefficients of the country-specific +#' models. For the default specification all prior means are set to zero and the diagonal elements of +#' the inverse prior variance-covariance matrix are set to 1 for coefficients corresponding to non-deterministic +#' terms. For deterministic coefficients the prior variances are set to 10 via \code{v_i_det = 0.1}. +#' The variances need to be specified as precisions, i.e. as inverses of the variances. +#' For further specifications see 'Details'. +#' @param coint a named list of prior specifications for cointegration coefficients of +#' country-specific VEC models. See 'Details'. +#' @param sigma a named list of prior specifications for the error variance-covariance matrix +#' of the country models. For the default specification of an inverse Wishart distribution +#' the prior degrees of freedom are set to the number of endogenous variables and +#' the prior variances to 1. See 'Details'. +#' @param ssvs a named list of prior specifications for the SSVS algorithm. See 'Details'. +#' @param bvs a named list of prior specifications for the BVS algorithm. See 'Details'. +#' @param ... further arguments passed to or from other methods. +#' +#' @details The argument \code{coef} can contain the following elements +#' \describe{ +#' \item{\code{v_i}}{a numeric specifying the prior precision of the coefficients. Default is 1.} +#' \item{\code{v_i_det}}{a numeric specifying the prior precision of coefficients corresponding to deterministic terms. Default is 0.1.} +#' \item{\code{coint_var}}{a logical specifying whether the prior mean of the first own lag of an +#' endogenous variable in a VAR model should be set to 1. Default is \code{FALSE}.} +#' \item{\code{const}}{a numeric or character specifying the prior mean of coefficients, which correspond +#' to the intercept. If a numeric is provided, all prior means are set to this value. +#' If \code{const = "mean"}, the means of the series of endogenous variables are used as prior means. +#' If \code{const = "first"}, the first values of the series of endogenous variables are used as prior means.} +#' \item{\code{minnesota}}{a list of length 4 containing parameters for the calculation of +#' the Minnesota prior, where the element names are \code{kappa0}, \code{kappa1}, \code{kappa2} and \code{kappa3}. +#' For the endogenous variable \eqn{i} the prior variance of the \eqn{l}th +#' lag of regressor \eqn{j} is obtained as +#' \deqn{ \frac{\kappa_{0}}{l^2} \textrm{ for own lags of endogenous variables,}} +#' \deqn{ \frac{\kappa_{0} \kappa_{1}}{l^2} \frac{\sigma_{i}^2}{\sigma_{j}^2} \textrm{ for endogenous variables other than own lags,}} +#' \deqn{ \frac{\kappa_{0} \kappa_{2}}{(1 + l)^2} \frac{\sigma_{i}^2}{\sigma_{j}^2} \textrm{ for foreign and global exogenous variables,}} +#' \deqn{ \kappa_{0} \kappa_{3} \sigma_{i}^2 \textrm{ for deterministic terms,}} +#' where \eqn{\sigma_{i}} is the residual standard deviation of variable \eqn{i} of an unrestricted +#' LS estimate. For exogenous variables \eqn{\sigma_{i}} is the sample standard deviation. +#' If \code{kappa2 = NULL}, \eqn{\kappa_{0} \kappa_{3} \sigma_{i}^2} will be used for foreign +#' and global exogenous variables instead. +#' +#' For VEC models the function only provides priors for the non-cointegration part of the model. The +#' residual standard errors \eqn{\sigma_i} are based on an unrestricted LS regression of the +#' endogenous variables on the error correction term and the non-cointegration regressors.} +#' \item{\code{max_var}}{a numeric specifying the maximum prior variance that is allowed for +#' non-deterministic coefficients.} +#' \item{\code{shape}}{an integer specifying the prior degrees of freedom of the error term of the state equation. Default is 3.} +#' \item{\code{rate}}{a numeric specifying the prior error variance of the state equation. Default is 0.0001.} +#' \item{\code{rate_det}}{a numeric specifying the prior error variance of the state equation corresponding to deterministic terms. Default is 0.0001.} +#' } +#' If \code{minnesota} is specified, \code{v_i} and \code{v_i_det} are ignored. +#' +#' The argument \code{coint} can contain the following elements +#' \describe{ +#' \item{\code{coint_v_i}}{numeric between 0 and 1 specifying the shrinkage of the cointegration space prior. Default is 0.} +#' \item{\code{coint_p_tau_i}}{numeric of the diagonal elements of the inverse prior matrix of +#' the central location of the cointegration space \eqn{sp(\beta)}. Default is 1.} +#' } +#' +#' Argument \code{sigma} can contain the following elements: +#' \describe{ +#' \item{\code{df}}{an integer or character specifying the prior degrees of freedom of the error term. Only +#' used, if the prior is inverse Wishart. +#' Default is \code{"k"}, which indicates the amount of endogenous variables in the respective model. +#' \code{"k + 3"} can be used to set the prior to the amount of endogenous variables plus 3. If an integer +#' is provided, the degrees of freedom are set to this value in all models. +#' In all cases the rank \eqn{r} of the cointegration matrix is automatically added.} +#' \item{\code{scale}}{a numeric specifying the prior error variance of endogenous variables. +#' Default is 1.} +#' \item{\code{shape}}{a numeric or character specifying the prior shape parameter of the error term. Only +#' used, if the prior is inverse gamma or if time varying volatilities are estimated. +#' For models with constant volatility the default is \code{"k"}, which indicates the amount of endogenous +#' variables in the respective country model. \code{"k + 3"} can be used to set the prior to the amount of +#' endogenous variables plus 3. If a numeric is provided, the shape parameters are set to this value in all +#' models. For models with stochastic volatility this prior refers to the error variance of the state +#' equation.} +#' \item{\code{rate}}{a numeric specifying the prior rate parameter of the error term. Only used, if the +#' prior is inverse gamma or if time varying volatilities are estimated. For models with stochastic +#' volatility this prior refers to the error variance of the state equation.} +#' \item{\code{mu}}{numeric of the prior mean of the initial state of the log-volatilities. +#' Only used for models with time varying volatility.} +#' \item{\code{v_i}}{numeric of the prior precision of the initial state of the log-volatilities. +#' Only used for models with time varying volatility.} +#' \item{\code{sigma_h}}{numeric of the initial draw for the variance of the log-volatilities. +#' Only used for models with time varying volatility.} +#' \item{\code{constant}}{numeric of the constant, which is added before taking the log of the squared errors. +#' Only used for models with time varying volatility.} +#' \item{\code{covar}}{logical indicating whether error covariances should be estimated. Only used +#' in combination with an inverse gamma prior or stochastic volatility, for which \code{shape} and +#' \code{rate} must be specified.} +#' } +#' \code{df} and \code{scale} must be specified for an inverse Wishart prior. \code{shape} and \code{rate} +#' are required for an inverse gamma prior. For structural models or models with stochastic volatility +#' only a gamma prior specification is allowed. +#' +#' Argument \code{ssvs} can contain the following elements: +#' \describe{ +#' \item{\code{inprior}}{a numeric between 0 and 1 specifying the prior probability +#' of a variable to be included in the model.} +#' \item{\code{tau}}{a numeric vector of two elements containing the prior standard errors +#' of restricted variables (\eqn{\tau_0}) as its first element and unrestricted variables (\eqn{\tau_1}) +#' as its second.} +#' \item{\code{semiautomatic}}{an numeric vector of two elements containing the +#' factors by which the standard errors associated with an unconstrained least squares +#' estimate of the model are multiplied to obtain the prior standard errors +#' of restricted (\eqn{\tau_0}) and unrestricted (\eqn{\tau_1}) variables, respectively. +#' This is the semiautomatic approach described in George et al. (2008).} +#' \item{\code{covar}}{logical indicating if SSVS should also be applied to the error covariance matrix +#' as in George et al. (2008).} +#' \item{\code{exclude_det}}{logical indicating if deterministic terms should be excluded from the SSVS algorithm.} +#' \item{\code{minnesota}}{a numeric vector of length 4 containing parameters for the calculation of +#' the Minnesota-like inclusion priors. See below.} +#' } +#' Either \code{tau} or \code{semiautomatic} must be specified. +#' +#' The argument \code{bvs} can contain the following elements +#' \describe{ +#' \item{\code{inprior}}{a numeric between 0 and 1 specifying the prior probability +#' of a variable to be included in the model.} +#' \item{\code{covar}}{logical indicating if BVS should also be applied to the error covariance matrix.} +#' \item{\code{exclude_det}}{logical indicating if deterministic terms should be excluded from the BVS algorithm.} +#' \item{\code{minnesota}}{a numeric vector of length 4 containing parameters for the calculation of +#' the Minnesota-like inclusion priors. See below.} +#' } +#' +#' If either \code{ssvs$minnesota} or \code{bvs$minnesota} is specified, prior inclusion probabilities +#' are calculated in a Minnesota-like fashion as +#' \tabular{cl}{ +#' \eqn{\frac{\kappa_1}{l}} \tab for own lags of endogenous variables, \cr +#' \eqn{\frac{\kappa_2}{l}} \tab for other endogenous variables, \cr +#' \eqn{\frac{\kappa_3}{1 + l}} \tab for exogenous variables, \cr +#' \eqn{\kappa_{4}} \tab for deterministic variables, +#' } +#' for lag \eqn{l} with \eqn{\kappa_1}, \eqn{\kappa_2}, \eqn{\kappa_3}, \eqn{\kappa_4} as the first, second, +#' third and forth element in \code{ssvs$minnesota} or \code{bvs$minnesota}, respectively. +#' +#' @return A list of country models. +#' +#' @references +#' +#' Chan, J., Koop, G., Poirier, D. J., & Tobias J. L. (2019). \emph{Bayesian econometric methods} +#' (2nd ed.). Cambridge: Cambridge University Press. +#' +#' George, E. I., Sun, D., & Ni, S. (2008). Bayesian stochastic search for VAR model +#' restrictions. \emph{Journal of Econometrics, 142}(1), 553--580. +#' \doi{10.1016/j.jeconom.2007.08.017} +#' +#' Koop, G., León-González, R., & Strachan R. W. (2010). Efficient posterior +#' simulation for cointegrated models with priors on the cointegration space. +#' \emph{Econometric Reviews, 29}(2), 224--242. +#' \doi{10.1080/07474930903382208} +#' +#' Koop, G., León-González, R., & Strachan R. W. (2011). Bayesian inference in +#' a time varying cointegration model. \emph{Journal of Econometrics, 165}(2), 210--220. +#' \doi{10.1016/j.jeconom.2011.07.007} +#' +#' Korobilis, D. (2013). VAR forecasting using Bayesian variable selection. +#' \emph{Journal of Applied Econometrics, 28}(2), 204--230. \doi{10.1002/jae.1271} +#' +#' Lütkepohl, H. (2006). \emph{New introduction to multiple time series analysis} (2nd ed.). Berlin: Springer. +#' +#' +#' @export +add_priors.gvecsubmodels <- function(object, ..., + coef = list(v_i = 1, v_i_det = 0.1, shape = 3, rate = 0.0001, rate_det = 0.01), + coint = list(v_i = 0, p_tau_i = 1, shape = 3, rate = 0.0001, rho = 0.999), + sigma = list(df = "k", scale = 1, mu = 0, v_i = 0.01, sigma_h = 0.05), + ssvs = NULL, + bvs = NULL){ + + # Checks - Coefficient priors ---- + if (!is.null(coef)) { + if (!is.null(coef[["v_i"]])) { + if (coef[["v_i"]] < 0) { + stop("Argument 'v_i' must be at least 0.") + } + # Define "v_i_det" if not specified (needed for a check later) + if (is.null(coef[["v_i_det"]])) { + coef[["v_i_det"]] <- coef[["v_i"]] + } + } else { + if (!any(c("minnesota", "ssvs") %in% names(coef))) { + stop("If 'coef$v_i' is not specified, at least 'coef$minnesota' or 'coef$ssvs' must be specified.") + } + } + } + + if (!is.null(coef[["const"]])) { + if ("character" %in% class(coef[["const"]])) { + if (!coef[["const"]] %in% c("first", "mean")) { + stop("Invalid specificatin of coef$const.") + } + } + } + + if (length(coint) < 2) { + stop("Argument 'coint' must be at least of length 2.") + } + + # Checks - Error priors ---- + if (length(sigma) < 2) { + stop("Argument 'sigma' must be at least of length 2.") + } else { + error_prior <- NULL + if (any(unlist(lapply(object, function(x) {x[["model"]][["sv"]]})))) { # Check for SV + if (any(!c("mu", "v_i", "shape", "rate") %in% names(sigma))) { + stop("Missing prior specifications for stochastic volatility prior.") + } + if ("character" %in% class(sigma[["shape"]])) { + stop("Argument sigma$shape may not be a character when used with SV.") + } + error_prior <- "sv" + } else { + + if (all(c("df", "scale", "shape", "rate") %in% names(sigma))) { + error_prior <- "both" + } else { + if (all(c("shape", "rate") %in% names(sigma))) { + error_prior <- "gamma" + } + if (all(c("df", "scale") %in% names(sigma))) { + error_prior <- "wishart" + } + } + + if (is.null(error_prior)) { + stop("Invalid specification for argument 'sigma'.") + } + if (error_prior == "wishart" & any(unlist(lapply(object, function(x) {x$model$structural})))) { + stop("Structural models may not use a Wishart prior. Consider specifying argument 'sigma$shape' instead.") + } + + if (error_prior == "wishart") { + if (sigma[["df"]] < 0) { + stop("Argument 'sigma$df' must be at least 0.") + } + if (sigma[["scale"]] <= 0) { + stop("Argument 'sigma$scale' must be larger than 0.") + } + } + if (error_prior == "gamma") { + if (sigma[["shape"]] < 0) { + stop("Argument 'sigma$shape' must be at least 0.") + } + if (sigma[["rate"]] <= 0) { + stop("Argument 'sigma$rate' must be larger than 0.") + } + } + } + } + + # Check Minnesota ---- + minnesota <- FALSE # Minnesota prior? + if (!is.null(coef[["minnesota"]])) { + minnesota <- TRUE + if (!"list" %in% class(coef[["minnesota"]])) { + stop("Argument coeff$minnesota must be a list.") + } + if (!all(c("kappa0", "kappa1", "kappa3") %in% names(coef[["minnesota"]]))) { + stop("Argument coeff$minnesota must contain at least the elements 'kappa0', 'kappa1' and 'kappa3'.") + } + } + + # Check SSVS ---- + use_ssvs <- FALSE + use_ssvs_error <- FALSE + use_ssvs_semi <- FALSE + if (!is.null(ssvs)) { + + if (error_prior == "sv") { + stop("SSVS is not supported for stochastic volatility models. You could use BVS instead.") + } + + if (is.null(ssvs[["inprior"]])) { + stop("Argument 'ssvs$inprior' must be specified for SSVS.") + } + if (is.null(ssvs[["tau"]]) & is.null(ssvs[["semiautomatic"]])) { + stop("Either argument 'ssvs$tau' or 'ssvs$semiautomatic' must be specified for SSVS.") + } + if (is.null(ssvs[["exclude_det"]])) { + ssvs[["exclude_det"]] <- FALSE + } + # In case ssvs is specified, check if the semi-automatic approach of + # George et al. (2008) should be used + if (!is.null(ssvs[["semiautomatic"]])) { + use_ssvs_semi <- TRUE + } + + use_ssvs <- TRUE + if (minnesota) { + minnesota <- FALSE + warning("Minnesota prior specification overwritten by SSVS.") + } + + if (!is.null(ssvs[["covar"]])) { + if (ssvs[["covar"]]) { + if (error_prior == "wishart") { + stop("If SSVS should be applied to error covariances, argument 'sigma$shape' and 'sigma$rate' must be specified.") + } + use_ssvs_error <- TRUE + } + if (is.null(ssvs[["tau"]])) { + stop("If SSVS should be applied to error covariances, argument 'ssvs$tau' must be specified.") + } + } + } + + # Prior a la Korobilis 2013 + use_bvs <- FALSE + use_bvs_error <- FALSE + if (!is.null(bvs)) { + use_bvs <- TRUE + if (is.null(bvs[["inprior"]])) { + stop("If BVS should be applied, argument 'bvs$inprior' must be specified.") + } + if (is.null(bvs[["exclude_det"]])) { + bvs[["exclude_det"]] <- FALSE + } + if (!is.null(bvs[["covar"]])) { + if (bvs[["covar"]]) { + if (error_prior == "wishart") { + stop("If BVS should be applied to error covariances, argument 'sigma$shape' must be specified.") + } + use_bvs_error <- TRUE + } + } + if (coef[["v_i"]] == 0 | (coef[["v_i_det"]] == 0 & !bvs[["exclude_det"]])) { + warning("Using BVS with an uninformative prior is not recommended.") + } + } + + if (use_ssvs & use_bvs) { + stop("SSVS and BVS cannot be applied at the same time.") + } + + if (error_prior == "wishart" & (use_ssvs_error | use_bvs_error)) { + stop("Wishart prior not allowed when BVS or SSVS are applied to covariance matrix.") + } + + varsel_covar <- use_ssvs_error | use_bvs_error + + # Generate priors for each country ---- + for (i in 1:length(object)) { + + # Get model specs to obtain total number of coeffs + k_domestic <- length(object[[i]][["model"]][["domestic"]][["variables"]]) + p_domestic <- object[[i]][["model"]][["domestic"]][["lags"]] + + if (k_domestic == 1 & (use_ssvs_error | use_bvs_error)) { + stop("BVS or SSVS cannot be applied to covarianc matrix when there is only one endogenous variable.") + } + + k_foreign <- length(object[[i]][["model"]][["foreign"]][["variables"]]) + p_foreign <- object[[i]][["model"]][["foreign"]][["lags"]] + + global <- !is.null(object[[i]][["model"]][["global"]]) + if (global) { + k_global <- length(object[[i]][["model"]][["global"]][["variables"]]) + p_global <- object[[i]][["model"]][["global"]][["lags"]] + } else { + k_global <- 0 + p_global <- 0 + } + + # Substract lag from domestic model for VEC + p_domestic <- p_domestic - 1 + + # Total # of non-deterministic coefficients + tot_par <- k_domestic * (k_domestic * p_domestic + k_foreign * p_foreign + k_global * p_global) + + # Add number of unrestricted deterministic terms + n_det <- 0 + if (!is.null(object[[i]][["model"]][["deterministic"]])){ + if (!is.null(object[[i]][["model"]][["deterministic"]][["unrestricted"]])) { + n_det <- length(object[[i]][["model"]][["deterministic"]][["unrestricted"]]) * k_domestic + } + } + + tot_par <- tot_par + n_det + + covar <- FALSE + if (!is.null(sigma[["covar"]])) { + covar <- sigma[["covar"]] + } + structural <- object[[i]][["model"]][["structural"]] + if (covar & structural) { + stop("Error covariances and structural coefficients cannot be estimated at the same time.") + } + sv <- object[[i]][["model"]][["sv"]] + n_struct <- 0 + n_z <- NCOL(object[[i]][["data"]][["SUR"]]) + if (object[[i]][["model"]][["rank"]] > 0) { + n_w <- NCOL(object[[i]][["data"]][["W"]]) + n_z <- n_z - k_domestic * n_w + } + + if (structural & k_domestic > 1) { + n_struct <- (k_domestic - 1) * k_domestic / 2 + tot_par <- tot_par + n_struct + } + + # Cointegration (constant) ---- + n_ect <- k_domestic * (k_domestic + k_foreign) + if (global) { + n_ect <- n_ect + k_domestic * k_global + } + if (!is.null(object[[i]][["model"]][["deterministic"]][["restricted"]])) { + n_ect <- n_ect + k_domestic * length(object[[i]][["model"]][["deterministic"]][["restricted"]]) + } + + r_temp <- object[[i]][["model"]][["rank"]] + if (r_temp > 0) { + + n_alpha <- r_temp * k_domestic + n_beta <- r_temp * n_ect / k_domestic + + if (object[[i]][["model"]][["tvp"]]) { + rho <- coint[["rho"]] + if (is.null(rho)) { + rho <- .999 + } else { + if (rho >= 1) { + stop("rho must be smaller than 1.") + } + if (rho < .8) { + warning("Value of rho appears very small.") + } + } + object[[i]][["priors"]][["cointegration"]][["alpha"]] <- list(mu = matrix(0, n_alpha), + v_i = diag(1 / (1 - rho^2), n_alpha), + shape = matrix(coint[["shape"]], n_alpha), + rate = matrix(coint[["rate"]], n_alpha)) + object[[i]][["priors"]][["cointegration"]][["beta"]] <- list(mu = matrix(0, n_beta), + v_i = diag(1 - rho^2, n_beta)) + } else { + object[[i]][["priors"]][["cointegration"]] <- list(v_i = coint[["v_i"]], + p_tau_i = diag(coint[["p_tau_i"]], n_ect / k_domestic)) + } + } else { + if (r_temp == 0 & tot_par == 0) { + warning("Model with zero cointegration rank and no non-cointegration regressors is skipped.") + } + } + + #### Non-cointegration (constant) #### + # Generate prior matrices ---- + if (tot_par > 0) { + + # Zero prior means + mu <- matrix(rep(0, tot_par - n_struct), k_domestic) + + # Prior for intercept terms + if (n_det > 0) { + if (!is.null(coef[["const"]])) { + + pos <- which(dimnames(object[[i]][["data"]][["X"]])[[2]] == "const") + + if (length(pos) == 1) { + if ("character" %in% class(coef[["const"]])) { + if (coef[["const"]] == "first") { + mu[, pos] <- object[[i]][["data"]][["Y"]][1, ] + } + if (coef[["const"]] == "mean") { + mu[, pos] <- colMeans(object[[i]][["data"]][["Y"]]) + } + } + if ("numeric" %in% class(coef[["const"]])) { + mu[, pos] <- coef[["const"]] + } + } + } + } + + mu <- matrix(mu) + + if (structural) { + mu <- rbind(mu, matrix(0, n_struct)) + } + + object[[i]][["priors"]][["noncointegration"]] <- list(mu = mu) + + # Minnesota prior here + # Obtain OLS estimates for the calculation of the + # Minnesota prior or the semi-automatic SSVS approach + if (minnesota | use_ssvs_semi) { + # Get data + y <- t(object[[i]][["data"]][["Y"]]) + x <- t(cbind(object[[i]][["data"]][["W"]], object[[i]][["data"]][["X"]])) + k_ect <- ncol(object[[i]][["data"]][["W"]]) + + tt <- ncol(y) + ols_sigma <- y %*% (diag(1, tt) - t(x) %*% solve(tcrossprod(x)) %*% x) %*% t(y) / (tt - nrow(x)) + + if (minnesota) { + # Determine positions of deterministic terms for calculation of sigma + pos_det <- NULL + if (!is.null(object[[i]][["model"]][["deterministic"]])) { + if (!is.null(object[[i]][["model"]][["deterministic"]][["restricted"]])) { + pos_det <- c(pos_det, k_domestic + k_foreign + k_global + 1:length(object[[i]][["model"]][["deterministic"]][["restricted"]])) + } + if (!is.null(object[[i]][["model"]][["deterministic"]][["unrestricted"]])) { + pos_det <- c(pos_det, k_ect + k_domestic * p_domestic + k_foreign * p_foreign + k_global * p_global + 1:length(object[[i]][["model"]][["deterministic"]][["unrestricted"]])) + } + } + + # Obtain sigmas for V_i via estimation of AR model for each endogenous variable + s_domestic <- diag(0, k_domestic) + for (j in 1:k_domestic) { + if (p_domestic > 0) { + pos <- c(j, k_ect + j + k_domestic * ((1:p_domestic) - 1), pos_det) + } else { + pos <- c(j, pos_det) + } + y_temp <- matrix(y[j, ], 1) + x_temp <- matrix(x[pos, ], length(pos)) + s_domestic[j, j] <- y_temp %*% (diag(1, tt) - t(x_temp) %*% solve(tcrossprod(x_temp)) %*% x_temp) %*% t(y_temp) / (tt - length(pos)) + } + s_domestic <- sqrt(diag(s_domestic)) # Residual standard deviations (OLS) + + # Generate prior matrices ---- + + # Minnesota prior ---- + V <- matrix(rep(NA, tot_par - n_struct), k_domestic) # Set up matrix for variances + + # Endogenous variables + if (p_domestic > 0) { + for (r in 1:p_domestic) { + for (l in 1:k_domestic) { + for (j in 1:k_domestic) { + if (l == j) { + V[l, (r - 1) * k_domestic + j] <- coef$minnesota$kappa0 / r^2 + } else { + V[l, (r - 1) * k_domestic + j] <- coef$minnesota$kappa0 * coef$minnesota$kappa1 / r^2 * s_domestic[l]^2 / s_domestic[j]^2 + } + } + } + } + } + + # Weakly exogenous variables + s_foreign <- sqrt(apply(matrix(x[k_ect + k_domestic * p_domestic + 1:k_foreign, ], k_foreign), 1, stats::var)) + + for (r in 1:p_foreign) { + for (l in 1:k_domestic) { + for (j in 1:k_foreign) { + # Note that r starts at 1, so that this is equivalent to l + 1 + if (is.null(coef$minnesota$kappa2)) { + V[l, p_domestic * k_domestic + (r - 1) * k_foreign + j] <- coef$minnesota$kappa0 * coef$minnesota$kappa3 * s_domestic[l]^2 + } else { + V[l, p_domestic * k_domestic + (r - 1) * k_foreign + j] <- coef$minnesota$kappa0 * coef$minnesota$kappa2 / r^2 * s_domestic[l]^2 / s_foreign[j]^2 + } + } + } + } + + # Global variables + if (global) { + s_global <- sqrt(apply(matrix(x[k_ect + k_domestic * p_domestic + k_foreign * p_foreign + 1:k_global, ], k_global), 1, stats::var)) + for (r in 1:p_global) { + for (l in 1:k_domestic) { + for (j in 1:k_global) { + # Note that r starts at 1, so that this is equivalent to l + 1 + if (is.null(coef$minnesota$kappa2)) { + V[l, k_domestic * p_domestic + k_foreign * p_foreign + (r - 1) * k_global + j] <- coef$minnesota$kappa0 * coef$minnesota$kappa3 * s_domestic[l]^2 + } else { + V[l, k_domestic * p_domestic + k_foreign * p_foreign + (r - 1) * k_global + j] <- coef$minnesota$kappa0 * coef$minnesota$kappa2 / r^2 * s_domestic[l]^2 / s_global[j]^2 + } + } + } + } + } + + # Restrict prior variances + if (!is.null(coef[["max_var"]])) { + if (any(stats::na.omit(c(V)) > coef[["max_var"]])) { + V[which(V > coef[["max_var"]])] <- coef[["max_var"]] + } + } + + # Deterministic variables + if (!is.null(object[[i]][["data"]][["deterministic"]][["unrestricted"]])){ + V[, -(1:(k_domestic * p_domestic + k_foreign * p_foreign + k_global * p_global))] <- coef$minnesota$kappa0 * coef$minnesota$kappa3 * s_domestic^2 + } + + V <- matrix(V) + + # Structural parameters + if (structural & k_domestic > 1) { + V_struct <- matrix(NA, k_domestic, k_domestic) + for (j in 1:(k_domestic - 1)) { + V_struct[(j + 1):k_domestic, j] <- coef$minnesota$kappa0 * coef$minnesota$kappa1 * s_domestic[(j + 1):k_domestic]^2 / s_domestic[j]^2 + } + V_struct <- matrix(V_struct[lower.tri(V_struct)]) + V <- rbind(V, V_struct) + } + + v_i <- diag(c(1 / V)) + + object[[i]][["priors"]][["noncointegration"]][["v_i"]] <- v_i + } # End of minnesota condition + } # Ende of minnesota/ssvs_semi condition + + # Inclusion priors + if (use_ssvs | use_bvs) { + inprior <- matrix(NA, k_domestic, (tot_par - n_struct) / k_domestic) + include <- matrix(1:tot_par) + if (use_ssvs) { + prob <- ssvs[["inprior"]] + kappa <- ssvs[["minnesota"]] + exclude_det <- ssvs[["exclude_det"]] + } + if (use_bvs) { + prob <- bvs[["inprior"]] + kappa <- bvs[["minnesota"]] + exclude_det <- bvs[["exclude_det"]] + } + + # For Minnesota-like inclusion parameters + if (!is.null(kappa)) { + # Domestic + if (p_domestic > 0) { + for (r in 1:p_domestic) { + inprior[, (r - 1) * k_domestic + 1:k_domestic] <- kappa[2] / r + if (k_domestic > 1) { + diag(inprior[, (r - 1) * k_domestic + 1:k_domestic]) <- kappa[1] / r + } else { + inprior[, (r - 1) * k_domestic + 1] <- kappa[1] / r + } + } + } + + # Foreign + if (k_foreign > 0) { + inprior[, p_domestic * k_domestic + 1:k_foreign] <- kappa[3] + if (p_foreign > 1) { + for (r in 1:(p_foreign - 1)) { + inprior[, k_domestic * p_domestic + k_foreign + (r - 1) * k_foreign + 1:k_foreign] <- kappa[3] / (1 + r) + } + } + } + + # Global + if (global) { + inprior[, k_domestic * p_domestic + k_foreign * p_foreign + 1:k_global] <- kappa[3] + if (p_global > 1) { + for (r in 1:(p_global - 1)) { + inprior[, k_domestic * p_domestic + k_foreign * p_foreign + k_global + (r - 1) * k_global + 1:k_global] <- kappa[3] / (1 + r) + } + } + } + + if (n_det > 0) { + inprior[, k_domestic * p_domestic + k_foreign * p_foreign + k_global * p_global + 1:(n_det / k_domestic)] <- kappa[4] + } + } else { + inprior[,] <- prob + } + + inprior <- matrix(inprior) + if (structural & k_domestic > 1) { + inprior <- rbind(inprior, matrix(prob, n_struct)) + } + + if (n_det > 0 & exclude_det) { + pos_det <- tot_par - n_det - n_struct + 1:n_det + include <- matrix(include[-pos_det]) + } + + # SSVS + if (use_ssvs) { + if (use_ssvs_semi) { + # Semiautomatic approach + cov_ols <- kronecker(solve(tcrossprod(x)), ols_sigma) + se_ols <- sqrt(diag(cov_ols)) # OLS standard errors + + se_ols <- se_ols[-(1:k_ect)] + se_ols <- matrix(se_ols) + + tau0 <- se_ols * ssvs[["semiautomatic"]][1] # Prior if excluded + tau1 <- se_ols * ssvs[["semiautomatic"]][2] # Prior if included + + if (structural & k_domestic > 1) { + warning("Semiautomatic approach for SSVS not available for structural variables. Using values of argument 'ssvs$tau' instead.") + tau0 <- rbind(tau0, matrix(ssvs[["tau"]][1], n_struct)) + tau1 <- rbind(tau1, matrix(ssvs[["tau"]][2], n_struct)) + } + } else { + tau0 <- matrix(rep(ssvs[["tau"]][1], tot_par)) + tau1 <- matrix(rep(ssvs[["tau"]][2], tot_par)) + } + + object[[i]][["model"]][["varselect"]] <- "SSVS" + + + object[[i]][["priors"]][["noncointegration"]][["v_i"]] <- diag(1 / tau1[, 1]^2, tot_par) + object[[i]][["priors"]][["noncointegration"]][["ssvs"]][["inprior"]] <- inprior + object[[i]][["priors"]][["noncointegration"]][["ssvs"]][["include"]] <- include + object[[i]][["priors"]][["noncointegration"]][["ssvs"]][["tau0"]] <- tau0 + object[[i]][["priors"]][["noncointegration"]][["ssvs"]][["tau1"]] <- tau1 + + } + } + + # Regular prior ---- + if (!minnesota & !use_ssvs) { + v_i <- diag(coef[["v_i"]], tot_par) + # Add priors for deterministic terms if they were specified + if (n_det > 0 & !is.null(coef[["v_i_det"]])) { + diag(v_i)[tot_par - n_struct - n_det + 1:n_det] <- coef[["v_i_det"]] + } + + object[[i]][["priors"]][["noncointegration"]][["v_i"]] <- v_i + } + + # BVS + if (use_bvs) { + object[[i]][["model"]][["varselect"]] <- "BVS" + + object[[i]][["priors"]][["noncointegration"]][["bvs"]][["inprior"]] <- inprior + object[[i]][["priors"]][["noncointegration"]][["bvs"]][["include"]] <- include + } + + ## TVP prior - variances of the state equations ---- + if (object[[i]][["model"]][["tvp"]]) { + object[[i]][["priors"]][["noncointegration"]][["shape"]] <- matrix(coef[["shape"]], tot_par) + object[[i]][["priors"]][["noncointegration"]][["rate"]] <- matrix(coef[["rate"]], tot_par) + if (n_det > 0 & !is.null(coef[["rate_det"]])) { + object[[i]][["priors"]][["noncointegration"]][["rate"]][tot_par - n_struct - n_det + 1:n_det, ] <- coef[["rate_det"]] + } + } + } + + ## Covar priors ---- + if (!structural & covar & k_domestic > 1) { + + n_covar <- k_domestic * (k_domestic - 1) / 2 + object[[i]][["priors"]][["psi"]][["mu"]] <- matrix(0, n_covar) + object[[i]][["priors"]][["psi"]][["v_i"]] <- diag(coef[["v_i"]], n_covar) + if (object[[i]][["model"]][["tvp"]]) { + object[[i]][["priors"]][["psi"]][["shape"]] <- matrix(coef[["shape"]], n_covar) + object[[i]][["priors"]][["psi"]][["rate"]] <- matrix(coef[["rate"]], n_covar) + } + + if (use_ssvs_error) { + object[[i]][["priors"]][["psi"]][["ssvs"]][["inprior"]] <- matrix(ssvs[["inprior"]], n_covar) + object[[i]][["priors"]][["psi"]][["ssvs"]][["include"]] <- matrix(1:n_covar) + object[[i]][["priors"]][["psi"]][["ssvs"]][["tau0"]] <- matrix(ssvs[["tau"]][1], n_covar) + object[[i]][["priors"]][["psi"]][["ssvs"]][["tau1"]] <- matrix(ssvs[["tau"]][2], n_covar) + } + + if (use_bvs_error) { + object[[i]][["priors"]][["psi"]][["bvs"]][["inprior"]] <- matrix(bvs[["inprior"]], n_covar) + object[[i]][["priors"]][["psi"]][["bvs"]][["include"]] <- matrix(1:n_covar) + } + + } + + # Error term ---- + if (sv) { + + object[[i]][["priors"]][["sigma"]][["mu"]] <- matrix(sigma[["mu"]], k_domestic) + object[[i]][["priors"]][["sigma"]][["v_i"]] <- diag(sigma[["v_i"]], k_domestic) + object[[i]][["priors"]][["sigma"]][["shape"]] <- matrix(sigma[["shape"]], k_domestic) + object[[i]][["priors"]][["sigma"]][["rate"]] <- matrix(sigma[["rate"]], k_domestic) + + } else { + + if (error_prior == "wishart" | (error_prior == "both" & !structural)) { + object[[i]][["priors"]][["sigma"]][["type"]] <- "wishart" + help_df <- sigma[["df"]] + object[[i]][["priors"]][["sigma"]][["df"]] <- NA_real_ + object[[i]][["priors"]][["sigma"]][["scale"]] = diag(sigma[["scale"]], k_domestic) + } + + if (error_prior == "gamma" | (error_prior == "both" & structural)) { + object[[i]][["priors"]][["sigma"]][["type"]] <- "gamma" + help_df <- sigma[["shape"]] + object[[i]][["priors"]][["sigma"]][["shape"]] <- NA_real_ + object[[i]][["priors"]][["sigma"]][["rate"]] = matrix(sigma[["rate"]], k_domestic) + } + + if ("character" %in% class(help_df)) { + if (grepl("k", help_df)) { + k <- k_domestic + # Transform character specification to expression and evaluate + help_df <- eval(parse(text = help_df)) + rm(k) + } else { + stop("Use no other letter than 'k' in 'sigma$df' to indicate the number of endogenous variables.") + } + } + + if (help_df < 0) { + stop("Current specification implies a negative prior degree of\nfreedom or shape parameter of the error term.") + } + + # Add rank to degrees of freedom for cointegration model + if (!is.na(object[[i]][["model"]][["rank"]])) { + help_df <- help_df + r_temp + } + if (error_prior == "wishart" | (error_prior == "both" & !structural)) { + object[[i]][["priors"]][["sigma"]][["df"]] <- help_df + } + if (error_prior == "gamma" | (error_prior == "both" & structural)) { + object[[i]][["priors"]][["sigma"]][["shape"]] <- matrix(help_df, k_domestic) + } + } + + # Initial values ---- + y <- t(object[[i]][["data"]][["Y"]]) + # If there are enough observations obtain OLS estimates for first draw... + if (tot_par < NCOL(y)) { + z <- object[[i]][["data"]][["SUR"]] + ols <- solve(crossprod(z)) %*% crossprod(z, matrix(y)) + object[[i]][["initial"]][["noncointegration"]][["gamma"]] <- matrix(ols[-c(1:n_ect),]) + u <- matrix(matrix(y) - z %*% ols, NROW(y)) # Residuals for initial var-covar matrix draw + } else { + # ... if not, use the prior mean. + if (tot_par > 0) { + object[[i]][["initial"]][["noncointegration"]][["gamma"]] <- mu + } + u <- y - matrix(apply(y, 1, mean), nrow = NROW(y), ncol = NCOL(y)) # Residuals for initial var-covar matrix draw + } + if (r_temp > 0) { + beta <- matrix(0, n_ect / k_domestic, r_temp) + beta[1:r_temp, 1:r_temp] <- diag(1, r_temp) + object[[i]][["initial"]][["cointegration"]][["alpha"]] <- matrix(0, n_alpha) + object[[i]][["initial"]][["cointegration"]][["beta"]] <- beta + if (object[[i]][["model"]][["tvp"]]) { + object[[i]][["initial"]][["cointegration"]][["rho"]] <- rho + object[[i]][["initial"]][["cointegration"]][["sigma_alpha_i"]] <- diag(c(1 / object[[i]][["priors"]][["cointegration"]][["alpha"]][["rate"]])) + } + } + if (object[[i]][["model"]][["tvp"]]) { + object[[i]][["initial"]][["noncointegration"]][["sigma_gamma_i"]] <- diag(c(1 / object[[i]][["priors"]][["noncointegration"]][["rate"]]), tot_par) + } + if (covar) { + y_covar <- kronecker(-t(u), diag(1, k_domestic)) + pos <- NULL + for (j in 1:k_domestic) {pos <- c(pos, (j - 1) * k_domestic + 1:j)} + y_covar <- y_covar[, -pos] + psi <- solve(crossprod(y_covar)) %*% crossprod(y_covar, matrix(u)) + object[[i]][["initial"]][["psi"]][["psi"]] <- psi + object[[i]][["initial"]][["psi"]][["sigma_psi_i"]] <- diag(1 / coef[["rate"]], nrow(psi)) + Psi <- diag(1, k_domestic) + for (j in 2:k_domestic) { + Psi[j, 1:(j - 1)] <- t(psi[((j - 2) * (j - 1) / 2) + 1:(j - 1), 1]) + } + u <- Psi %*% u + } + u <- apply(u, 1, stats::var) + if (object[[i]][["model"]][["sv"]]) { + object[[i]][["initial"]][["sigma"]][["h"]] <- log(matrix(u, nrow = NCOL(y), ncol = NROW(y), byrow = TRUE)) + object[[i]][["initial"]][["sigma"]][["sigma_h"]] <- matrix(sigma[["sigma_h"]], NROW(y)) + object[[i]][["initial"]][["sigma"]][["constant"]] <- matrix(sigma[["constant"]], NROW(y)) + } else { + object[[i]][["initial"]][["sigma"]][["sigma_i"]] <- diag(1 / u, NROW(y)) + dimnames(object[[i]][["initial"]][["sigma"]][["sigma_i"]]) <- list(dimnames(y)[[2]], dimnames(y)[[2]]) + } + + } + return(object) +} diff --git a/R/bvecxpost.R b/R/bvecxpost.R new file mode 100644 index 0000000..9342fe9 --- /dev/null +++ b/R/bvecxpost.R @@ -0,0 +1,359 @@ +#' Posterior Simulation Country-Specific VECX Models of a GVAR Model +#' +#' Produces draws from the posterior distributions of Bayesian VECX models. +#' +#' @param object an object of class \code{"gvecsubmodels"}, usually, a result of a call to \code{\link{create_models}} +#' in combination with \code{\link{add_priors}}. +#' +#' @details The function implements posterior simulation algorithms proposed in Koop et al. (2010) +#' and Koop et al. (2011), which place identifying restrictions on the cointegration space. +#' Both algorithms are able to employ Bayesian variable selection (BVS) as proposed in Korobilis (2013). +#' The algorithm of Koop et al. (2010) is also able to employ stochastic search variable selection (SSVS) +#' as proposed by Geroge et al. (2008). +#' Both SSVS and BVS can also be applied to the covariances of the error term. However, the algorithms +#' cannot be applied to cointegration related coefficients, i.e. to the loading matrix \eqn{\alpha} or +#' the cointegration matrix \eqn{beta}. +#' +#' The implementation primarily follows the description in Koop et al. (2010). Chan et al. (2019), +#' George et al. (2008) and Korobilis (2013) were used to implement the variable selection algorithms. +#' For all approaches the SUR form of a VEC model is used to obtain posterior draws. The algorithm is implemented +#' in C++ to reduce calculation time. +#' +#' The function also supports structural BVEC models, where the structural coefficients are estimated from +#' contemporary endogenous variables, which corresponds to the so-called (A-model). Currently, only +#' specifications are supported, where the structural matrix contains ones on its diagonal and all lower +#' triangular elements are freely estimated. Since posterior draws are obtained based on the SUR form of +#' the VEC model, the structural coefficients are drawn jointly with the other coefficients. No identifying +#' restrictions are made regarding the cointegration matrix. +#' +#' @return An object of class \code{"bvar"}. +#' +#' @references +#' +#' Chan, J., Koop, G., Poirier, D. J., & Tobias J. L. (2019). \emph{Bayesian econometric methods} +#' (2nd ed.). Cambridge: Cambridge University Press. +#' +#' George, E. I., Sun, D., & Ni, S. (2008). Bayesian stochastic search for VAR model +#' restrictions. \emph{Journal of Econometrics, 142}(1), 553--580. +#' \doi{10.1016/j.jeconom.2007.08.017} +#' +#' Koop, G., León-González, R., & Strachan R. W. (2010). Efficient posterior +#' simulation for cointegrated models with priors on the cointegration space. +#' \emph{Econometric Reviews, 29}(2), 224--242. +#' \doi{10.1080/07474930903382208} +#' +#' Koop, G., León-González, R., & Strachan R. W. (2011). Bayesian inference in +#' a time varying cointegration model. \emph{Journal of Econometrics, 165}(2), 210--220. +#' \doi{10.1016/j.jeconom.2011.07.007} +#' +#' Korobilis, D. (2013). VAR forecasting using Bayesian variable selection. +#' \emph{Journal of Applied Econometrics, 28}(2), 204--230. \doi{10.1002/jae.1271} +#' +#' @export +bvecxpost <- function(object) { + + # Use C++ functions to draw posteriors + if (object[["model"]][["tvp"]]) { + object <- .bgvectvpalg(object) + } else { + object <- .bgvecalg(object) + } + + k_dom <- NCOL(object[["data"]][["Y"]]) + k_for <- length(object[["model"]][["foreign"]][["variables"]]) + if (is.null(object[["model"]][["global"]])) { + k_glo <- 0 + } else { + k_glo <- length(object[["model"]][["global"]][["variables"]]) + } + if (is.null(object[["model"]][["deterministic"]][["restricted"]])) { + k_det_r <- 0 + } else { + k_det_r <- length(object[["model"]][["deterministic"]][["restricted"]]) + } + if (is.null(object[["model"]][["deterministic"]][["unrestricted"]])) { + k_det_ur <- 0 + } else { + k_det_ur <- length(object[["model"]][["deterministic"]][["unrestricted"]]) + } + tt <- NROW(object[["data"]][["Y"]]) + + tvp_a0 <- FALSE + tvp_alpha <- FALSE + tvp_beta_dom <- FALSE + tvp_beta_for <- FALSE + tvp_beta_glo <- FALSE + tvp_beta_det <- FALSE + tvp_gamma_dom <- FALSE + tvp_gamma_for <- FALSE + tvp_gamma_glo <- FALSE + tvp_gamma_det <- FALSE + tvp_sigma <- FALSE + structural <- FALSE + + mc_spec <- c(object[["model"]][["burnin"]] + 1, object[["model"]][["burnin"]] + object[["model"]][["iterations"]], 1) + + if (!is.null(object[["posteriors"]][["sigma"]][["lambda"]])) { + sigma_lambda <- matrix(diag(NA_real_, k_dom), k_dom * k_dom, object[["model"]][["iterations"]]) + sigma_lambda[which(lower.tri(diag(1, k_dom))), ] <- object[["posteriors"]][["sigma"]][["lambda"]] + sigma_lambda[which(upper.tri(diag(1, k_dom))), ] <- object[["posteriors"]][["sigma"]][["lambda"]] + object[["posteriors"]][["sigma"]][["lambda"]] <- sigma_lambda + rm(sigma_lambda) + } + + #### Combine coefficients #### + + A0 <- NULL + if (object[["model"]][["structural"]]) { + + pos <- which(lower.tri(diag(1, k_dom))) + draws <- object[["model"]][["iterations"]] + + if (is.list(object[["posteriors"]][["a0"]])) { + + if ("coeffs" %in% names(object[["posteriors"]][["a0"]])) { + if (object[["model"]][["tvp"]]) { + A0[["coeffs"]] <- matrix(diag(1, k_dom), k_dom * k_dom * tt, draws) + A0[["coeffs"]][rep(0:(tt - 1) * k_dom * k_dom, each = length(pos)) + rep(pos, tt), ] <- object[["posteriors"]][["a0"]][["coeffs"]] + } else { + A0[["coeffs"]] <- matrix(diag(1, k_dom), k_dom * k_dom, draws) + A0[["coeffs"]][rep(0:(draws - 1) * k_dom * k_dom, each = length(pos)) + pos ] <- object[["posteriors"]][["a0"]][["coeffs"]] + } + } + + if ("sigma" %in% names(object[["posteriors"]][["a0"]])) { + A0[["sigma"]] <- matrix(0, k_dom * k_dom, draws) + A0[["sigma"]][pos, ] <- object[["posteriors"]][["a0"]][["sigma"]] + } + + if ("lambda" %in% names(object[["posteriors"]][["a0"]])) { + A0[["lambda"]] <- matrix(diag(1, k_dom), k_dom * k_dom, draws) + A0[["lambda"]][pos, ] <- object[["posteriors"]][["a0"]][["lambda"]] + A0[["lambda"]][-pos, ] <- NA_real_ + } + + } else { + A0 <- matrix(diag(1, k_dom), k_dom * k_dom, draws) + A0[pos, ] <- object[["posteriors"]][["a0"]] + } + + object[["posteriors"]][["a0"]] <- NULL + + if(!is.null(A0)) { + if (is.list(A0)) { + if ("coeffs" %in% names(A0)) { + n_a0 <- nrow(A0[["coeffs"]]) + } + } else { + n_a0 <- nrow(A0) + } + if (n_a0 / tt >= 1) { + tvp_a0 <- TRUE + n_a0 <- n_a0 / tt + } + if (n_a0 %% (k_dom * k_dom) != 0) { + stop("Row number of coefficient draws of 'A0' is not k^2 or multiples thereof.") + } + structural <- TRUE + + object[["posteriors"]] <- c(object[["posteriors"]], .country_fill_helper(A0, tvp_a0, n_a0, tt, "a0")) + } + } + + # alpha ---- + if(!is.null(object[["posteriors"]][["alpha"]])) { + + if ("coeffs" %in% names(object[["posteriors"]][["alpha"]])) { + n_alpha <- nrow(object[["posteriors"]][["alpha"]][["coeffs"]]) + } + if ((n_alpha / tt) %% k_dom == 0 & n_alpha / tt >= 1) { + tvp_alpha <- TRUE + n_alpha <- n_alpha / tt + } + + alpha <- object[["posteriors"]][["alpha"]] + object[["posteriors"]][["alpha"]] <- NULL + object[["posteriors"]] <- c(object[["posteriors"]], .country_fill_helper(alpha, tvp_alpha, n_alpha, tt, "alpha")) + rm(alpha) + } + + # Beta domestic ---- + if(!is.null(object[["posteriors"]][["beta_dom"]])) { + + if ("coeffs" %in% names(object[["posteriors"]][["beta_dom"]])) { + n_beta_dom <- nrow(object[["posteriors"]][["beta_dom"]][["coeffs"]]) + } + if ((n_beta_dom / tt) %% k_dom == 0 & n_beta_dom / tt >= 1) { + tvp_beta_dom <- TRUE + n_beta_dom <- n_beta_dom / tt + } + + beta_dom <- object[["posteriors"]][["beta_dom"]] + object[["posteriors"]][["beta_dom"]] <- NULL + object[["posteriors"]] <- c(object[["posteriors"]], .country_fill_helper(beta_dom, tvp_beta_dom, n_beta_dom, tt, "beta_domestic")) + rm(beta_dom) + } else { + object[["posteriors"]][["beta_dom"]] <- NULL + } + + # Beta foreign ---- + if(!is.null(object[["posteriors"]][["beta_for"]])) { + + if ("coeffs" %in% names(object[["posteriors"]][["beta_for"]])) { + n_beta_for <- nrow(object[["posteriors"]][["beta_for"]][["coeffs"]]) + } + if ((n_beta_for / tt) %% k_for == 0 & n_beta_for / tt >= 1) { + tvp_beta_for <- TRUE + n_beta_for <- n_beta_for / tt + } + + beta_for <- object[["posteriors"]][["beta_for"]] + object[["posteriors"]][["beta_for"]] <- NULL + object[["posteriors"]] <- c(object[["posteriors"]], .country_fill_helper(beta_for, tvp_beta_for, n_beta_for, tt, "beta_foreign")) + rm(beta_for) + } else { + object[["posteriors"]][["beta_for"]] <- NULL + } + + # Beta global ---- + if(!is.null(object[["posteriors"]][["beta_glo"]])) { + + if ("coeffs" %in% names(object[["posteriors"]][["beta_glo"]])) { + n_beta_glo <- nrow(object[["posteriors"]][["beta_glo"]][["coeffs"]]) + } + if ((n_beta_glo / tt) %% k_glo == 0 & n_beta_glo / tt >= 1) { + tvp_beta_glo <- TRUE + n_beta_glo <- n_beta_glo / tt + } + + beta_glo <- object[["posteriors"]][["beta_glo"]] + object[["posteriors"]][["beta_glo"]] <- NULL + object[["posteriors"]] <- c(object[["posteriors"]], .country_fill_helper(beta_glo, tvp_beta_glo, n_beta_glo, tt, "beta_global")) + rm(beta_glo) + } else { + object[["posteriors"]][["beta_glo"]] <- NULL + } + + # Beta deterministic ---- + if(!is.null(object[["posteriors"]][["beta_det"]])) { + + if ("coeffs" %in% names(object[["posteriors"]][["beta_det"]])) { + n_beta_det <- nrow(object[["posteriors"]][["beta_det"]][["coeffs"]]) + } + if ((n_beta_det / tt) %% k_det_r == 0 & n_beta_det / tt >= 1) { + tvp_beta_det <- TRUE + n_beta_det <- n_beta_det / tt + } + + beta_d <- object[["posteriors"]][["beta_det"]] + object[["posteriors"]][["beta_det"]] <- NULL + object[["posteriors"]] <- c(object[["posteriors"]], .country_fill_helper(beta_d, tvp_beta_det, n_beta_det, tt, "beta_deterministic")) + rm(beta_d) + } else { + object[["posteriors"]][["beta_det"]] <- NULL + } + + # Gamma domestic ---- + if(!is.null(object[["posteriors"]][["gamma_dom"]])) { + if ("coeffs" %in% names(object[["posteriors"]][["gamma_dom"]])) { + n_gamma_dom <- nrow(object[["posteriors"]][["gamma_dom"]][["coeffs"]]) + } + if ((n_gamma_dom / tt) %% k_dom == 0 & n_gamma_dom / tt >= 1) { + tvp_gamma_dom <- TRUE + n_gamma_dom <- n_gamma_dom / tt + } + + gamma_dom <- object[["posteriors"]][["gamma_dom"]] + object[["posteriors"]][["gamma_dom"]] <- NULL + object[["posteriors"]] <- c(object[["posteriors"]], .country_fill_helper(gamma_dom, tvp_gamma_dom, n_gamma_dom, tt, "gamma_domestic")) + rm(gamma_dom) + } else { + object[["posteriors"]][["gamma_dom"]] <- NULL + } + + # Gamma foreign ---- + if(!is.null(object[["posteriors"]][["gamma_for"]])) { + if ("coeffs" %in% names(object[["posteriors"]][["gamma_for"]])) { + n_gamma_for <- nrow(object[["posteriors"]][["gamma_for"]][["coeffs"]]) + } + if ((n_gamma_for / tt) %% k_for == 0 & n_gamma_for / tt >= 1) { + tvp_gamma_for <- TRUE + n_gamma_for <- n_gamma_for / tt + } + + gamma_for <- object[["posteriors"]][["gamma_for"]] + object[["posteriors"]][["gamma_for"]] <- NULL + object[["posteriors"]] <- c(object[["posteriors"]], .country_fill_helper(gamma_for, tvp_gamma_for, n_gamma_for, tt, "gamma_foreign")) + rm(gamma_for) + } + + # Gamma global ---- + if(!is.null(object[["posteriors"]][["gamma_glo"]])) { + if ("coeffs" %in% names(object[["posteriors"]][["gamma_glo"]])) { + n_gamma_glo <- nrow(object[["posteriors"]][["gamma_glo"]][["coeffs"]]) + } + if ((n_gamma_glo / tt) %% k_glo == 0 & n_gamma_glo / tt >= 1) { + tvp_gamma_glo <- TRUE + n_gamma_glo <- n_gamma_glo / tt + } + + gamma_glo <- object[["posteriors"]][["gamma_glo"]] + object[["posteriors"]][["gamma_glo"]] <- NULL + object[["posteriors"]] <- c(object[["posteriors"]], .country_fill_helper(gamma_glo, tvp_gamma_glo, n_gamma_glo, tt, "gamma_global")) + rm(gamma_glo) + } else { + object[["posteriors"]][["gamma_glo"]] <- NULL + } + + # Gamma deterministic ---- + if(!is.null(object[["posteriors"]][["gamma_det"]])) { + if ("coeffs" %in% names(object[["posteriors"]][["gamma_det"]])) { + n_gamma_det <- nrow(object[["posteriors"]][["gamma_det"]][["coeffs"]]) + } + if ((n_gamma_det / tt) %% k_det_ur == 0 & n_gamma_det / tt >= 1) { + tvp_gamma_det <- TRUE + n_gamma_det <- n_gamma_det / tt + } + + gamma_det <- object[["posteriors"]][["gamma_det"]] + object[["posteriors"]][["gamma_det"]] <- NULL + object[["posteriors"]] <- c(object[["posteriors"]], .country_fill_helper(gamma_det, tvp_gamma_det, n_gamma_det, tt, "gamma_deterministic")) + rm(gamma_det) + } else { + object[["posteriors"]][["gamma_det"]] <- NULL + } + + if(!is.null(object[["posteriors"]][["sigma"]])) { + if ("coeffs" %in% names(object[["posteriors"]][["sigma"]])) { + n_sigma <- nrow(object[["posteriors"]][["sigma"]][["coeffs"]]) + } + if ((n_sigma / tt) %% k_dom == 0 & n_sigma / tt >= 1) { + tvp_sigma <- TRUE + n_sigma <- n_sigma / tt + } + if (n_sigma %% (k_dom * k_dom) != 0) { + stop("Row number of coefficient draws of 'Sigma' is not k^2 or multiples thereof.") + } + + sigma <- object[["posteriors"]][["sigma"]] + object[["posteriors"]][["sigma"]] <- NULL + object[["posteriors"]] <- c(object[["posteriors"]], .country_fill_helper(sigma, tvp_sigma, n_sigma, tt, "sigma")) + rm(sigma) + } + + vars <- c("a0", + "alpha", + "beta_domestic", "beta_foreign", "beta_global", "beta_determinisitc", + "gamma_domestic", "gamma_foreign", "gamma_global", "gamma_determinisitc", + "sigma") + + for (i in vars) { + if (!is.null(object[["posteriors"]][[i]])) { + attr(object[["posteriors"]][[i]], "mcpar") <- mc_spec + } + } + + class(object) <- c("ctryvecest", "list") + + return(object) +} \ No newline at end of file diff --git a/R/draw_posterior.gvecsubmodels.R b/R/draw_posterior.gvecsubmodels.R new file mode 100644 index 0000000..5b6509a --- /dev/null +++ b/R/draw_posterior.gvecsubmodels.R @@ -0,0 +1,89 @@ +#' Posterior simulation +#' +#' Estimates the country models of a Bayesian GVAR model. +#' +#' @param object a list of data and model specifications for country-specific +#' models, which should be passed on to function \code{FUN}. Usually, the +#' output of a call to \code{\link{create_models}} in combination with \code{\link{add_priors}}. +#' @param FUN the function to be applied to each country model in argument \code{object}. +#' If \code{NULL} (default), the internal functions \code{\link{bvecxpost}} is used. +#' @param mc.cores the number of cores to use, i.e. at most how many child +#' processes will be run simultaneously. The option is initialized from +#' environment variable MC_CORES if set. Must be at least one, and +#' parallelization requires at least two cores. +#' @param ctry character vector specifying for which countries posterior distributions +#' should be drawn. If \code{NULL} (default), draws are generated for all country models. +#' @param ... further arguments passed to or from other methods. +#' +#' @return An object of class \code{bgvecest}, which contains a list of data, +#' model specifications, priors, coefficient draws and information +#' criteria for each estimated country model. +#' +#' @export +draw_posterior.gvecsubmodels <- function(object, ..., FUN = NULL, mc.cores = NULL, ctry = NULL){ + + # If 'ctry' is specified, reduce list to relevant elements + if (!is.null(ctry)) { + pos <- which(names(object) %in% ctry) + temp <- list() + for (i in 1:length(pos)) { + temp[[i]] <- object[[pos[i]]] + names(temp)[i] <- names(object)[pos[i]] + } + rm(object) + object <- temp + rm(temp) + names_obj <- names(object) + } + + names_obj <- names(object) + names_temp <- names_obj + if (length(unique(names_temp)) != length(names_temp)) { + for (i in unique(names_temp)) { + pos_temp <- which(names_temp == i) + temp <- names_temp[pos_temp] + id <- paste0("0000", 1:length(temp)) + id <- substring(id, nchar(id) - 3, nchar(id)) + temp <- paste0(id, "-", temp) + names_temp[pos_temp] <- temp + } + } + names(object) <- names_temp + + # Print estimation information + cat("Estimating submodels...\n") + if (is.null(mc.cores)) { + object <- lapply(object, .posterior_gvecsubmodels, use = FUN) + } else { + object <- parallel::mclapply(object, .posterior_gvecsubmodels, use = FUN, + mc.cores = mc.cores, mc.preschedule = FALSE) + } + + names(object) <- names_obj + + class(object) <- c("bgvecest", "list") + + return(object) +} + +# Helper function to implement try() functionality +.posterior_gvecsubmodels <- function(object, use) { + + # Save specs in case Gibbs sampler fails + model <- list(data = object[["data"]], + model = object[["model"]]) + + if (is.null(use)) { + object <- try(bvecxpost(object)) + } else { + # Apply own function + object <- try(use(object)) + } + + # Produce something if estimation fails + if (inherits(object, "try-error")) { + object <- c(object, c(model, list(error = TRUE))) + } + + return(object) +} \ No newline at end of file diff --git a/R/plot.bgvecest.R b/R/plot.bgvecest.R new file mode 100644 index 0000000..2472264 --- /dev/null +++ b/R/plot.bgvecest.R @@ -0,0 +1,30 @@ +#' Plotting Draws of a VECX Submodel of a GVAR Model +#' +#' A plot function for objects of class \code{"bgvecest"}. +#' +#' @param x an object of class \code{"bgvecest"}, usually, a result of a call to \code{\link{draw_posterior}}. +#' @param ci interval used to calculate credible bands for time-varying parameters. +#' @param type either \code{"hist"} (default) for histograms, \code{"trace"} for a trace plot, +#' or \code{"boxplot"} for a boxplot. Only used for parameter draws of constant coefficients. +#' @param ctry character. Name of the element in argument \code{x}, for which posterior draws +#' should be plotted. If \code{NULL} (default), all submodels are used. +#' @param ... further graphical parameters. +#' +#' @export +plot.bgvecest <- function(x, ci = 0.95, type = "hist", ctry = NULL, ...) { + + pos <- 1:length(x) + if (!is.null(ctry)) { + pos <- which(names(x) %in% ctry) + if (length(pos) == 0) { + stop("There is no output for the specified country.") + } + } + + for (i in pos) { + plot(x[[i]], ci = ci, type = type, ctry = names(x)[i], ...) + } + +} + + diff --git a/R/plot.ctryvecest.R b/R/plot.ctryvecest.R new file mode 100644 index 0000000..ea4bb87 --- /dev/null +++ b/R/plot.ctryvecest.R @@ -0,0 +1,485 @@ +#' Plotting Draws of a VECX Submodel of a GVAR Model +#' +#' A plot function for objects of class \code{"ctryvecest"} for visual inspection +#' of posterior draws. +#' +#' @param x an object of class \code{"ctryvecest"}, usually, a result of a call to \code{\link{draw_posterior}}. +#' @param ci interval used to calculate credible bands for time-varying parameters. +#' @param type either \code{"hist"} (default) for histograms, \code{"trace"} for a trace plot, +#' or \code{"boxplot"} for a boxplot. Only used for parameter draws of constant coefficients. +#' @param variables character vector of variables that should be plotted. Default is \code{"all"}. +#' Other options are \code{"domestic"}, \code{"foreign"}, \code{"global"}, \code{"deterministic"} +#' and \code{"sigma"}. +#' @param ctry character (optional). Name of the country, which will be shown in the title of the output. +#' @param ... further graphical parameters. +#' +#' @export +plot.ctryvecest <- function(x, ci = 0.95, type = "hist", variables = "all", ctry = NULL, ...) { + + if (!type %in% c("hist", "trace", "boxplot")) { + stop("Argument 'type' must be 'hist', 'trace' or 'boxplot'.") + } + + if (!variables %in% c("all", "domestic", "foreign", "deterministic", "sigma")) { + stop("Invalid specification of argument 'variables'.") + } + + x <- .create_pi_matrices(x) + + names_domestic <- x[["model"]][["domestic"]][["variables"]] + k_domestic <- length(x[["model"]][["domestic"]][["variables"]]) + p_domestic <- x[["model"]][["domestic"]][["lags"]] + names_foreign <- x[["model"]][["foreign"]][["variables"]] + k_foreign <- length(x[["model"]][["foreign"]][["variables"]]) + p_foreign <- x[["model"]][["foreign"]][["lags"]] + global <- !is.null(x[["model"]][["global"]][["variables"]]) + if (global) { + names_global <- x[["model"]][["global"]][["variables"]] + k_global <- length(x[["model"]][["global"]][["variables"]]) + s_global <- x[["model"]][["global"]][["lags"]] + } else { + k_global <- 0 + s_global <- 0 + } + names_det_r <- x[["model"]][["deterministic"]][["restricted"]] + k_det_r <- length(names_det_r) + names_det_ur <- x[["model"]][["deterministic"]][["unrestricted"]] + k_det_ur <- length(names_det_ur) + r <- x[["model"]][["rank"]] + + tt <- NROW(x[["data"]][["Y"]]) + tsp_info <- stats::tsp(x[["data"]][["Y"]]) + structural <- x[["model"]][["structural"]] + ci_low <- (1 - ci) / 2 + ci_high <- 1 - ci_low + y_names <- paste0("d.", dimnames(x[["data"]][["domestic"]])[[2]]) + x_names <- .get_regressor_names_vec(x) + lab_size <- .05 + mar_orig <- graphics::par("mar") + + # Calculate number of columns in final layout + n_tot <- 0 + if ("pi_domestic" %in% names(x[["posteriors"]])) { + n_tot <- n_tot + k_domestic + } + if ("pi_foreign" %in% names(x[["posteriors"]])) { + n_tot <- n_tot + k_foreign + } + if ("pi_global" %in% names(x[["posteriors"]])) { + n_tot <- n_tot + k_global + } + if ("pi_deterministic" %in% names(x[["posteriors"]])) { + n_tot <- n_tot + length(names_det_r) + } + if ("gamma_domestic" %in% names(x[["posteriors"]])) { + n_tot <- n_tot + k_domestic * (p_domestic - 1) + } + if ("gamma_foreign" %in% names(x[["posteriors"]])) { + n_tot <- n_tot + k_foreign * p_foreign + } + if ("gamma_global" %in% names(x[["posteriors"]])) { + n_tot <- n_tot + k_global * s_global + } + if ("gamma_deterministic" %in% names(x[["posteriors"]])) { + n_tot <- n_tot + length(names_det_ur) + } + if (!is.null(x[["posteriors"]][["a0"]])) { + n_tot <- n_tot + k_domestic + } + if ("sigma" %in% names(x[["posteriors"]])) { + n_tot <- n_tot + k_domestic + } + + if (variables == "domestic") { + n_tot <- k_domestic * p_domestic + if (p_domestic > 1) { + x_names <- x_names[c(1:k_domestic, k_domestic + k_foreign + k_global + k_det_r + 1:(k_domestic * (p_domestic - 1)))] + } else { + x_names <- x_names[1:k_domestic] + } + } + if (variables == "foreign") { + n_tot <- k_foreign * (p_foreign + 1) + x_names <- x_names[c(k_domestic + 1:k_foreign, k_domestic * p_domestic + k_foreign + k_global + k_det_r + 1:(k_foreign * (p_foreign)))] + } + if (variables == "global") { + n_tot <- k_global * (s_global + 1) + x_names <- x_names[c(k_domestic + k_foreign + 1:k_global, + k_domestic * p_domestic + k_foreign * (p_foreign + 1) + k_global + k_det_r + 1:(k_global * (s_global)))] + } + if (variables == "deterministic") { + n_tot <- k_det_r + k_det_ur + pos_temp <- NULL + if (k_det_r > 0) { + pos_temp <- k_domestic + k_foreign + k_global + 1:k_det_r + } + if (k_det_ur > 0) { + pos_temp <- c(pos_temp, k_domestic * p_domestic + k_foreign * (p_foreign + 1) + k_global * (s_global + 1) + k_det_r + 1:k_det_ur) + } + x_names <- x_names[pos_temp] + } + if (variables == "sigma") { + n_tot <- k_domestic + } + + # Define layout + mat <- matrix(NA_integer_, + k_domestic + 2 , # k_domestic + title row + column names + n_tot + 1) # n_tot + row names + mat[1, ] <- 1 + mat[-1, 1] <- c(0, 2:(k_domestic + 1)) + mat[2, -1] <- (k_domestic + 1) + 1:n_tot + mat[-(1:2), -1] <- matrix(1:(k_domestic * n_tot) + k_domestic + n_tot + 1, k_domestic, n_tot) + graphics::layout(mat, + widths = c(lab_size, rep((1 - lab_size) / n_tot, n_tot)), + heights = c(.07, lab_size, rep((1 - lab_size) / k_domestic, k_domestic))) + + # Title + title_text <- "Bayesian " + tvp <- x[["model"]][["tvp"]] + if (tvp) { + title_text <- paste0(title_text, "TVP-") + } + sv <- x[["model"]][["sv"]] + if (sv) { + title_text <- paste0(title_text, "SV-") + } + if (x[["model"]][["structural"]]) { + title_text <- paste0(title_text, "S") + } + title_text <- paste0(title_text, "VECX submodel (r = ", r, ")") + if (!is.null(ctry)) { + title_text <- paste0(title_text, " - ", ctry) + } + + graphics::par(mar = c(0, 0, 0, 0)) + graphics::plot.new(); graphics::text(0.5, 0.5, labels = title_text, cex = 1.5) + # Fill rows + graphics::par(mar = c(3, 0, 0, 0)) + for (j in y_names) { + graphics::plot.new(); graphics::text(0.5, 0.5, labels = j, adj = 0.5) + } + # Fill columns + graphics::par(mar = c(0, 0, 0, 0)) + if (variables %in% c("all", "domestic", "foreign", "global", "deterministic")) { + for (j in x_names) { + graphics::plot.new(); graphics::text(0.5, 0.5, labels = j, adj = 0.5) + } + } + if (variables %in% c("all", "sigma")) { + for (j in y_names) { + graphics::plot.new(); graphics::text(0.5, 0.5, labels = paste0("Sigma\n", j), adj = 0.5) + } + } + + graphics::par(mar = c(3, 2.1, .5, 1)) + + if (variables %in% c("all", "domestic")) { + if (!is.null(x[["posteriors"]][["pi_domestic"]])) { + var_pos <- 1:(k_domestic * k_domestic) + if (tvp) { + for (j in var_pos) { + draws <- .tvpribbon(x[["posteriors"]][["pi_domestic"]], j, ci_low, ci_high) + stats::tsp(draws) <- tsp_info + stats::ts.plot(draws, xlab = "") + } + } else { + for (j in var_pos) { + if (all(x[["posteriors"]][["pi_domestic"]][, j] == x[["posteriors"]][["pi_domestic"]][1, j])) { + graphics::plot.new(); graphics::text(0.5, 0.5, labels = x[["posteriors"]][["pi_domestic"]][1, j], adj = 0.5) + } else { + if (type == "hist") { + graphics::hist(x[["posteriors"]][["pi_domestic"]][, j], plot = TRUE, main = NA) + } + if (type == "trace") { + stats::ts.plot(x[["posteriors"]][["pi_domestic"]][, j], xlab = "") + } + if (type == "boxplot") { + graphics::boxplot(x[["posteriors"]][["pi_domestic"]][, j]) + } + } + } + } + } + } + + if (variables %in% c("all", "foreign")) { + if (!is.null(x[["posteriors"]][["pi_foreign"]])) { + var_pos <- 1:(k_domestic * k_foreign) + if (tvp) { + for (j in var_pos) { + draws <- .tvpribbon(x[["posteriors"]][["pi_foreign"]], j, ci_low, ci_high) + stats::tsp(draws) <- tsp_info + stats::ts.plot(draws, xlab = "") + } + } else { + for (j in var_pos) { + if (all(x[["posteriors"]][["pi_foreign"]][, j] == x[["posteriors"]][["pi_foreign"]][1, j])) { + graphics::plot.new(); graphics::text(0.5, 0.5, labels = x[["posteriors"]][["pi_foreign"]][1, j], adj = 0.5) + } else { + if (type == "hist") { + graphics::hist(x[["posteriors"]][["pi_foreign"]][, j], plot = TRUE, main = NA) + } + if (type == "trace") { + stats::ts.plot(x[["posteriors"]][["pi_foreign"]][, j], xlab = "") + } + if (type == "boxplot") { + graphics::boxplot(x[["posteriors"]][["pi_foreign"]][, j]) + } + } + } + } + } + } + + if (variables %in% c("all", "global")) { + if (!is.null(x[["posteriors"]][["pi_global"]])) { + var_pos <- 1:(k_domestic * k_global) + if (tvp) { + for (j in var_pos) { + draws <- .tvpribbon(x[["posteriors"]][["pi_global"]], j, ci_low, ci_high) + stats::tsp(draws) <- tsp_info + stats::ts.plot(draws, xlab = "") + } + } else { + for (j in var_pos) { + if (all(x[["posteriors"]][["pi_global"]][, j] == x[["posteriors"]][["pi_global"]][1, j])) { + graphics::plot.new(); graphics::text(0.5, 0.5, labels = x[["posteriors"]][["pi_global"]][1, j], adj = 0.5) + } else { + if (type == "hist") { + graphics::hist(x[["posteriors"]][["pi_global"]][, j], plot = TRUE, main = NA) + } + if (type == "trace") { + stats::ts.plot(x[["posteriors"]][["pi_global"]][, j], xlab = "") + } + if (type == "boxplot") { + graphics::boxplot(x[["posteriors"]][["pi_global"]][, j]) + } + } + } + } + } + } + + if (variables %in% c("all", "deterministic")) { + if (!is.null(x[["posteriors"]][["pi_deterministic"]])) { + var_pos <- 1:(k_domestic * length(names_det_r)) + if (tvp) { + for (j in var_pos) { + draws <- .tvpribbon(x[["posteriors"]][["pi_deterministic"]], j, ci_low, ci_high) + stats::tsp(draws) <- tsp_info + stats::ts.plot(draws, xlab = "") + } + } else { + for (j in var_pos) { + if (all(x[["posteriors"]][["pi_deterministic"]][, j] == x[["posteriors"]][["pi_deterministic"]][1, j])) { + graphics::plot.new(); graphics::text(0.5, 0.5, labels = x[["posteriors"]][["pi_deterministic"]][1, j], adj = 0.5) + } else { + if (type == "hist") { + graphics::hist(x[["posteriors"]][["pi_deterministic"]][, j], plot = TRUE, main = NA) + } + if (type == "trace") { + stats::ts.plot(x[["posteriors"]][["pi_deterministic"]][, j], xlab = "") + } + if (type == "boxplot") { + graphics::boxplot(x[["posteriors"]][["pi_deterministic"]][, j]) + } + } + } + } + } + } + + if (variables %in% c("all", "domestic")) { + if (!is.null(x[["posteriors"]][["gamma_domestic"]])) { + if (p_domestic > 1) { + for (i in 1:(p_domestic - 1)) { + var_pos <- ((i - 1) * k_domestic * k_domestic) + 1:(k_domestic * k_domestic) + if (tvp) { + for (j in var_pos) { + draws <- .tvpribbon(x[["posteriors"]][["gamma_domestic"]], j, ci_low, ci_high) + stats::tsp(draws) <- tsp_info + stats::ts.plot(draws, xlab = "") + } + } else { + for (j in var_pos) { + if (all(x[["posteriors"]][["gamma_domestic"]][, j] == x[["posteriors"]][["gamma_domestic"]][1, j])) { + graphics::plot.new(); graphics::text(0.5, 0.5, labels = x[["posteriors"]][["gamma_domestic"]][1, j], adj = 0.5) + } else { + if (type == "hist") { + graphics::hist(x[["posteriors"]][["gamma_domestic"]][, j], plot = TRUE, main = NA) + } + if (type == "trace") { + stats::ts.plot(x[["posteriors"]][["gamma_domestic"]][, j], xlab = "") + } + if (type == "boxplot") { + graphics::boxplot(x[["posteriors"]][["gamma_domestic"]][, j]) + } + } + } + } + } + } + } + } + + if (variables %in% c("all", "foreign")) { + if (!is.null(x[["posteriors"]][["gamma_foreign"]])) { + for (i in 1:p_foreign) { + var_pos <- ((i - 1) * k_domestic * k_foreign) + 1:(k_domestic * k_foreign) + if (tvp) { + for (j in var_pos) { + draws <- .tvpribbon(x[["posteriors"]][["gamma_foreign"]], j, ci_low, ci_high) + stats::tsp(draws) <- tsp_info + stats::ts.plot(draws, xlab = "") + } + } else { + for (j in var_pos) { + if (all(x[["posteriors"]][["gamma_foreign"]][, j] == x[["posteriors"]][["gamma_foreign"]][1, j])) { + graphics::plot.new(); graphics::text(0.5, 0.5, labels = x[["posteriors"]][["gamma_foreign"]][1, j], adj = 0.5) + } else { + if (type == "hist") { + graphics::hist(x[["posteriors"]][["gamma_foreign"]][, j], plot = TRUE, main = NA) + } + if (type == "trace") { + stats::ts.plot(x[["posteriors"]][["gamma_foreign"]][, j], xlab = "") + } + if (type == "boxplot") { + graphics::boxplot(x[["posteriors"]][["gamma_foreign"]][, j]) + } + } + } + } + } + } + } + + if (variables %in% c("all", "global")) { + if (!is.null(x[["posteriors"]][["gamma_global"]])) { + for (i in 1:s_global) { + var_pos <- ((i - 1) * k_domestic * k_global) + 1:(k_domestic * k_global) + if (tvp) { + for (j in var_pos) { + draws <- .tvpribbon(x[["posteriors"]][["gamma_global"]], j, ci_low, ci_high) + stats::tsp(draws) <- tsp_info + stats::ts.plot(draws, xlab = "") + } + } else { + for (j in var_pos) { + if (all(x[["posteriors"]][["gamma_global"]][, j] == x[["posteriors"]][["gamma_global"]][1, j])) { + graphics::plot.new(); graphics::text(0.5, 0.5, labels = x[["posteriors"]][["gamma_global"]][1, j], adj = 0.5) + } else { + if (type == "hist") { + graphics::hist(x[["posteriors"]][["gamma_global"]][, j], plot = TRUE, main = NA) + } + if (type == "trace") { + stats::ts.plot(x[["posteriors"]][["gamma_global"]][, j], xlab = "") + } + if (type == "boxplot") { + graphics::boxplot(x[["posteriors"]][["gamma_global"]][, j]) + } + } + } + } + } + } + } + + if (variables %in% c("all", "deterministic")) { + if (!is.null(x[["posteriors"]][["gamma_deterministic"]])) { + if (tvp) { + for (j in 1:NCOL(x[["posteriors"]][["gamma_deterministic"]][[1]])) { + draws <- .tvpribbon(x[["posteriors"]][["gamma_deterministic"]], j, ci_low, ci_high) + stats::tsp(draws) <- tsp_info + stats::ts.plot(draws, xlab = "") + } + } else { + for (j in 1:NCOL(x[["posteriors"]][["gamma_deterministic"]])) { + if (all(x[["posteriors"]][["gamma_deterministic"]][, j] == x[["posteriors"]][["gamma_deterministic"]][1, j])) { + graphics::plot.new(); graphics::text(0.5, 0.5, labels = x[["posteriors"]][["gamma_deterministic"]][1, j], adj = 0.5) + } else { + if (type == "hist") { + graphics::hist(x[["posteriors"]][["gamma_deterministic"]][, j], plot = TRUE, main = NA) + } + if (type == "trace") { + stats::ts.plot(x[["posteriors"]][["gamma_deterministic"]][, j], xlab = "") + } + if (type == "boxplot") { + graphics::boxplot(x[["posteriors"]][["gamma_deterministic"]][, j]) + } + } + } + } + } + } + + if (variables %in% c("all")) { + if (!is.null(x[["posteriors"]][["a0"]])) { + if (tvp) { + for (j in 1:NCOL(x[["posteriors"]][["a0"]][[1]])) { + draws <- .tvpribbon(x[["posteriors"]][["a0"]], j, ci_low, ci_high) + if (all(draws[, 1] == draws[1, 1])) { + graphics::plot.new(); graphics::text(0.5, 0.5, labels = draws[1, 2], adj = 0.5) + } else { + stats::tsp(draws) <- tsp_info + stats::ts.plot(draws, xlab = "") + } + } + } else { + for (j in 1:NCOL(x[["posteriors"]][["a0"]])) { + if (all(x[["posteriors"]][["a0"]][, j] == x[["posteriors"]][["a0"]][1, j])) { + graphics::plot.new(); graphics::text(0.5, 0.5, labels = x[["posteriors"]][["a0"]][1, j], adj = 0.5) + } else { + if (type == "hist") { + graphics::hist(x[["posteriors"]][["a0"]][, j], plot = TRUE, main = NA) + } + if (type == "trace") { + stats::ts.plot(x[["posteriors"]][["a0"]][, j], xlab = "") + } + if (type == "boxplot") { + graphics::boxplot(x[["posteriors"]][["a0"]][, j]) + } + } + } + } + } + } + + if (variables %in% c("all", "sigma")) { + if (!is.null(x[["posteriors"]][["sigma"]])) { + var_pos <- 1:(k_domestic * k_domestic) + if (sv) { + for (j in var_pos) { + draws <- .tvpribbon(x[["posteriors"]][["sigma"]], j, ci_low, ci_high) + if (all(draws[, 1] == draws[1, 1])) { + graphics::plot.new(); graphics::text(0.5, 0.5, labels = draws[1, 2], adj = 0.5) + } else { + stats::tsp(draws) <- tsp_info + stats::ts.plot(draws, xlab = "") + } + } + } else { + for (j in var_pos) { + if (all(x[["posteriors"]][["sigma"]][, j] == x[["posteriors"]][["sigma"]][1, j])) { + graphics::plot.new(); graphics::text(0.5, 0.5, labels = x[["posteriors"]][["sigma"]][1, j], adj = 0.5) + } else { + if (type == "hist") { + graphics::hist(x[["posteriors"]][["sigma"]][, j], plot = TRUE, main = NA) + } + if (type == "trace") { + stats::ts.plot(x[["posteriors"]][["sigma"]][, j], xlab = "") + } + if (type == "boxplot") { + graphics::boxplot(x[["posteriors"]][["sigma"]][, j]) + } + } + } + } + } + } + + graphics::par(mar = mar_orig) + graphics::layout(matrix(1)) +} + + diff --git a/R/print.summary.bgvecest.R b/R/print.summary.bgvecest.R new file mode 100644 index 0000000..b9c0d69 --- /dev/null +++ b/R/print.summary.bgvecest.R @@ -0,0 +1,15 @@ +#' @include summary.bgvecest.R +#' +#' @export +#' @rdname summary.bgvecest +print.summary.bgvecest <- function(x, digits = max(3L, getOption("digits") - 3L), ...){ + for (i in 1:length(x)) { + if (i != 1) { + cat("--------------------------------------------------------------------------------\n") + cat("--------------------------------------------------------------------------------\n\n") + + } + cat("Submodel: ", names(x)[i], "\n") + print(x[[i]], digits = digits, ...) + } +} diff --git a/R/print.summary.ctryvecest.R b/R/print.summary.ctryvecest.R new file mode 100644 index 0000000..b54d28e --- /dev/null +++ b/R/print.summary.ctryvecest.R @@ -0,0 +1,113 @@ +#' @include summary.ctryvecest.R +#' +#' @export +#' @rdname summary.ctryvecest +print.summary.ctryvecest <- function(x, digits = max(3L, getOption("digits") - 3L), ...){ + + # Title + title_text <- "\nBayesian " + tvp <- x[["model"]][["tvp"]] + if (tvp) { + title_text <- paste0(title_text, "TVP-") + } + if (x[["model"]][["sv"]]) { + title_text <- paste0(title_text, "SV-") + } + if (x[["model"]][["structural"]]) { + title_text <- paste0(title_text, "S") + } + title_text <- paste0(title_text, "VECX submodel") + + cat(title_text, "\n") + + # Model + + if (is.null(x[["coefficients"]][["means"]])) { + regressors <- "nothing" + use_a <- FALSE + } else { + regressors <- paste(dimnames(x[["coefficients"]][["means"]])[[2]], collapse = " + ") + use_a <- TRUE + } + + cat("\nModel:\n\n", + paste("y ~ ", regressors, + sep = ""), "\n", sep = "") + + if (!is.null(x[["model"]][["period"]])) { + cat("\nPeriod:", x[["model"]][["period"]], "\n") + } + + k <- length(x[["model"]][["domestic"]][["variables"]]) + y_names <- x[["model"]][["domestic"]][["variables"]] + + # Coefficients per endogenous variable + + if (use_a) { + for (i in 1:k) { + temp <- cbind(x[["coefficients"]][["means"]][i, ], + x[["coefficients"]][["sd"]][i, ], + x[["coefficients"]][["naivesd"]][i, ], + x[["coefficients"]][["tssd"]][i, ], + x[["coefficients"]][["q_lower"]][i, ], + x[["coefficients"]][["median"]][i, ], + x[["coefficients"]][["q_upper"]][i, ]) + + dim_names_1 <- dimnames(x[["coefficients"]][["means"]])[[2]] + dim_names_2 <- c("Mean", "SD", "Naive SD", "Time-series SD", + x[["model"]][["ci"]][1], "50%", x[["model"]][["ci"]][2]) + + if ("lambda" %in% names(x[["coefficients"]])) { + temp <- cbind(temp, x[["coefficients"]][["lambda"]][i, ]) + dim_names_2 <- c(dim_names_2, "Incl. prob.") + } + + dimnames(temp)[[1]] <- dim_names_1 + dimnames(temp)[[2]] <- dim_names_2 + + cat("\nVariable:", y_names[i], "\n\n") + print.default(temp, quote = FALSE, right = TRUE, digits = digits) + } + } else { + cat("\n\nNo regressors.\n\n") + } + + # Error covariance matrix + + if (!is.null(x[["sigma"]])) { + x_names <- NULL + for (i in 1:k) { + x_names <- c(x_names , paste(dimnames(x[["sigma"]][["means"]])[[1]][i], dimnames(x[["sigma"]][["means"]])[[1]], sep = "_")) + } + + temp <- cbind(matrix(x[["sigma"]][["means"]]), + matrix(x[["sigma"]][["sd"]]), + matrix(x[["sigma"]][["naivesd"]]), + matrix(x[["sigma"]][["tssd"]]), + matrix(x[["sigma"]][["q_lower"]]), + matrix(x[["sigma"]][["median"]]), + matrix(x[["sigma"]][["q_upper"]])) + + dim_names_1 <- x_names + dim_names_2 <- c("Mean", "SD", "Naive SD", "Time-series SD", + x[["model"]][["ci"]][1], "50%", x[["model"]][["ci"]][2]) + + if ("lambda" %in% names(x[["sigma"]])) { + temp <- cbind(temp, matrix(x[["sigma"]][["lambda"]])) + dim_names_2 <- c(dim_names_2, "Incl. prob.") + } + + dimnames(temp) <- list(dim_names_1, + dim_names_2) + + if (k == 1) { + cat("\nVariance:\n\n") + } else { + cat("\nVariance-covariance matrix:\n\n") + } + print.default(temp, quote = FALSE, right = TRUE, digits = digits) + } + + cat("\n") + invisible(x) +} diff --git a/R/summary.bgvecest.R b/R/summary.bgvecest.R new file mode 100644 index 0000000..7ffcb49 --- /dev/null +++ b/R/summary.bgvecest.R @@ -0,0 +1,45 @@ +#' Summarising Country-Specific VECX Submodels of a GVAR Model +#' +#' summary method for class \code{"bgvecest"}. +#' +#' @param object an object of class \code{"bgvecest"}, usually, a result of a call to +#' \code{\link{draw_posterior}}. +#' @param ci a numeric between 0 and 1 specifying the probability of the credible band. +#' Defaults to 0.95. +#' @param period integer. Index of the period, for which the summary statistics should be generated. +#' Only used for TVP or SV models. Default is \code{NULL}, so that the posterior draws of the last time period +#' are used. +#' @param ctry character. Name of the element in argument \code{object}, for which summary +#' statistics should be obtained. If \code{NULL} (default), all country submodels are used. +#' @param x an object of class \code{"summary.bgvecest"}, usually, a result of a call to +#' \code{\link{summary.bgvecest}}. +#' @param digits the number of significant digits to use when printing. +#' @param ... further arguments passed to or from other methods. +#' +#' @return \code{summary.bgvecest} returns a list of class \code{"summary.bgvecest"}, +#' which contains summary statistics of country-specific submodels of class +#' \code{\link{summary.ctryvecest}}. +#' +#' @export +summary.bgvecest <- function(object, ci = .95, period = NULL, ctry = NULL, ...){ + + pos <- 1:length(object) + if (!is.null(ctry)) { + pos <- which(names(object) %in% ctry) + } + + if (length(pos) == 0) { + stop("Specified countries not available.") + } + + result <- NULL + for (i in pos) { + temp <- summary.ctryvecest(object[[i]], ci = ci, period = period, ...) + result <- c(result, list(temp)) + rm(temp) + } + names(result) <- names(object)[pos] + + class(result) <- c("summary.bgvecest", "list") + return(result) +} diff --git a/R/summary.ctryvecest.R b/R/summary.ctryvecest.R new file mode 100644 index 0000000..9ace33d --- /dev/null +++ b/R/summary.ctryvecest.R @@ -0,0 +1,215 @@ +#' Summarising Country-Specific VECX Submodels of a GVAR Model +#' +#' summary method for class \code{"ctryvecest"}. +#' +#' @param object an object of class \code{"ctryvecest"}, usually, a result of a call to \code{\link{draw_posterior}}. +#' @param ci a numeric between 0 and 1 specifying the probability of the credible band. +#' Defaults to 0.95. +#' @param period integer. Index of the period, for which the summary statistics should be generated. +#' Only used for TVP or SV models. Default is \code{NULL}, so that the posterior draws of the last time period +#' are used. +#' @param x an object of class \code{"summary.ctryvecest"}, usually, a result of a call to +#' \code{\link{summary.ctryvecest}}. +#' @param digits the number of significant digits to use when printing. +#' @param ... further arguments passed to or from other methods. +#' +#' @return \code{summary.ctryvecest} returns a list of class \code{"summary.ctryvecest"}, +#' which contains the following components: +#' \item{coefficients}{A list of various summary statistics of the posterior +#' draws of the VEC coefficients.} +#' \item{sigma}{A list of various summary statistics of the posterior +#' draws of the variance-covariance matrix.} +#' \item{model}{a list containing information on the model specifications.} +#' +#' @export +summary.ctryvecest <- function(object, ci = .95, period = NULL, ...){ + + # Get model specs ---- + names_domestic <- object[["model"]][["domestic"]][["variables"]] + k_domestic <- length(object[["model"]][["domestic"]][["variables"]]) + p_domestic <- object[["model"]][["domestic"]][["lags"]] + names_foreign <- object[["model"]][["foreign"]][["variables"]] + k_foreign <- length(object[["model"]][["foreign"]][["variables"]]) + p_foreign <- object[["model"]][["foreign"]][["lags"]] + names_global <- object[["model"]][["global"]][["variables"]] + global <- !is.null(object[["model"]][["global"]][["variables"]]) + r <- object[["model"]][["rank"]] + if (global) { + names_global <- object[["model"]][["global"]][["variables"]] + k_global <- length(object[["model"]][["global"]][["variables"]]) + s_global <- object[["model"]][["global"]][["lags"]] + } else { + k_global <- 0 + s_global <- 0 + } + k_det_r <- 0 + names_det_r <- object[["model"]][["deterministic"]][["restricted"]] + if (!is.null(names_det_r) & r > 0) { + k_det_r <- length(object[["model"]][["deterministic"]][["restricted"]]) + } + k_det_ur <- 0 + names_det_ur <- object[["model"]][["deterministic"]][["unrestricted"]] + if (!is.null(names_det_ur)) { + k_det_ur <- length(object[["model"]][["deterministic"]][["unrestricted"]]) + } + + tt <- NROW(object[["data"]][["Y"]]) + tvp <- object[["model"]][["tvp"]] + sv <- object[["model"]][["sv"]] + if (tvp | sv) { + if (is.null(period)) { + period <- tt + } else { + if (period > tt | period < 1) { + stop("Implausible specification of argument 'period'.") + } + } + period_long <- stats::time(object[["data"]][["Y"]])[period] + } else { + period_long <- NULL + } + + # Obtain variable names + x_names <- .get_regressor_names_vec(object) + dim_names <- list(names_domestic, x_names) + + # Non-error coefficients + means <- NULL + median <- NULL + sds <- NULL + naive_sd <- NULL + ts_sd <- NULL + + ci_low <- (1 - ci) / 2 + ci_high <- 1 - ci_low + q_low <- NULL + q_high <- NULL + + use_incl <- FALSE + if (any(grepl("lambda", names(object[["posteriors"]])))) { + use_incl <- TRUE + incl <- NULL + } + + object <- .create_pi_matrices(object) + + vars <- c("pi_domestic", "pi_foreign", "pi_global", "pi_deterministic", + "gamma_domestic", "gamma_foreign", "gamma_global", "gamma_deterministic", "a0") + for (i in vars) { + if (!is.null(object[["posteriors"]][[i]])) { + if (tvp) { + temp <- summary(object[["posteriors"]][[i]][[period]], quantiles = c(ci_low, .5, ci_high)) + } else { + temp <- summary(object[["posteriors"]][[i]], quantiles = c(ci_low, .5, ci_high)) + } + if ("numeric" %in% class(temp[["statistics"]])) { + means <- cbind(means, matrix(temp[["statistics"]]["Mean"], k_domestic)) + sds <- cbind(sds, matrix(temp[["statistics"]]["SD"], k_domestic)) + naive_sd <- cbind(naive_sd, matrix(temp[["statistics"]]["Naive SE"], k_domestic)) + ts_sd <- cbind(ts_sd, matrix(temp[["statistics"]]["Time-series SE"], k_domestic)) + q_low <- cbind(q_low, matrix(temp[["quantiles"]][1], k_domestic)) + median <- cbind(median, matrix(temp[["quantiles"]][2], k_domestic)) + q_high <- cbind(q_high, matrix(temp[["quantiles"]][3], k_domestic)) + } else { + means <- cbind(means, matrix(temp[["statistics"]][, "Mean"], k_domestic)) + sds <- cbind(sds, matrix(temp[["statistics"]][, "SD"], k_domestic)) + naive_sd <- cbind(naive_sd, matrix(temp[["statistics"]][, "Naive SE"], k_domestic)) + ts_sd <- cbind(ts_sd, matrix(temp[["statistics"]][, "Time-series SE"], k_domestic)) + q_low <- cbind(q_low, matrix(temp[["quantiles"]][, 1], k_domestic)) + median <- cbind(median, matrix(temp[["quantiles"]][, 2], k_domestic)) + q_high <- cbind(q_high, matrix(temp[["quantiles"]][, 3], k_domestic)) + } + if (use_incl) { + var_temp <- paste0(i, "_lambda") + if (var_temp %in% names(object[["posteriors"]])) { + incl <- cbind(incl, matrix(colMeans(object[["posteriors"]][[var_temp]]), k_domestic)) + } else { + incl <- cbind(incl, matrix(rep(NA_real_, ncol(object[["posteriors"]][[i]])), k_domestic)) + } + } + } + } + + if (!is.null(means)) { + dimnames(means) <- dim_names + dimnames(sds) <- dim_names + dimnames(naive_sd) <- dim_names + dimnames(ts_sd) <- dim_names + dimnames(q_low) <- dim_names + dimnames(median) <- dim_names + dimnames(q_high) <- dim_names + if (use_incl) { + dimnames(incl) <- dim_names + } + } + + result <- list(coefficients = list(means = means, + median = median, + sd = sds, + naivesd = naive_sd, + tssd = ts_sd, + q_lower = q_low, + q_upper = q_high)) + + if (use_incl) { + result[["coefficients"]][["lambda"]] = incl + } + + dim_names <- list(dim_names[[1]], dim_names[[1]]) + + # Error coefficients + if (!is.null(object[["posteriors"]][["sigma"]])) { + if (sv) { + temp <- summary(object[["posteriors"]][["sigma"]][[period]], quantiles = c(ci_low, .5, ci_high)) + } else { + temp <- summary(object[["posteriors"]][["sigma"]], quantiles = c(ci_low, .5, ci_high)) + } + if (k_domestic == 1) { + means <- matrix(temp[["statistics"]]["Mean"], k_domestic) + sds <- matrix(temp[["statistics"]]["SD"], k_domestic) + naive_sd <- matrix(temp[["statistics"]]["Naive SE"], k_domestic) + ts_sd <- matrix(temp[["statistics"]]["Time-series SE"], k_domestic) + q_low <- matrix(temp[["quantiles"]][1], k_domestic) + median <- matrix(temp[["quantiles"]][2], k_domestic) + q_high <- matrix(temp[["quantiles"]][3], k_domestic) + } else { + means <- matrix(temp[["statistics"]][, "Mean"], k_domestic) + sds <- matrix(temp[["statistics"]][, "SD"], k_domestic) + naive_sd <- matrix(temp[["statistics"]][, "Naive SE"], k_domestic) + ts_sd <- matrix(temp[["statistics"]][, "Time-series SE"], k_domestic) + q_low <- matrix(temp[["quantiles"]][, 1], k_domestic) + median <- matrix(temp[["quantiles"]][, 2], k_domestic) + q_high <- matrix(temp[["quantiles"]][, 3], k_domestic) + } + + + dimnames(means) <- dim_names + dimnames(sds) <- dim_names + dimnames(naive_sd) <- dim_names + dimnames(ts_sd) <- dim_names + dimnames(q_low) <- dim_names + dimnames(median) <- dim_names + dimnames(q_high) <- dim_names + + result[["sigma"]] <- list(means = means, + median = median, + sd = sds, + naivesd = naive_sd, + tssd = ts_sd, + q_lower = q_low, + q_upper = q_high) + + if ("sigma_lambda" %in% names(object[["posteriors"]])) { + incl <- matrix(colMeans(object[["posteriors"]][["sigma_lambda"]]), k_domestic) + dimnames(incl) <- dim_names + result[["sigma"]][["lambda"]] = incl + } + } + + result[["model"]] <- object[["model"]] + result[["model"]][["ci"]] <- paste(c(ci_low, ci_high) * 100, "%", sep = "") + result[["model"]][["period"]] <- period_long + + class(result) <- c("summary.ctryvecest", "list") + return(result) +} diff --git a/R/thin.bgvecest.R b/R/thin.bgvecest.R new file mode 100644 index 0000000..55dc8c7 --- /dev/null +++ b/R/thin.bgvecest.R @@ -0,0 +1,74 @@ +#' Thinning Posterior Draws +#' +#' Thins the MCMC posterior draws in an object of class \code{"bgvecest"}. +#' +#' @param x an object of class \code{"bgvecest"}. +#' @param thin an integer specifying the thinning interval between successive values of posterior draws. +#' @param ... further arguments passed to or from other methods. +#' +#' @return An object of class \code{"bgvecest"}. +#' +#' @export +thin.bgvecest <- function(x, thin = 10, ...) { + + vars <- c("a0", "a0_lambda", "a0_sigma", + "alpha", "alpha_lambda", "alpha_sigma", + "beta_domestic", "beta_domestic_lambda", "beta_domestic_sigma", + "beta_foreign", "beta_foreign_lambda", "beta_foreign_sigma", + "beta_global", "beta_global_lambda", "beta_global_sigma", + "beta_deterministic", "beta_deterministic_lambda", "beta_deterministic_sigma", + "gamma_domestic", "gamma_domestic_lambda", "gamma_domestic_sigma", + "gamma_foreign", "gamma_foreign_lambda", "gamma_foreign_sigma", + "gamma_global", "gamma_global_lambda", "gamma_global_sigma", + "gamma_deterministic", "gamma_deterministic_lambda", "gamma_deterministic_sigma", + "sigma", "sigma_lambda") + + draws <- NA + for (j in 1:length(x)) { + + if (!is.null(x[[j]][["error"]])) { + if (x[[j]][["error"]]) { + next + } + } + + for (i in vars) { + if (is.na(draws)) { + if (!is.null(x[[j]][["posteriors"]][[i]])) { + if (is.list(x[[j]][["posteriors"]][[i]])) { + draws <- nrow(x[[j]][["posteriors"]][[i]][[1]]) + } else { + draws <- nrow(x[[j]][["posteriors"]][[i]]) + } + } + } + } + } + + pos_thin <- seq(from = thin, to = draws, by = thin) + start <- pos_thin[1] + end <- pos_thin[length(pos_thin)] + + for (j in 1:length(x)) { + + if (!is.null(x[[j]][["error"]])) { + if (x[[j]][["error"]]) { + next + } + } + + for (i in vars) { + if (!is.null(x[[j]][["posteriors"]][[i]])) { + if (is.list(x[[j]][["posteriors"]][[i]])) { + for (k in 1:length(x[[j]][["posteriors"]][[i]])) { + x[[j]][["posteriors"]][[i]][[k]] <- coda::mcmc(as.matrix(x[[j]][["posteriors"]][[i]][[k]][pos_thin,]), start = start, end = end, thin = thin) + } + } else { + x[[j]][["posteriors"]][[i]] <- coda::mcmc(as.matrix(x[[j]][["posteriors"]][[i]][pos_thin,]), start = start, end = end, thin = thin) + } + } + } + } + + return(x) +} \ No newline at end of file diff --git a/man/add_priors.gvecsubmodels.Rd b/man/add_priors.gvecsubmodels.Rd new file mode 100644 index 0000000..516e88c --- /dev/null +++ b/man/add_priors.gvecsubmodels.Rd @@ -0,0 +1,188 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/add_priors.gvecsubmodels.R +\name{add_priors.gvecsubmodels} +\alias{add_priors.gvecsubmodels} +\title{Add Priors to Country-Specific VECX Models of a GVAR Model} +\usage{ +\method{add_priors}{gvecsubmodels}( + object, + ..., + coef = list(v_i = 1, v_i_det = 0.1, shape = 3, rate = 1e-04, rate_det = 0.01), + coint = list(v_i = 0, p_tau_i = 1, shape = 3, rate = 1e-04, rho = 0.999), + sigma = list(df = "k", scale = 1, mu = 0, v_i = 0.01, sigma_h = 0.05), + ssvs = NULL, + bvs = NULL +) +} +\arguments{ +\item{object}{a named list, usually, the output of a call to \code{\link{create_models}}.} + +\item{...}{further arguments passed to or from other methods.} + +\item{coef}{a named list of prior specifications for the coefficients of the country-specific +models. For the default specification all prior means are set to zero and the diagonal elements of +the inverse prior variance-covariance matrix are set to 1 for coefficients corresponding to non-deterministic +terms. For deterministic coefficients the prior variances are set to 10 via \code{v_i_det = 0.1}. +The variances need to be specified as precisions, i.e. as inverses of the variances. +For further specifications see 'Details'.} + +\item{coint}{a named list of prior specifications for cointegration coefficients of +country-specific VEC models. See 'Details'.} + +\item{sigma}{a named list of prior specifications for the error variance-covariance matrix +of the country models. For the default specification of an inverse Wishart distribution +the prior degrees of freedom are set to the number of endogenous variables and +the prior variances to 1. See 'Details'.} + +\item{ssvs}{a named list of prior specifications for the SSVS algorithm. See 'Details'.} + +\item{bvs}{a named list of prior specifications for the BVS algorithm. See 'Details'.} +} +\value{ +A list of country models. +} +\description{ +Adds prior specifications to a list of country models, which was produced by +the function \code{\link{create_models}}. +} +\details{ +The argument \code{coef} can contain the following elements +\describe{ + \item{\code{v_i}}{a numeric specifying the prior precision of the coefficients. Default is 1.} + \item{\code{v_i_det}}{a numeric specifying the prior precision of coefficients corresponding to deterministic terms. Default is 0.1.} + \item{\code{coint_var}}{a logical specifying whether the prior mean of the first own lag of an + endogenous variable in a VAR model should be set to 1. Default is \code{FALSE}.} + \item{\code{const}}{a numeric or character specifying the prior mean of coefficients, which correspond + to the intercept. If a numeric is provided, all prior means are set to this value. + If \code{const = "mean"}, the means of the series of endogenous variables are used as prior means. + If \code{const = "first"}, the first values of the series of endogenous variables are used as prior means.} + \item{\code{minnesota}}{a list of length 4 containing parameters for the calculation of + the Minnesota prior, where the element names are \code{kappa0}, \code{kappa1}, \code{kappa2} and \code{kappa3}. + For the endogenous variable \eqn{i} the prior variance of the \eqn{l}th + lag of regressor \eqn{j} is obtained as + \deqn{ \frac{\kappa_{0}}{l^2} \textrm{ for own lags of endogenous variables,}} + \deqn{ \frac{\kappa_{0} \kappa_{1}}{l^2} \frac{\sigma_{i}^2}{\sigma_{j}^2} \textrm{ for endogenous variables other than own lags,}} + \deqn{ \frac{\kappa_{0} \kappa_{2}}{(1 + l)^2} \frac{\sigma_{i}^2}{\sigma_{j}^2} \textrm{ for foreign and global exogenous variables,}} + \deqn{ \kappa_{0} \kappa_{3} \sigma_{i}^2 \textrm{ for deterministic terms,}} + where \eqn{\sigma_{i}} is the residual standard deviation of variable \eqn{i} of an unrestricted + LS estimate. For exogenous variables \eqn{\sigma_{i}} is the sample standard deviation. + If \code{kappa2 = NULL}, \eqn{\kappa_{0} \kappa_{3} \sigma_{i}^2} will be used for foreign + and global exogenous variables instead. + + For VEC models the function only provides priors for the non-cointegration part of the model. The + residual standard errors \eqn{\sigma_i} are based on an unrestricted LS regression of the + endogenous variables on the error correction term and the non-cointegration regressors.} + \item{\code{max_var}}{a numeric specifying the maximum prior variance that is allowed for + non-deterministic coefficients.} + \item{\code{shape}}{an integer specifying the prior degrees of freedom of the error term of the state equation. Default is 3.} + \item{\code{rate}}{a numeric specifying the prior error variance of the state equation. Default is 0.0001.} + \item{\code{rate_det}}{a numeric specifying the prior error variance of the state equation corresponding to deterministic terms. Default is 0.0001.} +} +If \code{minnesota} is specified, \code{v_i} and \code{v_i_det} are ignored. + +The argument \code{coint} can contain the following elements +\describe{ + \item{\code{coint_v_i}}{numeric between 0 and 1 specifying the shrinkage of the cointegration space prior. Default is 0.} + \item{\code{coint_p_tau_i}}{numeric of the diagonal elements of the inverse prior matrix of +the central location of the cointegration space \eqn{sp(\beta)}. Default is 1.} +} + +Argument \code{sigma} can contain the following elements: +\describe{ + \item{\code{df}}{an integer or character specifying the prior degrees of freedom of the error term. Only + used, if the prior is inverse Wishart. + Default is \code{"k"}, which indicates the amount of endogenous variables in the respective model. + \code{"k + 3"} can be used to set the prior to the amount of endogenous variables plus 3. If an integer + is provided, the degrees of freedom are set to this value in all models. + In all cases the rank \eqn{r} of the cointegration matrix is automatically added.} + \item{\code{scale}}{a numeric specifying the prior error variance of endogenous variables. + Default is 1.} + \item{\code{shape}}{a numeric or character specifying the prior shape parameter of the error term. Only + used, if the prior is inverse gamma or if time varying volatilities are estimated. + For models with constant volatility the default is \code{"k"}, which indicates the amount of endogenous + variables in the respective country model. \code{"k + 3"} can be used to set the prior to the amount of + endogenous variables plus 3. If a numeric is provided, the shape parameters are set to this value in all + models. For models with stochastic volatility this prior refers to the error variance of the state + equation.} + \item{\code{rate}}{a numeric specifying the prior rate parameter of the error term. Only used, if the + prior is inverse gamma or if time varying volatilities are estimated. For models with stochastic + volatility this prior refers to the error variance of the state equation.} + \item{\code{mu}}{numeric of the prior mean of the initial state of the log-volatilities. + Only used for models with time varying volatility.} + \item{\code{v_i}}{numeric of the prior precision of the initial state of the log-volatilities. + Only used for models with time varying volatility.} + \item{\code{sigma_h}}{numeric of the initial draw for the variance of the log-volatilities. + Only used for models with time varying volatility.} + \item{\code{constant}}{numeric of the constant, which is added before taking the log of the squared errors. + Only used for models with time varying volatility.} + \item{\code{covar}}{logical indicating whether error covariances should be estimated. Only used + in combination with an inverse gamma prior or stochastic volatility, for which \code{shape} and + \code{rate} must be specified.} +} +\code{df} and \code{scale} must be specified for an inverse Wishart prior. \code{shape} and \code{rate} +are required for an inverse gamma prior. For structural models or models with stochastic volatility +only a gamma prior specification is allowed. + +Argument \code{ssvs} can contain the following elements: +\describe{ + \item{\code{inprior}}{a numeric between 0 and 1 specifying the prior probability + of a variable to be included in the model.} + \item{\code{tau}}{a numeric vector of two elements containing the prior standard errors + of restricted variables (\eqn{\tau_0}) as its first element and unrestricted variables (\eqn{\tau_1}) + as its second.} + \item{\code{semiautomatic}}{an numeric vector of two elements containing the + factors by which the standard errors associated with an unconstrained least squares + estimate of the model are multiplied to obtain the prior standard errors + of restricted (\eqn{\tau_0}) and unrestricted (\eqn{\tau_1}) variables, respectively. + This is the semiautomatic approach described in George et al. (2008).} + \item{\code{covar}}{logical indicating if SSVS should also be applied to the error covariance matrix + as in George et al. (2008).} + \item{\code{exclude_det}}{logical indicating if deterministic terms should be excluded from the SSVS algorithm.} + \item{\code{minnesota}}{a numeric vector of length 4 containing parameters for the calculation of + the Minnesota-like inclusion priors. See below.} +} +Either \code{tau} or \code{semiautomatic} must be specified. + +The argument \code{bvs} can contain the following elements +\describe{ + \item{\code{inprior}}{a numeric between 0 and 1 specifying the prior probability + of a variable to be included in the model.} + \item{\code{covar}}{logical indicating if BVS should also be applied to the error covariance matrix.} + \item{\code{exclude_det}}{logical indicating if deterministic terms should be excluded from the BVS algorithm.} + \item{\code{minnesota}}{a numeric vector of length 4 containing parameters for the calculation of + the Minnesota-like inclusion priors. See below.} +} + +If either \code{ssvs$minnesota} or \code{bvs$minnesota} is specified, prior inclusion probabilities +are calculated in a Minnesota-like fashion as +\tabular{cl}{ +\eqn{\frac{\kappa_1}{l}} \tab for own lags of endogenous variables, \cr +\eqn{\frac{\kappa_2}{l}} \tab for other endogenous variables, \cr +\eqn{\frac{\kappa_3}{1 + l}} \tab for exogenous variables, \cr +\eqn{\kappa_{4}} \tab for deterministic variables, +} +for lag \eqn{l} with \eqn{\kappa_1}, \eqn{\kappa_2}, \eqn{\kappa_3}, \eqn{\kappa_4} as the first, second, +third and forth element in \code{ssvs$minnesota} or \code{bvs$minnesota}, respectively. +} +\references{ +Chan, J., Koop, G., Poirier, D. J., & Tobias J. L. (2019). \emph{Bayesian econometric methods} +(2nd ed.). Cambridge: Cambridge University Press. + +George, E. I., Sun, D., & Ni, S. (2008). Bayesian stochastic search for VAR model +restrictions. \emph{Journal of Econometrics, 142}(1), 553--580. +\doi{10.1016/j.jeconom.2007.08.017} + +Koop, G., León-González, R., & Strachan R. W. (2010). Efficient posterior +simulation for cointegrated models with priors on the cointegration space. +\emph{Econometric Reviews, 29}(2), 224--242. +\doi{10.1080/07474930903382208} + +Koop, G., León-González, R., & Strachan R. W. (2011). Bayesian inference in +a time varying cointegration model. \emph{Journal of Econometrics, 165}(2), 210--220. +\doi{10.1016/j.jeconom.2011.07.007} + +Korobilis, D. (2013). VAR forecasting using Bayesian variable selection. +\emph{Journal of Applied Econometrics, 28}(2), 204--230. \doi{10.1002/jae.1271} + +Lütkepohl, H. (2006). \emph{New introduction to multiple time series analysis} (2nd ed.). Berlin: Springer. +} diff --git a/man/bvecxpost.Rd b/man/bvecxpost.Rd new file mode 100644 index 0000000..d67a037 --- /dev/null +++ b/man/bvecxpost.Rd @@ -0,0 +1,60 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/bvecxpost.R +\name{bvecxpost} +\alias{bvecxpost} +\title{Posterior Simulation Country-Specific VECX Models of a GVAR Model} +\usage{ +bvecxpost(object) +} +\arguments{ +\item{object}{an object of class \code{"gvecsubmodels"}, usually, a result of a call to \code{\link{create_models}} +in combination with \code{\link{add_priors}}.} +} +\value{ +An object of class \code{"bvar"}. +} +\description{ +Produces draws from the posterior distributions of Bayesian VECX models. +} +\details{ +The function implements posterior simulation algorithms proposed in Koop et al. (2010) +and Koop et al. (2011), which place identifying restrictions on the cointegration space. +Both algorithms are able to employ Bayesian variable selection (BVS) as proposed in Korobilis (2013). +The algorithm of Koop et al. (2010) is also able to employ stochastic search variable selection (SSVS) +as proposed by Geroge et al. (2008). +Both SSVS and BVS can also be applied to the covariances of the error term. However, the algorithms +cannot be applied to cointegration related coefficients, i.e. to the loading matrix \eqn{\alpha} or +the cointegration matrix \eqn{beta}. + +The implementation primarily follows the description in Koop et al. (2010). Chan et al. (2019), +George et al. (2008) and Korobilis (2013) were used to implement the variable selection algorithms. +For all approaches the SUR form of a VEC model is used to obtain posterior draws. The algorithm is implemented +in C++ to reduce calculation time. + +The function also supports structural BVEC models, where the structural coefficients are estimated from +contemporary endogenous variables, which corresponds to the so-called (A-model). Currently, only +specifications are supported, where the structural matrix contains ones on its diagonal and all lower +triangular elements are freely estimated. Since posterior draws are obtained based on the SUR form of +the VEC model, the structural coefficients are drawn jointly with the other coefficients. No identifying +restrictions are made regarding the cointegration matrix. +} +\references{ +Chan, J., Koop, G., Poirier, D. J., & Tobias J. L. (2019). \emph{Bayesian econometric methods} +(2nd ed.). Cambridge: Cambridge University Press. + +George, E. I., Sun, D., & Ni, S. (2008). Bayesian stochastic search for VAR model +restrictions. \emph{Journal of Econometrics, 142}(1), 553--580. +\doi{10.1016/j.jeconom.2007.08.017} + +Koop, G., León-González, R., & Strachan R. W. (2010). Efficient posterior +simulation for cointegrated models with priors on the cointegration space. +\emph{Econometric Reviews, 29}(2), 224--242. +\doi{10.1080/07474930903382208} + +Koop, G., León-González, R., & Strachan R. W. (2011). Bayesian inference in +a time varying cointegration model. \emph{Journal of Econometrics, 165}(2), 210--220. +\doi{10.1016/j.jeconom.2011.07.007} + +Korobilis, D. (2013). VAR forecasting using Bayesian variable selection. +\emph{Journal of Applied Econometrics, 28}(2), 204--230. \doi{10.1002/jae.1271} +} diff --git a/man/draw_posterior.gvecsubmodels.Rd b/man/draw_posterior.gvecsubmodels.Rd new file mode 100644 index 0000000..461d2b1 --- /dev/null +++ b/man/draw_posterior.gvecsubmodels.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/draw_posterior.gvecsubmodels.R +\name{draw_posterior.gvecsubmodels} +\alias{draw_posterior.gvecsubmodels} +\title{Posterior simulation} +\usage{ +\method{draw_posterior}{gvecsubmodels}(object, ..., FUN = NULL, mc.cores = NULL, ctry = NULL) +} +\arguments{ +\item{object}{a list of data and model specifications for country-specific +models, which should be passed on to function \code{FUN}. Usually, the +output of a call to \code{\link{create_models}} in combination with \code{\link{add_priors}}.} + +\item{...}{further arguments passed to or from other methods.} + +\item{FUN}{the function to be applied to each country model in argument \code{object}. +If \code{NULL} (default), the internal functions \code{\link{bvecxpost}} is used.} + +\item{mc.cores}{the number of cores to use, i.e. at most how many child +processes will be run simultaneously. The option is initialized from +environment variable MC_CORES if set. Must be at least one, and +parallelization requires at least two cores.} + +\item{ctry}{character vector specifying for which countries posterior distributions +should be drawn. If \code{NULL} (default), draws are generated for all country models.} +} +\value{ +An object of class \code{bgvecest}, which contains a list of data, +model specifications, priors, coefficient draws and information +criteria for each estimated country model. +} +\description{ +Estimates the country models of a Bayesian GVAR model. +} diff --git a/man/plot.bgvecest.Rd b/man/plot.bgvecest.Rd new file mode 100644 index 0000000..de00856 --- /dev/null +++ b/man/plot.bgvecest.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot.bgvecest.R +\name{plot.bgvecest} +\alias{plot.bgvecest} +\title{Plotting Draws of a VECX Submodel of a GVAR Model} +\usage{ +\method{plot}{bgvecest}(x, ci = 0.95, type = "hist", ctry = NULL, ...) +} +\arguments{ +\item{x}{an object of class \code{"bgvecest"}, usually, a result of a call to \code{\link{draw_posterior}}.} + +\item{ci}{interval used to calculate credible bands for time-varying parameters.} + +\item{type}{either \code{"hist"} (default) for histograms, \code{"trace"} for a trace plot, +or \code{"boxplot"} for a boxplot. Only used for parameter draws of constant coefficients.} + +\item{ctry}{character. Name of the element in argument \code{x}, for which posterior draws +should be plotted. If \code{NULL} (default), all submodels are used.} + +\item{...}{further graphical parameters.} +} +\description{ +A plot function for objects of class \code{"bgvecest"}. +} diff --git a/man/plot.ctryvecest.Rd b/man/plot.ctryvecest.Rd new file mode 100644 index 0000000..4cb0fcb --- /dev/null +++ b/man/plot.ctryvecest.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot.ctryvecest.R +\name{plot.ctryvecest} +\alias{plot.ctryvecest} +\title{Plotting Draws of a VECX Submodel of a GVAR Model} +\usage{ +\method{plot}{ctryvecest}(x, ci = 0.95, type = "hist", variables = "all", ctry = NULL, ...) +} +\arguments{ +\item{x}{an object of class \code{"ctryvecest"}, usually, a result of a call to \code{\link{draw_posterior}}.} + +\item{ci}{interval used to calculate credible bands for time-varying parameters.} + +\item{type}{either \code{"hist"} (default) for histograms, \code{"trace"} for a trace plot, +or \code{"boxplot"} for a boxplot. Only used for parameter draws of constant coefficients.} + +\item{variables}{character vector of variables that should be plotted. Default is \code{"all"}. +Other options are \code{"domestic"}, \code{"foreign"}, \code{"global"}, \code{"deterministic"} +and \code{"sigma"}.} + +\item{ctry}{character (optional). Name of the country, which will be shown in the title of the output.} + +\item{...}{further graphical parameters.} +} +\description{ +A plot function for objects of class \code{"ctryvecest"} for visual inspection +of posterior draws. +} diff --git a/man/summary.bgvecest.Rd b/man/summary.bgvecest.Rd new file mode 100644 index 0000000..8b669be --- /dev/null +++ b/man/summary.bgvecest.Rd @@ -0,0 +1,40 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/summary.bgvecest.R, R/print.summary.bgvecest.R +\name{summary.bgvecest} +\alias{summary.bgvecest} +\alias{print.summary.bgvecest} +\title{Summarising Country-Specific VECX Submodels of a GVAR Model} +\usage{ +\method{summary}{bgvecest}(object, ci = 0.95, period = NULL, ctry = NULL, ...) + +\method{print}{summary.bgvecest}(x, digits = max(3L, getOption("digits") - 3L), ...) +} +\arguments{ +\item{object}{an object of class \code{"bgvecest"}, usually, a result of a call to +\code{\link{draw_posterior}}.} + +\item{ci}{a numeric between 0 and 1 specifying the probability of the credible band. +Defaults to 0.95.} + +\item{period}{integer. Index of the period, for which the summary statistics should be generated. +Only used for TVP or SV models. Default is \code{NULL}, so that the posterior draws of the last time period +are used.} + +\item{ctry}{character. Name of the element in argument \code{object}, for which summary +statistics should be obtained. If \code{NULL} (default), all country submodels are used.} + +\item{...}{further arguments passed to or from other methods.} + +\item{x}{an object of class \code{"summary.bgvecest"}, usually, a result of a call to +\code{\link{summary.bgvecest}}.} + +\item{digits}{the number of significant digits to use when printing.} +} +\value{ +\code{summary.bgvecest} returns a list of class \code{"summary.bgvecest"}, +which contains summary statistics of country-specific submodels of class +\code{\link{summary.ctryvecest}}. +} +\description{ +summary method for class \code{"bgvecest"}. +} diff --git a/man/summary.ctryvecest.Rd b/man/summary.ctryvecest.Rd new file mode 100644 index 0000000..b685651 --- /dev/null +++ b/man/summary.ctryvecest.Rd @@ -0,0 +1,41 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/summary.ctryvecest.R, +% R/print.summary.ctryvecest.R +\name{summary.ctryvecest} +\alias{summary.ctryvecest} +\alias{print.summary.ctryvecest} +\title{Summarising Country-Specific VECX Submodels of a GVAR Model} +\usage{ +\method{summary}{ctryvecest}(object, ci = 0.95, period = NULL, ...) + +\method{print}{summary.ctryvecest}(x, digits = max(3L, getOption("digits") - 3L), ...) +} +\arguments{ +\item{object}{an object of class \code{"ctryvecest"}, usually, a result of a call to \code{\link{draw_posterior}}.} + +\item{ci}{a numeric between 0 and 1 specifying the probability of the credible band. +Defaults to 0.95.} + +\item{period}{integer. Index of the period, for which the summary statistics should be generated. +Only used for TVP or SV models. Default is \code{NULL}, so that the posterior draws of the last time period +are used.} + +\item{...}{further arguments passed to or from other methods.} + +\item{x}{an object of class \code{"summary.ctryvecest"}, usually, a result of a call to +\code{\link{summary.ctryvecest}}.} + +\item{digits}{the number of significant digits to use when printing.} +} +\value{ +\code{summary.ctryvecest} returns a list of class \code{"summary.ctryvecest"}, +which contains the following components: +\item{coefficients}{A list of various summary statistics of the posterior +draws of the VEC coefficients.} +\item{sigma}{A list of various summary statistics of the posterior +draws of the variance-covariance matrix.} +\item{model}{a list containing information on the model specifications.} +} +\description{ +summary method for class \code{"ctryvecest"}. +} diff --git a/man/thin.bgvecest.Rd b/man/thin.bgvecest.Rd new file mode 100644 index 0000000..1e6e676 --- /dev/null +++ b/man/thin.bgvecest.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/thin.bgvecest.R +\name{thin.bgvecest} +\alias{thin.bgvecest} +\title{Thinning Posterior Draws} +\usage{ +\method{thin}{bgvecest}(x, thin = 10, ...) +} +\arguments{ +\item{x}{an object of class \code{"bgvecest"}.} + +\item{thin}{an integer specifying the thinning interval between successive values of posterior draws.} + +\item{...}{further arguments passed to or from other methods.} +} +\value{ +An object of class \code{"bgvecest"}. +} +\description{ +Thins the MCMC posterior draws in an object of class \code{"bgvecest"}. +} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 9c47d36..8d541ad 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -33,6 +33,28 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// bgvecalg +Rcpp::List bgvecalg(Rcpp::List object); +RcppExport SEXP _bgvars_bgvecalg(SEXP objectSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Rcpp::List >::type object(objectSEXP); + rcpp_result_gen = Rcpp::wrap(bgvecalg(object)); + return rcpp_result_gen; +END_RCPP +} +// bgvectvpalg +Rcpp::List bgvectvpalg(Rcpp::List object); +RcppExport SEXP _bgvars_bgvectvpalg(SEXP objectSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< Rcpp::List >::type object(objectSEXP); + rcpp_result_gen = Rcpp::wrap(bgvectvpalg(object)); + return rcpp_result_gen; +END_RCPP +} // draw_forecast arma::mat draw_forecast(int& i, int& k, arma::mat& a0, arma::mat& a, Rcpp::Nullable& b_, Rcpp::Nullable& c_, arma::mat& sigma, arma::mat pred); RcppExport SEXP _bgvars_draw_forecast(SEXP iSEXP, SEXP kSEXP, SEXP a0SEXP, SEXP aSEXP, SEXP b_SEXP, SEXP c_SEXP, SEXP sigmaSEXP, SEXP predSEXP) { @@ -97,6 +119,8 @@ END_RCPP static const R_CallMethodDef CallEntries[] = { {"_bgvars_bgvaralg", (DL_FUNC) &_bgvars_bgvaralg, 1}, {"_bgvars_bgvartvpalg", (DL_FUNC) &_bgvars_bgvartvpalg, 1}, + {"_bgvars_bgvecalg", (DL_FUNC) &_bgvars_bgvecalg, 1}, + {"_bgvars_bgvectvpalg", (DL_FUNC) &_bgvars_bgvectvpalg, 1}, {"_bgvars_draw_forecast", (DL_FUNC) &_bgvars_draw_forecast, 8}, {"_bgvars_gir", (DL_FUNC) &_bgvars_gir, 4}, {"_bgvars_ir", (DL_FUNC) &_bgvars_ir, 5}, diff --git a/src/bgvecalg.cpp b/src/bgvecalg.cpp new file mode 100644 index 0000000..9895b60 --- /dev/null +++ b/src/bgvecalg.cpp @@ -0,0 +1,890 @@ +#include +#include +// [[Rcpp::depends(RcppArmadillo)]] +// [[Rcpp::depends(bvartools)]] + +// [[Rcpp::export(.bgvecalg)]] +Rcpp::List bgvecalg(Rcpp::List object) { + + // Initialise variables + Rcpp::List data = object["data"]; + const arma::mat y = arma::trans(Rcpp::as(data["Y"])); + const arma::mat yvec = arma::vectorise(y); + const arma::mat w = arma::trans(Rcpp::as(data["W"])); + Rcpp::Nullable x_test = data["X"]; + arma::mat x; + if (x_test.isNotNull()) { + x = arma::trans(Rcpp::as(data["X"])); + } + Rcpp::List initial = object["initial"]; + + // Model information + Rcpp::List model = object["model"]; + Rcpp::CharacterVector model_names = model.names(); + Rcpp::List endogen = model["domestic"]; + Rcpp::CharacterVector endo_names = Rcpp::as(endogen["variables"]); + + // Define useful variables + const int tt = y.n_cols; + const int k_dom = y.n_rows; + const int p_dom = Rcpp::as(endogen["lags"]); + int n_a0 = 0; + int n_dom = 0; + if (p_dom > 1) { + n_dom = k_dom * k_dom * (p_dom - 1); + } + int k_for = 0; + int p_for = 0; + int n_for = 0; + int k_glo = 0; + int p_glo = 0; + int n_glo = 0; + int n_x = 0; + if (x_test.isNotNull()) { + n_x = x.n_rows; + } + int n_r = 0; + int n_ur = 0; + int n_c_ur = 0; + int n_psi = 0; + const int n_sigma = k_dom * k_dom; + const int r = Rcpp::as(model["rank"]); + bool use_rr = false; + int n_w = 0; + if (r > 0) { + use_rr = true; + n_w = w.n_rows; + } + const int n_alpha = k_dom * r; + const int n_beta = n_w * r; + + const bool sv = Rcpp::as(model["sv"]); + const bool structural = Rcpp::as(model["structural"]); + if (structural) { + n_a0 = k_dom * (k_dom - 1) / 2; + } + + bool covar = false; + bool varsel = false; + bool psi_varsel = false; + bool ssvs = false; + bool psi_ssvs = false; + bool bvs = false; + bool psi_bvs = false; + + const int n_nonalpha = k_dom * n_x + n_a0; + const int n_tot = n_alpha + n_nonalpha; + const arma::mat diag_k = arma::eye(k_dom, k_dom); + const arma::mat diag_tot = arma::eye(n_tot, n_tot); + const arma::mat diag_tt = arma::eye(tt, tt); + arma::mat diag_r; + arma::mat diag_beta; + arma::mat z = arma::zeros(tt * k_dom, n_tot); + arma::mat BB_sqrt, z_beta, alpha, Alpha, beta, Beta; + arma::vec y_beta; + Rcpp::List init_coint; + if (use_rr) { + init_coint = initial["cointegration"]; + alpha = arma::zeros(k_dom, r); + beta = Rcpp::as(init_coint["beta"]); + diag_r = arma::eye(r, r); + diag_beta = arma::eye(n_beta, n_beta); + z.cols(0, n_alpha - 1) = arma::kron(arma::trans(arma::trans(beta) * w), diag_k); + if (n_x > 0) { + z.cols(n_alpha, n_tot - 1) = Rcpp::as(data["SUR"]).cols(k_dom * n_w, k_dom * (n_w + n_x) + n_a0 - 1); + } + z_beta = arma::zeros(k_dom * tt, n_beta); + } else { + z = Rcpp::as(data["SUR"]).cols(k_dom * w.n_rows, k_dom * (w.n_rows + n_x) + n_a0 - 1); + } + + // Foreign variables + Rcpp::CharacterVector foreign_names; + Rcpp::List foreign = model["foreign"]; + foreign_names = Rcpp::as(foreign["variables"]); + k_for = foreign_names.length(); + p_for = Rcpp::as(foreign["lags"]); + n_for = k_dom * k_for * p_for; + + // Global variables + Rcpp::CharacterVector global_names; + Rcpp::List global; + if (std::find(model_names.begin(), model_names.end(), "global") != model_names.end()) { + global = model["global"]; + global_names = Rcpp::as(global["variables"]); + k_glo = global_names.length(); + p_glo = Rcpp::as(global["lags"]); + n_glo = k_dom * k_glo * p_glo; + } + + // Deterministic terms + Rcpp::CharacterVector determ_names, det_names_r, det_names_ur; + Rcpp::List determ; + if (std::find(model_names.begin(), model_names.end(), "deterministic") != model_names.end()) { + determ = model["deterministic"]; + determ_names = determ.names(); + if (std::find(determ_names.begin(), determ_names.end(), "restricted") != determ_names.end()) { + det_names_r = Rcpp::as(determ["restricted"]); + n_r = det_names_r.length(); + } + if (std::find(determ_names.begin(), determ_names.end(), "unrestricted") != determ_names.end()) { + det_names_ur = Rcpp::as(determ["unrestricted"]); + n_ur = det_names_ur.length(); + n_c_ur = k_dom * n_ur; + } + } + + // Priors ---- + Rcpp::List priors = object["priors"]; + Rcpp::CharacterVector priors_names = priors.names(); + Rcpp::List init_sigma = initial["sigma"]; + + // Priors - Cointegration + Rcpp::List priors_cointegration; + Rcpp::CharacterVector prcoint_names; + double coint_v_i; + arma::mat beta_post_v, g_i, post_beta_mu, p_tau_i; + if (use_rr) { + priors_cointegration = priors["cointegration"]; + prcoint_names = priors_cointegration.names(); + } + + // Priors - Non-cointegration coefficients + Rcpp::List init_noncoint, priors_coefficients; + Rcpp::CharacterVector prcoeff_names; + arma::vec gamma_prior_mu; + arma::mat gamma_prior_vi; + if (std::find(priors_names.begin(), priors_names.end(), "noncointegration") != priors_names.end()) { + init_noncoint = initial["noncointegration"]; + priors_coefficients = priors["noncointegration"]; + prcoeff_names = priors_coefficients.names(); + } else { + if (!use_rr) { + Rcpp::stop("Cointegration rank is zero and no non-cointegration priors provided."); + } + } + + // Put cointegration and non-cointegration together + if (use_rr) { + gamma_prior_mu = arma::zeros(n_tot); + gamma_prior_vi = arma::zeros(n_tot, n_tot); + coint_v_i = Rcpp::as(priors_cointegration["v_i"]); + p_tau_i = Rcpp::as(priors_cointegration["p_tau_i"]); + if (n_x > 0) { + gamma_prior_mu.subvec(n_alpha, n_tot - 1) = Rcpp::as(priors_coefficients["mu"]); + gamma_prior_vi.submat(n_alpha, n_alpha, n_tot - 1, n_tot - 1) = Rcpp::as(priors_coefficients["v_i"]); + } + } else { + gamma_prior_mu = Rcpp::as(priors_coefficients["mu"]); + gamma_prior_vi = Rcpp::as(priors_coefficients["v_i"]); + } + + // Priors - Coefficient - Variables selection + arma::vec a_prior_incl, a_tau0, a_tau1, a_tau0sq, a_tau1sq; + Rcpp::List a_prior_varsel; + + if (std::find(prcoeff_names.begin(), prcoeff_names.end(), "ssvs") != prcoeff_names.end()) { + if (n_tot - n_alpha > 0) { + ssvs = true; + a_prior_varsel = priors_coefficients["ssvs"]; + a_prior_incl = arma::zeros(n_tot, 1); + a_prior_incl.submat(n_alpha, 0, n_tot - 1, 0) = Rcpp::as(a_prior_varsel["inprior"]); + a_tau0 = arma::zeros(n_tot); + a_tau1 = arma::zeros(n_tot); + a_tau0.subvec(n_alpha, n_tot - 1) = Rcpp::as(a_prior_varsel["tau0"]); + a_tau1.subvec(n_alpha, n_tot - 1) = Rcpp::as(a_prior_varsel["tau1"]); + a_tau0sq = arma::square(a_tau0); + a_tau1sq = arma::square(a_tau1); + } else { + Rcpp::stop("Model does not contain non-cointegration variables for SSVS."); + } + } + if (sv && ssvs) { + Rcpp::stop("Not allowed to use SSVS with stochastic volatility."); + } + + arma::vec a_bvs_lprior_0, a_bvs_lprior_1; + if (std::find(prcoeff_names.begin(), prcoeff_names.end(), "bvs") != prcoeff_names.end()) { + if (n_tot - n_alpha > 0) { + bvs = true; + a_prior_varsel = priors_coefficients["bvs"]; + a_bvs_lprior_0 = arma::zeros(n_tot); + a_bvs_lprior_1 = arma::zeros(n_tot); + a_bvs_lprior_0.submat(n_alpha, 0, n_tot - 1, 0) = arma::log(1 - Rcpp::as(a_prior_varsel["inprior"])); + a_bvs_lprior_1.submat(n_alpha, 0, n_tot - 1, 0) = arma::log(Rcpp::as(a_prior_varsel["inprior"])); + } else { + Rcpp::stop("Model does not contain non-cointegration variables for BVS."); + } + } + + varsel = ssvs || bvs; + + arma::mat gamma_post_mu = gamma_prior_mu * 0; + arma::mat gamma_post_v = gamma_prior_vi * 0; + arma::mat gamma; + int a_varsel_n, a_varsel_pos; + double a_lambda_draw, a_l0, a_l1, a_bayes, a_bayes_rand; + arma::vec post_a_incl, a_u0, a_u1, a_theta0_res, a_theta1_res, a_randu, a_varsel_include, a_varsel_include_draw;; + arma::mat a_AG, a_lambda, a_theta0, a_theta1, z_bvs; + if (varsel) { + a_varsel_include = Rcpp::as(a_prior_varsel["include"]) - 1 + n_alpha; + a_varsel_n = size(a_varsel_include)(0); + if (ssvs) { + a_lambda = arma::ones(n_tot, 1); + } + if (bvs) { + a_lambda = arma::eye(n_tot, n_tot); + a_l0 = 0; + a_l1 = 0; + a_bayes = 0; + a_bayes_rand = 0; + z_bvs = z; + } + } + + /////////////////////////////////////////////////////////////////////// + // Priors & initial values - Error variances + /////////////////////////////////////////////////////////////////////// + + // Empty objects + Rcpp::List sigma_pr = priors["sigma"]; + Rcpp::CharacterVector sigma_names = sigma_pr.names(); + bool use_gamma = false; + double sigma_post_df; + arma::vec h_const, h_init, h_init_post_mu, sigma_h, u_vec, sigma_post_scale, + sigma_post_shape, sigma_prior_mu, sigma_prior_rate; + arma::mat diag_omega_i, diag_sigma_i, diag_sigma_i_temp, h, h_init_post_v, + h_lag, sigma_h_i, sigma_i, sigma_prior_scale, sigma_prior_vi, sse, + omega_i, u; + u = y * 0; + diag_sigma_i = arma::zeros(k_dom * tt, k_dom * tt); + + if (sv) { + // Initial conditions of state equation + sigma_prior_mu = Rcpp::as(sigma_pr["mu"]); + sigma_prior_vi = Rcpp::as(sigma_pr["v_i"]); + // Variances of state equation + sigma_post_shape = Rcpp::as(sigma_pr["shape"]) + 0.5 * tt; + sigma_prior_rate = Rcpp::as(sigma_pr["rate"]); + + h = Rcpp::as(init_sigma["h"]); + h_lag = h * 0; + sigma_h = Rcpp::as(init_sigma["sigma_h"]); + h_init = arma::vectorise(h.row(0)); + sigma_i = arma::diagmat(1 / exp(h_init)); + h_const = Rcpp::as(init_sigma["constant"]); + } else { + if (std::find(sigma_names.begin(), sigma_names.end(), "df") != sigma_names.end()) { + sigma_post_df = Rcpp::as(sigma_pr["df"]) + tt; + sigma_prior_scale = Rcpp::as(sigma_pr["scale"]); + } + if (std::find(sigma_names.begin(), sigma_names.end(), "shape") != sigma_names.end()) { + use_gamma = true; + sigma_post_shape = Rcpp::as(sigma_pr["shape"]) + 0.5 * tt; + sigma_prior_rate = Rcpp::as(sigma_pr["rate"]); + } + + omega_i = Rcpp::as(init_sigma["sigma_i"]); + sigma_i = omega_i; + } + diag_sigma_i.diag() = arma::repmat(sigma_i.diag(), tt, 1); + if (covar || sv) { + diag_omega_i = diag_sigma_i; + } + g_i = sigma_i; + + /////////////////////////////////////////////////////////////////////// + // Covariances + /////////////////////////////////////////////////////////////////////// + + Rcpp::List psi_priors, psi_prior_varsel; + Rcpp::CharacterVector psi_priors_names; + arma::vec psi_prior_incl, psi_tau0, psi_tau1, psi_tau0sq, psi_tau1sq, psi_bvs_lprior_0, psi_bvs_lprior_1; + arma::mat psi_prior_mu, psi_prior_vi; + int psi_varsel_n, psi_varsel_pos; + double psi_bayes, psi_bayes_rand, psi_l0, psi_l1, psi_lambda_draw; + arma::vec psi_post_incl, psi_post_mu, psi_randu, psi_theta0_res, psi_theta1_res, psi_varsel_include, psi_varsel_include_draw, psi_u0, psi_u1, psi_y; + arma::mat diag_covar_omega_i, diag_Psi, psi, Psi, psi_AG, psi_lambda, psi_post_v, psi_theta0, psi_theta1, psi_z, psi_z_bvs; + if (std::find(priors_names.begin(), priors_names.end(), "psi") != priors_names.end()) { + covar = true; + psi_priors = priors["psi"]; + psi_priors_names = psi_priors.names(); + psi_prior_mu = Rcpp::as(psi_priors["mu"]); + psi_prior_vi = Rcpp::as(psi_priors["v_i"]); + + if (std::find(psi_priors_names.begin(), psi_priors_names.end(), "ssvs") != psi_priors_names.end()) { + psi_ssvs = true; + psi_prior_varsel = psi_priors["ssvs"]; + psi_prior_incl = Rcpp::as(psi_prior_varsel["inprior"]); + psi_tau0 = Rcpp::as(psi_prior_varsel["tau0"]); + psi_tau1 = Rcpp::as(psi_prior_varsel["tau1"]); + psi_tau0sq = arma::square(psi_tau0); + psi_tau1sq = arma::square(psi_tau1); + } + if (sv && psi_ssvs) { + Rcpp::stop("Not allowed to use SSVS with stochastic volatility."); + } + + if (std::find(psi_priors_names.begin(), psi_priors_names.end(), "bvs") != psi_priors_names.end()) { + psi_bvs = true; + psi_prior_varsel = psi_priors["bvs"]; + psi_bvs_lprior_0 = arma::log(1 - Rcpp::as(psi_prior_varsel["inprior"])); + psi_bvs_lprior_1 = arma::log(Rcpp::as(psi_prior_varsel["inprior"])); + } + + psi_varsel = psi_ssvs || psi_bvs; + + n_psi = k_dom * (k_dom - 1) / 2; + Psi = arma::eye(k_dom, k_dom); + psi_z = arma::zeros((k_dom - 1) * tt, n_psi); + if (psi_varsel) { + psi_varsel_include = Rcpp::as(psi_prior_varsel["include"]) - 1; + psi_varsel_n = size(psi_varsel_include)(0); + if (psi_ssvs) { + psi_lambda = arma::ones(n_psi, 1); + } + if (psi_bvs) { + psi_lambda = arma::eye(n_psi, n_psi); + psi_l0 = 0; + psi_l1 = 0; + psi_bayes = 0; + psi_bayes_rand = 0; + } + } + diag_covar_omega_i = arma::zeros(tt * (k_dom - 1), tt * (k_dom - 1)); + } + + + // Storage objects + const int iter = Rcpp::as(model["iterations"]); + const int burnin = Rcpp::as(model["burnin"]); + const int draws = iter + burnin; + int pos_draw; + const int alpha_pos_start = 0; + const int alpha_pos_end = n_alpha - 1; + const int dom_pos_start = n_alpha; + const int dom_pos_end = n_alpha + n_dom - 1; + const int for_pos_start = n_alpha + n_dom; + const int for_pos_end = n_alpha + n_dom + n_for - 1; + const int glo_pos_start = n_alpha + n_dom + n_for; + const int glo_pos_end = n_alpha + n_dom + n_for + n_glo - 1; + const int c_pos_start = n_alpha + n_dom + n_for + n_glo; + const int c_pos_end = n_alpha + n_dom + n_for + n_glo + n_c_ur - 1; + const int a0_pos_start = n_alpha + n_dom + n_for + n_glo + n_c_ur; + const int a0_pos_end = n_alpha + n_dom + n_for + n_glo + n_c_ur + n_a0 - 1; + + arma::mat draws_alpha, draws_beta; + if (use_rr) { + draws_alpha = arma::zeros(n_alpha, iter); + draws_beta = arma::zeros(n_beta, iter); + } + arma::mat draws_a0 = arma::zeros(n_a0, iter); + arma::mat draws_dom = arma::zeros(n_dom, iter); + arma::mat draws_for = arma::zeros(n_for, iter); + arma::mat draws_glo = arma::zeros(n_glo, iter); + arma::mat draws_c = arma::zeros(n_c_ur, iter); + arma::mat draws_sigma, draws_sigma_sigma; + if (sv) { + draws_sigma = arma::zeros(k_dom * k_dom * tt, iter); + draws_sigma_sigma = arma::zeros(k_dom * k_dom, iter); + } else { + draws_sigma = arma::zeros(k_dom * k_dom, iter); + } + + arma::vec lambda_vec, psi_lambda_vec; + arma::mat draws_lambda_a0, draws_lambda_dom, draws_lambda_for, draws_lambda_glo, draws_lambda_c; + if (varsel) { + if (structural) { + draws_lambda_a0 = arma::zeros(n_a0, iter); + } + draws_lambda_dom = arma::zeros(n_dom, iter); + draws_lambda_for = arma::zeros(n_for, iter); + draws_lambda_glo = arma::zeros(n_glo, iter); + draws_lambda_c = arma::zeros(n_c_ur, iter); + } + if (covar && psi_varsel) { + draws_lambda_a0 = arma::zeros(n_psi, iter); + } + + // Start Gibbs sampler + for (int draw = 0; draw < draws; draw++) { + + if (draw % 20 == 0) { // Check for user interuption ever now and then + Rcpp::checkUserInterrupt(); + } + + // Draw non-cointegration coefficients ---- + + if (use_rr) { // Update priors for alpha + gamma_prior_vi.submat(0, 0, n_alpha - 1, n_alpha - 1) = arma::kron(coint_v_i * (arma::trans(beta) * p_tau_i * beta), g_i); + if (bvs) { + z_bvs.cols(0, n_alpha - 1) = arma::kron(arma::trans(arma::trans(beta) * w), diag_k); + } else { + z.cols(0, n_alpha - 1) = arma::kron(arma::trans(arma::trans(beta) * w), diag_k); + } + } + if (bvs) { + z = z_bvs * a_lambda; + } + gamma_post_v = gamma_prior_vi + arma::trans(z) * diag_sigma_i * z; + gamma_post_mu = arma::solve(gamma_post_v, gamma_prior_vi * gamma_prior_mu + arma::trans(z) * diag_sigma_i * yvec); + gamma = gamma_post_mu + arma::solve(arma::chol(gamma_post_v), arma::randn(n_tot)); + + // Variables selection + if (varsel) { + + // Reorder positions of variable selection + a_varsel_include_draw = shuffle(a_varsel_include); + + if (ssvs) { + // Obtain inclusion posterior + a_u0 = 1 / a_tau0 % arma::exp(-(arma::square(gamma) / (2 * a_tau0sq))) % (1 - a_prior_incl); + a_u1 = 1 / a_tau1 % arma::exp(-(arma::square(gamma) / (2 * a_tau1sq))) % a_prior_incl; + post_a_incl = a_u1 / (a_u0 + a_u1); + + // Draw inclusion parameters in random order + for (int i = 0; i < a_varsel_n; i++){ + a_lambda_draw = Rcpp::as(Rcpp::rbinom(1, 1, post_a_incl(a_varsel_include_draw(i)))); + a_lambda(a_varsel_include_draw(i), 0) = a_lambda_draw; + if (a_lambda_draw == 0) { + gamma_prior_vi(a_varsel_include_draw(i), a_varsel_include_draw(i)) = 1 / a_tau0sq(a_varsel_include_draw(i)); + } else { + gamma_prior_vi(a_varsel_include_draw(i), a_varsel_include_draw(i)) = 1 / a_tau1sq(a_varsel_include_draw(i)); + } + } + lambda_vec = arma::vectorise(a_lambda); + } + + if (bvs) { + z = z_bvs; + a_AG = a_lambda * gamma; + for (int j = 0; j < a_varsel_n; j++){ + a_varsel_pos = a_varsel_include_draw(j); + a_randu = arma::log(arma::randu(1)); + if (a_lambda(a_varsel_pos, a_varsel_pos) == 1 && a_randu(0) >= a_bvs_lprior_1(a_varsel_pos)){continue;} + if (a_lambda(a_varsel_pos, a_varsel_pos) == 0 && a_randu(0) >= a_bvs_lprior_0(a_varsel_pos)){continue;} + if ((a_lambda(a_varsel_pos, a_varsel_pos) == 1 && a_randu(0) < a_bvs_lprior_1(a_varsel_pos)) || (a_lambda(a_varsel_pos, a_varsel_pos) == 0 && a_randu(0) < a_bvs_lprior_0(a_varsel_pos))){ + a_theta0 = a_AG; + a_theta1 = a_AG; + a_theta0.row(a_varsel_pos) = 0; + a_theta1.row(a_varsel_pos) = gamma.row(a_varsel_pos); + a_theta0_res = yvec - z * a_theta0; + a_theta1_res = yvec - z * a_theta1; + a_l0 = -arma::as_scalar(trans(a_theta0_res) * diag_sigma_i * a_theta0_res) / 2 + arma::as_scalar(a_bvs_lprior_0(a_varsel_pos)); + a_l1 = -arma::as_scalar(trans(a_theta1_res) * diag_sigma_i * a_theta1_res) / 2 + arma::as_scalar(a_bvs_lprior_1(a_varsel_pos)); + a_bayes = a_l1 - a_l0; + a_bayes_rand = arma::as_scalar(arma::log(arma::randu(1))); + if (a_bayes >= a_bayes_rand){ + a_lambda(a_varsel_pos, a_varsel_pos) = 1; + } else { + a_lambda(a_varsel_pos, a_varsel_pos) = 0; + } + } + } + gamma = a_lambda * gamma; + lambda_vec = a_lambda.diag(); + } + } + + if (n_x > 0) { + y_beta = yvec - z.cols(n_alpha, n_tot - 1) * gamma.rows(n_alpha, n_tot - 1); + } else { + y_beta = yvec; + } + + // Cointegration + if (use_rr) { + // Reparameterise alpha + alpha = arma::reshape(gamma.rows(0, n_alpha - 1), k_dom, r); + Alpha = alpha * arma::solve(arma::sqrtmat_sympd(alpha.t() * alpha), diag_r); + + // Draw Beta + for (int i = 0; i < tt; i++){ + z_beta.rows(i * k_dom, (i + 1) * k_dom - 1) = arma::kron(Alpha, arma::trans(w.col(i))); + } + beta_post_v = arma::kron(Alpha.t() * g_i * Alpha, coint_v_i * p_tau_i) + arma::trans(z_beta) * diag_sigma_i * z_beta; + post_beta_mu = arma::solve(beta_post_v, arma::trans(z_beta) * diag_sigma_i * y_beta); + Beta = arma::reshape(post_beta_mu + arma::solve(arma::chol(beta_post_v), arma::randn(n_beta)), n_w, r); + + // Final cointegration values + BB_sqrt = arma::sqrtmat_sympd(arma::trans(Beta) * Beta); + alpha = Alpha * BB_sqrt; + beta = Beta * arma::solve(BB_sqrt, diag_r); + + u_vec = y_beta - arma::vectorise(alpha * beta.t() * w); + } else { + u_vec = y_beta; + } + + u = arma::reshape(u_vec, k_dom, tt); + + /////////////////////////////////////////////////////////////////////// + // Draw covariances + /////////////////////////////////////////////////////////////////////// + if (covar) { + + // Prepare data + psi_y = arma::vectorise(u.rows(1, k_dom - 1)); + for (int i = 1; i < k_dom; i++) { + for (int j = 0; j < tt; j++) { + psi_z.submat(j * (k_dom - 1) + i - 1, + i * (i - 1) / 2, + j * (k_dom - 1) + i - 1, + (i + 1) * i / 2 - 1) = -arma::trans(u.submat(0, j, i - 1, j)); + + diag_covar_omega_i(j * (k_dom - 1) + i - 1, j * (k_dom - 1) + i - 1) = diag_omega_i(j * k_dom + i, j * k_dom + i); + } + } + + if (psi_bvs) { + psi_z_bvs = psi_z; + psi_z = psi_z_bvs * psi_lambda; + } + psi_post_v = psi_prior_vi + arma::trans(psi_z) * diag_covar_omega_i * psi_z; + psi_post_mu = arma::solve(psi_post_v, psi_prior_vi * psi_prior_mu + arma::trans(psi_z) * diag_covar_omega_i * psi_y); + psi = psi_post_mu + arma::solve(arma::chol(psi_post_v), arma::randn(n_psi)); + + if (psi_varsel) { + + // Reorder positions of variable selection + psi_varsel_include_draw = shuffle(psi_varsel_include); + + if (psi_ssvs) { + // Obtain inclusion posterior + psi_u0 = 1 / psi_tau0 % arma::exp(-(arma::square(psi) / (2 * psi_tau0sq))) % (1 - psi_prior_incl); + psi_u1 = 1 / psi_tau1 % arma::exp(-(arma::square(psi) / (2 * psi_tau1sq))) % psi_prior_incl; + psi_post_incl = psi_u1 / (psi_u0 + psi_u1); + + // Draw inclusion parameters in random order + for (int i = 0; i < psi_varsel_n; i++){ + psi_lambda_draw = Rcpp::as(Rcpp::rbinom(1, 1, psi_post_incl(psi_varsel_include_draw(i)))); + psi_lambda(psi_varsel_include_draw(i), 0) = psi_lambda_draw; + if (psi_lambda_draw == 0) { + psi_prior_vi(psi_varsel_include_draw(i), psi_varsel_include_draw(i)) = 1 / psi_tau0sq(psi_varsel_include_draw(i)); + } else { + psi_prior_vi(psi_varsel_include_draw(i), psi_varsel_include_draw(i)) = 1 / psi_tau1sq(psi_varsel_include_draw(i)); + } + } + psi_lambda_vec = arma::vectorise(psi_lambda); + } + + if (psi_bvs) { + psi_z = psi_z_bvs; + psi_AG = psi_lambda * psi; + for (int j = 0; j < psi_varsel_n; j++){ + psi_varsel_pos = psi_varsel_include_draw(j); + psi_randu = arma::log(arma::randu(1)); + if (psi_lambda(psi_varsel_pos, psi_varsel_pos) == 1 && psi_randu(0) >= psi_bvs_lprior_1(psi_varsel_pos)){continue;} + if (psi_lambda(psi_varsel_pos, psi_varsel_pos) == 0 && psi_randu(0) >= psi_bvs_lprior_0(psi_varsel_pos)){continue;} + if ((psi_lambda(psi_varsel_pos, psi_varsel_pos) == 1 && psi_randu(0) < psi_bvs_lprior_1(psi_varsel_pos)) || (psi_lambda(psi_varsel_pos, psi_varsel_pos) == 0 && psi_randu(0) < psi_bvs_lprior_0(psi_varsel_pos))){ + psi_theta0 = psi_AG; + psi_theta1 = psi_AG; + psi_theta0.row(psi_varsel_pos) = 0; + psi_theta1.row(psi_varsel_pos) = psi.row(psi_varsel_pos); + psi_theta0_res = psi_y - psi_z * psi_theta0; + psi_theta1_res = psi_y - psi_z * psi_theta1; + psi_l0 = -arma::as_scalar(trans(psi_theta0_res) * diag_covar_omega_i * psi_theta0_res) / 2 + arma::as_scalar(psi_bvs_lprior_0(psi_varsel_pos)); + psi_l1 = -arma::as_scalar(trans(psi_theta1_res) * diag_covar_omega_i * psi_theta1_res) / 2 + arma::as_scalar(psi_bvs_lprior_1(psi_varsel_pos)); + psi_bayes = psi_l1 - psi_l0; + psi_bayes_rand = arma::as_scalar(arma::log(arma::randu(1))); + if (psi_bayes >= psi_bayes_rand){ + psi_lambda(psi_varsel_pos, psi_varsel_pos) = 1; + } else { + psi_lambda(psi_varsel_pos, psi_varsel_pos) = 0; + } + } + } + psi = psi_lambda * psi; + psi_lambda_vec = psi_lambda.diag(); + } + } + + for (int i = 1; i < k_dom; i++) { + Psi.submat(i, 0, i, i - 1) = arma::trans(psi.submat(i * (i - 1) / 2, 0, (i + 1) * i / 2 - 1, 0)); + } + u = Psi * u; + } + + /////////////////////////////////////////////////////////////////////// + // Draw error variances + /////////////////////////////////////////////////////////////////////// + if (sv) { + + // Draw variances + for (int i = 0; i < k_dom; i++) { + h.col(i) = bvartools::stoch_vol(u.row(i).t(), h.col(i), sigma_h(i), h_init(i), h_const(i)); + } + 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; + } + + // Draw sigma_h + h_lag.row(0) = h_init.t(); + h_lag.rows(1, tt - 1) = h.rows(0, tt - 2); + h_lag = h - h_lag; + sigma_post_scale = 1 / (sigma_prior_rate + arma::trans(arma::sum(arma::pow(h_lag, 2))) * 0.5); + for (int i = 0; i < k_dom; 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_dom)); + + } else { + + sse = u * u.t(); + + if (use_gamma) { + for (int i = 0; i < k_dom; 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))); + } + 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 { + sigma_i = arma::wishrnd(arma::solve(sigma_prior_scale + sse, diag_k), sigma_post_df); + } + + diag_sigma_i = arma::kron(diag_tt, sigma_i); + + } + + // Update g_i + if (use_rr) { + g_i = diag_sigma_i.submat(0, 0, k_dom - 1, k_dom - 1); + if (sv) { + for (int i = 1; i < tt; i++) { + g_i = g_i + diag_sigma_i.submat(i * k_dom, i * k_dom, (i + 1) * k_dom - 1, (i + 1) * k_dom - 1); + } + g_i = g_i / tt; + } + } + + /////////////////////////////////////////////////////////////////////// + // Store draws + /////////////////////////////////////////////////////////////////////// + if (draw >= burnin) { + + pos_draw = draw - burnin; + + if (sv) { + for (int i = 0; i < tt; i ++) { + draws_sigma.submat(i * n_sigma, pos_draw, (i + 1) * n_sigma - 1, pos_draw) = arma::vectorise(arma::solve(diag_sigma_i.submat(i * k_dom, i * k_dom, (i + 1) * k_dom - 1, (i + 1) * k_dom - 1), diag_k)); + } + draws_sigma_sigma.col(pos_draw) = arma::vectorise(arma::diagmat(sigma_h)); + } else { + draws_sigma.col(pos_draw) = arma::vectorise(arma::solve(sigma_i, diag_k)); + } + + if (psi_varsel) { + draws_lambda_a0.col(pos_draw) = psi_lambda_vec; + } + + if (use_rr) { + draws_alpha.col(pos_draw) = arma::vectorise(gamma.rows(alpha_pos_start, alpha_pos_end)); + draws_beta.col(pos_draw) = arma::vectorise(beta.t()); + } + + if (n_dom > 0) { + draws_dom.col(pos_draw) = arma::vectorise(gamma.rows(dom_pos_start, dom_pos_end)); + if (varsel) { + draws_lambda_dom.col(pos_draw) = lambda_vec.subvec(dom_pos_start, dom_pos_end); + } + } + if (n_for > 0) { + draws_for.col(pos_draw) = arma::vectorise(gamma.rows(for_pos_start, for_pos_end)); + if (varsel) { + draws_lambda_for.col(pos_draw) = lambda_vec.subvec(for_pos_start, for_pos_end); + } + } + if (n_glo > 0) { + draws_glo.col(pos_draw) = arma::vectorise(gamma.rows(glo_pos_start, glo_pos_end)); + if (varsel) { + draws_lambda_glo.col(pos_draw) = lambda_vec.subvec(glo_pos_start, glo_pos_end); + } + } + if (n_c_ur > 0) { + draws_c.col(pos_draw) = arma::vectorise(gamma.rows(c_pos_start, c_pos_end)); + if (varsel) { + draws_lambda_c.col(pos_draw) = lambda_vec.subvec(c_pos_start, c_pos_end); + } + } + if (structural) { + draws_a0.col(pos_draw) = arma::vectorise(gamma.rows(a0_pos_start, a0_pos_end)); + if (varsel) { + draws_lambda_a0.col(pos_draw) = lambda_vec.subvec(a0_pos_start, a0_pos_end); + } + } + } + + } // End loop + + Rcpp::List posteriors = Rcpp::List::create(Rcpp::Named("a0") = R_NilValue, + Rcpp::Named("alpha") = R_NilValue, + Rcpp::Named("beta_dom") = R_NilValue, + Rcpp::Named("beta_for") = R_NilValue, + Rcpp::Named("beta_glo") = R_NilValue, + Rcpp::Named("beta_det") = R_NilValue, + Rcpp::Named("gamma_dom") = R_NilValue, + Rcpp::Named("gamma_for") = R_NilValue, + Rcpp::Named("gamma_glo") = R_NilValue, + Rcpp::Named("gamma_det") = R_NilValue, + Rcpp::Named("sigma") = R_NilValue); + + if (use_rr) { + posteriors["alpha"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_alpha)); + + // Reformat draws + for (int i = 0; i < iter; i ++) { + draws_beta.submat(0, i, r * k_dom - 1, i) = arma::vectorise(arma::trans(arma::reshape(draws_beta.submat(0, i, r * k_dom - 1, i), r, k_dom))); + draws_beta.submat(r * k_dom, i, r * (k_dom + k_for) - 1, i) = arma::vectorise(arma::trans(arma::reshape(draws_beta.submat(r * k_dom, i, r * (k_dom + k_for) - 1, i), r, k_for))); + if (k_glo > 0) { + draws_beta.submat(r * (k_dom + k_for), i, r * (k_dom + k_for + k_glo) - 1, i) = arma::vectorise(arma::trans(arma::reshape(draws_beta.submat(r * (k_dom + k_for), i, r * (k_dom + k_for + k_glo) - 1, i), r, k_glo))); + } + if (n_r > 0) { + draws_beta.submat(r * (k_dom + k_for + k_glo), i, r * (k_dom + k_for + k_glo + n_r) - 1, i) = arma::vectorise(arma::trans(arma::reshape(draws_beta.submat(r * (k_dom + k_for + k_glo), i, r * (k_dom + k_for + k_glo + n_r) - 1, i), r, n_r))); + } + } + + posteriors["beta_dom"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_beta.rows(0, r * k_dom - 1))); + if (k_for > 0) { + posteriors["beta_for"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_beta.rows(r * k_dom, r * (k_dom + k_for) - 1))); + } + if (k_glo > 0) { + posteriors["beta_glo"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_beta.rows(r * (k_dom + k_for), r * (k_dom + k_for + k_glo) - 1))); + } + if (n_r > 0) { + posteriors["beta_det"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_beta.rows(r * (k_dom + k_for + k_glo), r * (k_dom + k_for + k_glo + n_r) - 1))); + } + } + + if (n_dom > 0) { + if (varsel) { + posteriors["gamma_dom"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_dom, + Rcpp::Named("lambda") = draws_lambda_dom)); + } else { + posteriors["gamma_dom"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_dom)); + } + } + + if (n_for > 0) { + if (varsel) { + posteriors["gamma_for"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_for, + Rcpp::Named("lambda") = draws_lambda_for)); + } else { + posteriors["gamma_for"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_for)); + } + } + + if (n_glo > 0) { + if (varsel) { + posteriors["gamma_glo"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_glo, + Rcpp::Named("lambda") = draws_lambda_glo)); + } else { + posteriors["gamma_glo"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_glo)); + } + } + + if (n_ur > 0) { + if (varsel) { + posteriors["gamma_det"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_c, + Rcpp::Named("lambda") = draws_lambda_c)); + } else { + posteriors["gamma_det"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_c)); + } + } + + if (structural) { + if (varsel) { + posteriors["a0"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_a0, + Rcpp::Named("lambda") = draws_lambda_a0)); + } else { + posteriors["a0"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_a0)); + } + } + + if (psi_varsel) { + if (sv) { + posteriors["sigma"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_sigma, + Rcpp::Named("sigma") = draws_sigma_sigma, + Rcpp::Named("lambda") = draws_lambda_a0)); + } else { + posteriors["sigma"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_sigma, + Rcpp::Named("lambda") = draws_lambda_a0)); + } + } else { + if (sv) { + posteriors["sigma"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_sigma, + Rcpp::Named("sigma") = draws_sigma_sigma)); + } else { + posteriors["sigma"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_sigma)); + } + } + + return Rcpp::List::create(Rcpp::Named("data") = object["data"], + Rcpp::Named("model") = object["model"], + Rcpp::Named("priors") = object["priors"], + Rcpp::Named("posteriors") = posteriors); + +} + +/*** R + + +#library(bgvars) + +# Load data +data("gvar2019") + +# Create regions +temp <- create_regions(country_data = gvar2019$country_data, + weight_data = gvar2019$weight_data, + region_weights = gvar2019$region_weights, + regions = list(EA = c("AT", "BE", "DE", "ES", "FI", "FR", "IT", "NL")), period = 3) + +country_data <- temp$country_data +weight_data <- temp$weight_data +global_data = gvar2019$global_data + +# Create weights +weight_data <- create_weights(weight_data, period = 3, country_data = country_data) + +# Model ---- + +# Create general model specifications +model_specs <- create_specifications(country_data = country_data, + global_data = global_data, + domestic = list(variables = c("y", "Dp", "r"), lags = 1), + foreign = list(variables = c("y", "Dp", "r"), lags = 1), + global = list(variables = "poil", lags = 1), + deterministic = list("const" = "unrestricted"), + countries = c("EA", "US", "GB", "CA", "JP", "IN"), + type = "VEC", + tvp = FALSE, sv = TRUE, + iterations = 10, + burnin = 10) + +# Create all country models +country_models <- create_models(country_data = country_data, + weight_data = weight_data, + global_data = global_data, + model_specs = model_specs) + +# Add priors +temp_model <- add_priors(country_models, + coef = list(v_i = 1 / 9, v_i_det = 1 / 10), + sigma = list(shape = 3, rate = .0001, mu = 0, v_i = .1, sigma_h = 0.05, constant = .0001, covar = TRUE)) + +.bgvecalg(temp_model[[1]]) + +*/ \ No newline at end of file diff --git a/src/bgvectvpalg.cpp b/src/bgvectvpalg.cpp new file mode 100644 index 0000000..b58e2a1 --- /dev/null +++ b/src/bgvectvpalg.cpp @@ -0,0 +1,985 @@ +#include +#include +// [[Rcpp::depends(RcppArmadillo)]] +// [[Rcpp::depends(bvartools)]] + +// [[Rcpp::export(.bgvectvpalg)]] +Rcpp::List bgvectvpalg(Rcpp::List object) { + + // Initialise variables + Rcpp::List data = object["data"]; + const arma::mat y = arma::trans(Rcpp::as(data["Y"])); + const arma::vec yvec = arma::vectorise(y); + const arma::mat w = arma::trans(Rcpp::as(data["W"])); + Rcpp::Nullable x_test = data["X"]; + arma::mat x; + if (x_test.isNotNull()) { + x = arma::trans(Rcpp::as(data["X"])); + } + + // Model information + Rcpp::List model = object["model"]; + Rcpp::CharacterVector model_names = model.names(); + Rcpp::List endogen = model["domestic"]; + Rcpp::CharacterVector endo_names = Rcpp::as(endogen["variables"]); + + // Define useful variables + const int tt = y.n_cols; + const int k_dom = y.n_rows; + const int p_dom = Rcpp::as(endogen["lags"]); + int n_a0 = 0; + int n_dom = 0; + if (p_dom > 1) { + n_dom = k_dom * k_dom * (p_dom - 1); + } + int k_for = 0; + int p_for = 0; + int n_for = 0; + int k_glo = 0; + int p_glo = 0; + int n_glo = 0; + int k_det_r = 0; + int n_det_ur = 0; + int n_psi = 0; + const int n_sigma = k_dom * k_dom; + const int r = Rcpp::as(model["rank"]); + bool use_rr = false; + int k_w = 0; + int n_w = 0; + if (r > 0) { + use_rr = true; + k_w = w.n_rows; + n_w = w.n_rows * k_dom; + } + const int n_alpha = k_dom * r; + const int n_beta = k_w * r; + + const bool sv = Rcpp::as(model["sv"]); + const bool structural = Rcpp::as(model["structural"]); + if (structural) { + n_a0 = k_dom * (k_dom - 1) / 2; + } + + bool covar = false; + bool bvs = false; + bool psi_bvs = false; + + // Foreign variables + Rcpp::CharacterVector foreign_names; + Rcpp::List foreign = model["foreign"]; + foreign_names = Rcpp::as(foreign["variables"]); + k_for = foreign_names.length(); + p_for = Rcpp::as(foreign["lags"]); + n_for = k_dom * k_for * p_for; + + // Global variables + Rcpp::CharacterVector global_names; + Rcpp::List global; + if (std::find(model_names.begin(), model_names.end(), "global") != model_names.end()) { + global = model["global"]; + global_names = Rcpp::as(global["variables"]); + k_glo = global_names.length(); + p_glo = Rcpp::as(global["lags"]); + n_glo = k_dom * k_glo * p_glo; + } + + // Deterministic terms + Rcpp::CharacterVector determ_names, det_names_r, det_names_ur; + Rcpp::List determ; + if (std::find(model_names.begin(), model_names.end(), "deterministic") != model_names.end()) { + determ = model["deterministic"]; + determ_names = determ.names(); + if (std::find(determ_names.begin(), determ_names.end(), "restricted") != determ_names.end()) { + det_names_r = Rcpp::as(determ["restricted"]); + k_det_r = det_names_r.length(); + } + if (std::find(determ_names.begin(), determ_names.end(), "unrestricted") != determ_names.end()) { + det_names_ur = Rcpp::as(determ["unrestricted"]); + n_det_ur = det_names_ur.length() * k_dom; + } + } + + int n_nonalpha = n_dom + n_for + n_glo + n_det_ur + n_a0; + int n_tot = n_alpha + n_nonalpha; + + // Priors & initial values ---- + Rcpp::List priors = object["priors"]; + Rcpp::CharacterVector priors_names = priors.names(); + Rcpp::List initial = object["initial"]; + + /////////////////////////////////////////////////////////////////////// + // Priors & initial values - Cointegration + /////////////////////////////////////////////////////////////////////// + + Rcpp::List init_coint, priors_cointegration, priors_alpha, priors_beta; + Rcpp::CharacterVector priors_cointegration_names; + //arma::vec priors_rho; + arma::mat beta, beta_b, beta_mu_post, beta_mu_prior, beta_sigma, beta_t, beta_v_post, + beta_vinv_prior, beta_y, beta_z, pi, pi_temp; + arma::vec beta_init; + //arma::mat alpha, beta_help; + double rho; + + + + if (use_rr) { + + priors_cointegration = priors["cointegration"]; + priors_cointegration_names = priors_cointegration.names(); + init_coint = initial["cointegration"]; + + // alpha + priors_alpha = priors_cointegration["alpha"]; + + // beta + priors_beta = priors_cointegration["beta"]; + beta_mu_prior = Rcpp::as(priors_beta["mu"]); + beta_mu_post = beta_mu_prior * 0; + // post_beta_v = prior_beta_vinv * 0; + + // rho (future functionality) + // bool update_rho = false; + // priors_rho = Rcpp::as(priors_cointegration["rho"]); + // if (priors_rho.n_elem == 2) { + // update_rho = true; + // } + + //alpha = arma::zeros(n_alpha, tt); + beta = arma::mat(n_beta, tt); + for (int i = 0; i < tt; i++) { + beta.col(i) = arma::vectorise(Rcpp::as(init_coint["beta"])); + } + beta_y = y * 0; + beta_init = arma::vectorise(Rcpp::as(init_coint["beta"])); + beta_z = arma::zeros(k_dom * tt, n_beta); + beta_sigma = arma::eye(n_beta, n_beta); + rho = Rcpp::as(init_coint["rho"]); + beta_b = arma::eye(n_beta, n_beta) * rho; + pi = arma::zeros(n_w, tt); + + beta_vinv_prior = (1 - rho * rho) * arma::eye(n_beta, n_beta); + } + + /////////////////////////////////////////////////////////////////////// + // Priors & initial values - Non-cointegration + /////////////////////////////////////////////////////////////////////// + + Rcpp::List priors_noncointegration; + Rcpp::CharacterVector priors_noncointegration_names; + if (std::find(priors_names.begin(), priors_names.end(), "noncointegration") != priors_names.end()) { + priors_noncointegration = priors["noncointegration"]; + priors_noncointegration_names = priors_noncointegration.names(); + } else { + if (!use_rr) { + Rcpp::stop("Cointegration rank is zero and no non-cointegration priors provided."); + } + } + + Rcpp::List init_noncoint = initial["noncointegration"]; + + // Priors - BVS + Rcpp::List priors_bvs; + arma::vec gamma_bvs_lprior_0, gamma_bvs_lprior_1, gammma_prior_incl; + if (std::find(priors_noncointegration_names.begin(), priors_noncointegration_names.end(), "bvs") != priors_noncointegration_names.end()) { + if (n_nonalpha > 0) { + bvs = true; + priors_bvs = priors_noncointegration["bvs"]; + gamma_bvs_lprior_0 = arma::log(0.5 * arma::ones(n_tot)); + gamma_bvs_lprior_1 = arma::log(0.5 * arma::ones(n_tot)); + gamma_bvs_lprior_0.subvec(n_alpha, n_tot - 1) = arma::log(1 - Rcpp::as(priors_bvs["inprior"])); + gamma_bvs_lprior_1.subvec(n_alpha, n_tot - 1) = arma::log(Rcpp::as(priors_bvs["inprior"])); + } + } + + // gamma_0 + arma::mat prior_gamma_mu, prior_gamma_vinv, post_gamma_mu, post_gamma_v, post_gamma0_v; + if (use_rr) { + + // Generate empty prior matrices + prior_gamma_mu = arma::zeros(n_tot, 1); + prior_gamma_vinv = arma::zeros(n_tot, n_tot); + + // Add alpha priors + prior_gamma_mu.submat(0, 0, n_alpha - 1, 0) = Rcpp::as(priors_alpha["mu"]); + prior_gamma_vinv.submat(0, 0, n_alpha - 1, n_alpha - 1) = Rcpp::as(priors_alpha["v_i"]); + + // Add non-alpha priors + if (n_nonalpha > 0) { + prior_gamma_mu.submat(n_alpha, 0, n_tot - 1, 0) = Rcpp::as(priors_noncointegration["mu"]); + prior_gamma_vinv.submat(n_alpha, n_alpha, n_tot - 1, n_tot - 1) = Rcpp::as(priors_noncointegration["v_i"]); + } + } else { + // If r = 0, only use non-alpha priors + prior_gamma_mu = Rcpp::as(priors_noncointegration["mu"]); + prior_gamma_vinv = Rcpp::as(priors_noncointegration["v_i"]); + } + arma::vec prior_sigma_v_shape = arma::zeros(n_tot); + arma::vec gamma_sigma_v_prior_rate = arma::zeros(n_tot); + arma::mat gamma_sigma_v = arma::zeros(n_tot, n_tot); + arma::vec gamma_init = arma::zeros(n_tot); + if (use_rr) { + prior_sigma_v_shape.subvec(0, n_alpha - 1) = Rcpp::as(priors_alpha["shape"]); + gamma_sigma_v_prior_rate.subvec(0, n_alpha - 1) = Rcpp::as(priors_alpha["rate"]); + gamma_init.subvec(0, n_alpha - 1) = Rcpp::as(init_coint["alpha"]); + gamma_sigma_v.submat(0, 0, n_alpha - 1, n_alpha - 1) = Rcpp::as(init_coint["sigma_alpha_i"]); + } + if (n_nonalpha > 0) { + prior_sigma_v_shape.subvec(n_alpha, n_tot - 1) = Rcpp::as(priors_noncointegration["shape"]); + gamma_sigma_v_prior_rate.subvec(n_alpha, n_tot - 1) = Rcpp::as(priors_noncointegration["rate"]); + gamma_init.subvec(n_alpha, n_tot - 1) = Rcpp::as(init_noncoint["gamma"]); + gamma_sigma_v.submat(n_alpha, n_alpha, n_tot - 1, n_tot - 1) = Rcpp::as(init_noncoint["sigma_gamma_i"]); + } + gamma_sigma_v.diag() = 1 / gamma_sigma_v.diag(); + arma::vec gamma_sigma_v_post_shape = prior_sigma_v_shape + 0.5 * tt; + arma::vec gamma_sigma_v_post_scale; + arma::mat gamma_sigma_v_i = arma::eye(n_tot, n_tot); + gamma_sigma_v_i.diag() = 1 / gamma_sigma_v_prior_rate; + arma::mat gamma_v; + + arma::vec vec_tt = arma::ones(tt); // T vector + arma::mat diag_k = arma::eye(k_dom, k_dom); // K diag matrix + const arma::mat gamma_b = arma::eye(n_tot, n_tot); + arma::mat gamma = arma::zeros(n_tot, tt); + arma::mat gamma_lag = gamma; + + arma::mat z = arma::zeros(tt * k_dom, n_tot); + for (int i = 0; i < tt; i++) { + if (n_nonalpha > 0) { + z.submat(i * k_dom, n_alpha, (i + 1) * k_dom - 1, n_tot - 1) = Rcpp::as(data["SUR"]).submat(i * k_dom, n_w, (i + 1) * k_dom - 1, n_w + n_nonalpha - 1); + } + } + + // Variable selection + arma::mat gamma_AG, gamma_bvs_l0_res, gamma_bvs_l1_res, gamma_mat, gamma_theta0, gamma_theta1; + arma::vec gamma_randu, gamma_varsel_include, gamma_varsel_include_draw; + arma::mat zz_bvs, gamma_lambda; + int gamma_varsel_n, gamma_varsel_pos; + double gamma_l0, gamma_l1, gamma_bayes, gamma_bayes_rand; + if (bvs) { + gamma_bvs_l0_res = y * 0; + gamma_bvs_l1_res = y * 0; + gamma_varsel_include = Rcpp::as(priors_bvs["include"]) + n_alpha - 1; + gamma_varsel_n = size(gamma_varsel_include)(0); + gamma_lambda = arma::eye(n_tot, n_tot); + gamma_l0 = 0; + gamma_l1 = 0; + gamma_bayes = 0; + gamma_bayes_rand = 0; + zz_bvs = z; + } + + /////////////////////////////////////////////////////////////////////// + // Priors & initial values - Measurement covariances + /////////////////////////////////////////////////////////////////////// + + Rcpp::List init_psi, psi_priors, psi_prior_varsel; + Rcpp::CharacterVector psi_priors_names; + arma::vec psi_prior_incl, psi_tau0, psi_tau1, psi_tau0sq, psi_tau1sq, psi_bvs_lprior_0, psi_bvs_lprior_1; + arma::vec psi_sigma_v_post_scale, psi_sigma_v_post_shape, psi_sigma_v_prior_rate; + arma::mat psi_init_prior_mu, psi_init_prior_vi; + + if (std::find(priors_names.begin(), priors_names.end(), "psi") != priors_names.end()) { + covar = true; + psi_priors = priors["psi"]; + psi_priors_names = psi_priors.names(); + psi_init_prior_mu = Rcpp::as(psi_priors["mu"]); + psi_init_prior_vi = Rcpp::as(psi_priors["v_i"]); + psi_sigma_v_post_shape = Rcpp::as(psi_priors["shape"]) + 0.5 * tt; + psi_sigma_v_prior_rate = Rcpp::as(psi_priors["rate"]); + + if (std::find(psi_priors_names.begin(), psi_priors_names.end(), "bvs") != psi_priors_names.end()) { + psi_bvs = true; + psi_prior_varsel = psi_priors["bvs"]; + psi_bvs_lprior_0 = arma::log(1 - Rcpp::as(psi_prior_varsel["inprior"])); + psi_bvs_lprior_1 = arma::log(Rcpp::as(psi_prior_varsel["inprior"])); + } + } + + // Initial values + int psi_varsel_n, psi_varsel_pos; + double psi_bayes, psi_bayes_rand, psi_l0, psi_l1; + arma::vec psi_init, psi_init_post_mu, psi_post_incl, psi_post_mu, psi_randu, psi_varsel_include, psi_varsel_include_draw, psi_u0, psi_u1; + arma::mat psi, Psi, diag_Psi, psi_AG, psi_b, psi_init_post_v, psi_lag, psi_lambda, psi_mat, psi_sigma_v, psi_sigma_v_i, psi_theta0, psi_theta1, psi_theta0_res, psi_theta1_res, psi_v, psi_y, psi_z, psi_z_bvs; + if (covar) { + n_psi = k_dom * (k_dom - 1) / 2; + Psi = arma::zeros(k_dom * tt, k_dom); + for (int i = 0; i < tt; i++) { + Psi.rows(i * k_dom, (i + 1) * k_dom - 1) = arma::eye(k_dom, k_dom); + } + Rcpp::List init_psi = initial["psi"]; + psi_init = Rcpp::as(init_psi["psi"]); + psi_sigma_v_i = Rcpp::as(init_psi["sigma_psi_i"]); + psi_sigma_v = arma::solve(psi_sigma_v_i, arma::eye(n_psi, n_psi)); + psi_b = arma::eye(n_psi, n_psi); + psi_lag = arma::zeros(n_psi, tt); + psi_z = arma::zeros((k_dom - 1) * tt, n_psi); + if (psi_bvs) { + psi_varsel_include = Rcpp::as(psi_prior_varsel["include"]) - 1; + psi_varsel_n = size(psi_varsel_include)(0); + psi_lambda = arma::eye(n_psi, n_psi); + psi_l0 = 0; + psi_l1 = 0; + psi_theta0_res = arma::zeros(k_dom - 1, tt); + psi_theta1_res = arma::zeros(k_dom - 1, tt); + psi_bayes = 0; + psi_bayes_rand = 0; + } + } + + /////////////////////////////////////////////////////////////////////// + // Priors & initial values - Measurement error variances + /////////////////////////////////////////////////////////////////////// + + Rcpp::List sigma_pr = priors["sigma"]; + Rcpp::CharacterVector sigma_names = sigma_pr.names(); + Rcpp::List init_sigma = initial["sigma"]; + double sigma_post_df; + arma::vec sigma_post_shape, sigma_prior_rate, sigma_prior_mu; + arma::mat sigma_prior_scale, sigma_prior_vi; + bool use_gamma = false; + if (sv) { + sigma_prior_mu = Rcpp::as(sigma_pr["mu"]); + sigma_prior_vi = Rcpp::as(sigma_pr["v_i"]); + sigma_post_shape = Rcpp::as(sigma_pr["shape"]) + 0.5 * tt; + sigma_prior_rate = Rcpp::as(sigma_pr["rate"]); + } else { + if (std::find(sigma_names.begin(), sigma_names.end(), "df") != sigma_names.end()) { + sigma_post_df = Rcpp::as(sigma_pr["df"]) + tt; + sigma_prior_scale = Rcpp::as(sigma_pr["scale"]); + } + if (std::find(sigma_names.begin(), sigma_names.end(), "shape") != sigma_names.end()) { + use_gamma = true; + sigma_post_shape = Rcpp::as(sigma_pr["shape"]) + 0.5 * tt; + sigma_prior_rate = Rcpp::as(sigma_pr["rate"]); + } + } + + // Initial values + arma::vec h_const, h_init, h_init_post_mu, sigma_h, u_vec, sigma_post_scale; + arma::mat h_init_post_v, sigma_h_i, diag_sigma_i_temp, h, h_lag, sigma_u, sigma_u_i, sigma_u_temp, sse, omega_i, omega_psi; + arma::mat u = y * 0; + sigma_u = arma::zeros(k_dom * tt, k_dom); + sigma_u_i = arma::zeros(k_dom * tt, k_dom); + if (sv) { + h = Rcpp::as(init_sigma["h"]); + h_lag = h * 0; + sigma_h = Rcpp::as(init_sigma["sigma_h"]); + h_init = arma::vectorise(h.row(0)); + h_const = Rcpp::as(init_sigma["constant"]); + if (covar) { + omega_psi = arma::zeros((k_dom - 1) * tt, k_dom - 1); + } + for (int i = 0; i < tt; i++) { + sigma_u.rows(i * k_dom, (i + 1) * k_dom - 1) = arma::diagmat(exp(h_init)); + sigma_u_i.rows(i * k_dom, (i + 1) * k_dom - 1) = arma::solve(sigma_u.rows(i * k_dom, (i + 1) * k_dom - 1), diag_k); + if (covar) { + omega_psi.rows(i * (k_dom - 1), (i + 1) * (k_dom - 1) - 1) = sigma_u.submat(i * k_dom + 1, 1, (i + 1) * k_dom - 1, k_dom - 1); + } + } + omega_i = sigma_u_i; + } else { + if (use_gamma) { + omega_i = arma::zeros(k_dom, k_dom); + omega_i.diag() = Rcpp::as(init_sigma["sigma_i"]).diag(); + } else { + omega_i = Rcpp::as(init_sigma["sigma_i"]); + } + for (int i = 0; i < tt; i++) { + sigma_u.rows(i * k_dom, (i + 1) * k_dom - 1) = arma::solve(omega_i, diag_k); + sigma_u_i.rows(i * k_dom, (i + 1) * k_dom - 1) = omega_i; + } + } + + //////////////////////////////////////////// Storage objects /////////////////////////////////////////////////// + + const int iter = Rcpp::as(model["iterations"]); + const int burnin = Rcpp::as(model["burnin"]); + const int draws = iter + burnin; + int pos_draw; + const int alpha_pos_start = 0; + const int alpha_pos_end = n_alpha - 1; + const int dom_pos_start = n_alpha; + const int dom_pos_end = n_alpha + n_dom - 1; + const int for_pos_start = n_alpha + n_dom; + const int for_pos_end = n_alpha + n_dom + n_for - 1; + const int glo_pos_start = n_alpha + n_dom + n_for; + const int glo_pos_end = n_alpha + n_dom + n_for + n_glo - 1; + const int c_pos_start = n_alpha + n_dom + n_for + n_glo; + const int c_pos_end = n_alpha + n_dom + n_for + n_glo + n_det_ur - 1; + const int a0_pos_start = n_alpha + n_dom + n_for + n_glo + n_det_ur; + const int a0_pos_end = n_alpha + n_dom + n_for + n_glo + n_det_ur + n_a0 - 1; + + arma::mat draws_alpha, draws_beta_dom, draws_beta_for, draws_beta_glo, draws_beta_det; + if (use_rr) { + draws_alpha = arma::zeros(n_alpha * tt, iter); + draws_beta_dom = arma::zeros(k_dom * r * tt, iter); + draws_beta_for = arma::zeros(k_for * r * tt, iter); + draws_beta_glo = arma::zeros(k_glo * r * tt, iter); + draws_beta_det = arma::zeros(k_det_r * r * tt, iter); + } + + arma::mat draws_a0 = arma::zeros(n_a0 * tt, iter); + arma::mat draws_a0_sigma = arma::zeros(n_a0, iter); + arma::mat draws_dom = arma::zeros(n_dom * tt, iter); + arma::mat draws_dom_sigma = arma::zeros(n_dom, iter); + arma::mat draws_for = arma::zeros(n_for * tt, iter); + arma::mat draws_for_sigma = arma::zeros(n_for, iter); + arma::mat draws_glo = arma::zeros(n_glo * tt, iter); + arma::mat draws_glo_sigma = arma::zeros(n_glo, iter); + arma::mat draws_det_ur = arma::zeros(n_det_ur * tt, iter); + arma::mat draws_det_ur_sigma = arma::zeros(n_det_ur, iter); + arma::mat draws_sigma_u, draws_sigma_sigma; + if (sv || covar) { + draws_sigma_u = arma::zeros(k_dom * k_dom * tt, iter); + } else { + draws_sigma_u = arma::zeros(k_dom * k_dom, iter); + } + if (sv) { + draws_sigma_sigma = arma::zeros(k_dom * k_dom, iter); + } + + arma::vec gamma_lambda_vec, psi_lambda_vec; + arma::mat draws_a0_lambda, draws_dom_lambda, draws_for_lambda, draws_glo_lambda, draws_det_ur_lambda; + if (bvs) { + if (structural) { + draws_a0_lambda = arma::zeros(n_a0, iter); + } + draws_dom_lambda = arma::zeros(n_glo, iter); + draws_for_lambda = arma::zeros(n_for, iter); + draws_glo_lambda = arma::zeros(n_glo, iter); + draws_det_ur_lambda = arma::zeros(n_det_ur, iter); + } + if (covar && psi_bvs) { + draws_a0_lambda = arma::zeros(n_psi, iter); + } + + ///////////////////////////////////////////// Gibbs sampler /////////////////////////////////////////////////////// + for (int draw = 0; draw < draws; draw++) { + + if (draw % 20 == 0) { // Check for user interruption ever now and then + Rcpp::checkUserInterrupt(); + } + + /////////////////////////////////////////////////////////////////////// + // Draw non-cointegration parameters + /////////////////////////////////////////////////////////////////////// + + if (use_rr) { + // Update ECT + if (bvs) { + for (int i = 0; i < tt; i++) { + zz_bvs.submat(i * k_dom, 0, (i + 1) * k_dom - 1, n_alpha - 1) = arma::kron(arma::trans(arma::trans(arma::reshape(beta.col(i), k_w, r)) * w.col(i)), diag_k); + } + } else { + for (int i = 0; i < tt; i++) { + z.submat(i * k_dom, 0, (i + 1) * k_dom - 1, n_alpha - 1) = arma::kron(arma::trans(arma::trans(arma::reshape(beta.col(i), k_w, r)) * w.col(i)), diag_k); + } + } + } + + if (bvs) { + z = zz_bvs * gamma_lambda; // Update zz for BVS + } + + // Draw gamma + gamma = bvartools::kalman_dk(y, z, sigma_u, gamma_sigma_v, gamma_b, gamma_init, gamma_sigma_v); + gamma = gamma.cols(1, tt); + + // Draw sigma_v_i + gamma_lag.col(0) = gamma_init; + gamma_lag.cols(1, tt - 1) = gamma.cols(0, tt - 2); + gamma_v = arma::trans(gamma - gamma_lag); + gamma_sigma_v_post_scale = 1 / (gamma_sigma_v_prior_rate + arma::vectorise(arma::sum(arma::pow(gamma_v, 2))) * 0.5); + for (int i = 0; i < n_tot; i++) { + gamma_sigma_v_i(i, i) = arma::randg(arma::distr_param(gamma_sigma_v_post_shape(i), gamma_sigma_v_post_scale(i))); + gamma_sigma_v(i, i) = 1 / gamma_sigma_v_i(i, i); + } + + // Draw gamma_0 + if (use_rr) { + // Update alpha_0 prior + prior_gamma_vinv.submat(alpha_pos_start, alpha_pos_start, alpha_pos_end, alpha_pos_end) = 1 / (1 - rho * rho) * arma::eye(n_alpha, n_alpha); + } + post_gamma0_v = prior_gamma_vinv + gamma_sigma_v_i; + post_gamma_mu = arma::solve(post_gamma0_v, prior_gamma_vinv * prior_gamma_mu + gamma_sigma_v_i * gamma.submat(0, 0, n_tot - 1, 0)); + gamma_init = post_gamma_mu + arma::solve(arma::chol(post_gamma0_v), arma::randn(n_tot)); + + // Variable selection + if (n_nonalpha > 0 && bvs) { + z = zz_bvs; + gamma_AG = gamma_lambda * gamma; // Old selection + gamma_varsel_include_draw = shuffle(gamma_varsel_include); // Reorder positions of variable selection + for (int j = 0; j < gamma_varsel_n; j++){ // Repeat for each variable + gamma_varsel_pos = gamma_varsel_include_draw(j); + gamma_randu = arma::log(arma::randu(1)); + if (gamma_lambda(gamma_varsel_pos, gamma_varsel_pos) == 1 && gamma_randu(0) >= gamma_bvs_lprior_1(gamma_varsel_pos)){continue;} + if (gamma_lambda(gamma_varsel_pos, gamma_varsel_pos) == 0 && gamma_randu(0) >= gamma_bvs_lprior_0(gamma_varsel_pos)){continue;} + if ((gamma_lambda(gamma_varsel_pos, gamma_varsel_pos) == 1 && gamma_randu(0) < gamma_bvs_lprior_1(gamma_varsel_pos)) || (gamma_lambda(gamma_varsel_pos, gamma_varsel_pos) == 0 && gamma_randu(0) < gamma_bvs_lprior_0(gamma_varsel_pos))){ + // Candidate exclude + gamma_theta0 = gamma_AG; + gamma_theta0.row(gamma_varsel_pos) = arma::zeros(1, tt); + // Candidate include + gamma_theta1 = gamma_AG; + gamma_theta1.row(gamma_varsel_pos) = gamma.row(gamma_varsel_pos); + // Obtain errors + gamma_l0 = 0; + gamma_l1 = 0; + for (int i = 0; i < tt; i++) { + gamma_bvs_l0_res.col(i) = y.col(i) - z.rows(i * k_dom, (i + 1) * k_dom - 1) * gamma_theta0.col(i); + gamma_bvs_l1_res.col(i) = y.col(i) - z.rows(i * k_dom, (i + 1) * k_dom - 1) * gamma_theta1.col(i); + gamma_l0 = gamma_l0 + arma::as_scalar(arma::trans(gamma_bvs_l0_res.col(i)) * sigma_u_i.rows(i * k_dom, (i + 1) * k_dom - 1) * gamma_bvs_l0_res.col(i)); + gamma_l1 = gamma_l0 + arma::as_scalar(arma::trans(gamma_bvs_l1_res.col(i)) * sigma_u_i.rows(i * k_dom, (i + 1) * k_dom - 1) * gamma_bvs_l1_res.col(i)); + } + gamma_l0 = -gamma_l0 * 0.5 + arma::as_scalar(gamma_bvs_lprior_0(gamma_varsel_pos)); + gamma_l1 = -gamma_l1 * 0.5 + arma::as_scalar(gamma_bvs_lprior_1(gamma_varsel_pos)); + gamma_bayes = gamma_l1 - gamma_l0; + gamma_bayes_rand = arma::as_scalar(arma::log(arma::randu(1))); + if (gamma_bayes >= gamma_bayes_rand){ + gamma_lambda(gamma_varsel_pos, gamma_varsel_pos) = 1; + } else { + gamma_lambda(gamma_varsel_pos, gamma_varsel_pos) = 0; + } + } + } + gamma = gamma_lambda * gamma; + gamma_lambda_vec = gamma_lambda.diag(); + } + + //////// Draw cointegration parameters ///////////// + if (use_rr) { + + // Get y for beta + if (n_nonalpha > 0) { + for (int i = 0; i < tt; i++) { + beta_y.col(i) = y.col(i) - z.submat(i * k_dom, n_alpha, (i + 1) * k_dom - 1, n_tot - 1) * gamma.submat(n_alpha, i, n_tot - 1, i); + } + } else { + beta_y = y; + } + + for (int i = 0; i < tt; i++) { + beta_z.rows(i * k_dom, (i + 1) * k_dom - 1) = arma::kron(arma::reshape(gamma.submat(0, i, n_alpha - 1, i), k_dom, r), arma::trans(w.col(i))); + } + + beta = bvartools::kalman_dk(beta_y, beta_z, sigma_u, beta_sigma, beta_b, beta_init, beta_sigma); + beta = beta.cols(1, tt); + + for (int i = 0; i < tt; i++) { + pi_temp = arma::reshape(gamma.submat(0, i, n_alpha - 1, i), k_dom, r) * arma::trans(arma::reshape(beta.col(i), k_w, r)); + u.col(i) = y.col(i) - pi_temp * w.col(i); // Subtract cointegration part from y + if (n_nonalpha > 0) { // If estimated, subtract non-cointegration effects + u.col(i) = u.col(i) - z.submat(i * k_dom, n_alpha, (i + 1) * k_dom - 1, n_tot - 1) * gamma.submat(n_alpha, i, n_tot - 1, i); + } + pi.col(i) = arma::vectorise(pi_temp); // Update final Pi matrix + } + + /////// Draw beta_0 ////////// + beta_v_post = beta_vinv_prior + beta_sigma; + beta_mu_post = arma::solve(beta_v_post, beta_vinv_prior * beta_mu_prior + beta_sigma * beta.col(0)); + beta_init = beta_mu_post + arma::solve(arma::chol(beta_v_post), arma::randn(n_beta)); + + } else { + if (n_nonalpha > 0) { + for (int i = 0; i < tt; i++) { + u.col(i) = y.col(i) - z.rows(i * k_dom, (i + 1) * k_dom - 1) * gamma.col(i); + } + } else { + u = y; + } + } + + // Covariances + if (covar) { + + // Prepare data + for (int j = 0; j < tt; j++) { + for (int i = 1; i < k_dom; i++) { + psi_z.submat(j * (k_dom - 1) + i - 1, + i * (i - 1) / 2, + j * (k_dom - 1) + i - 1, + (i + 1) * i / 2 - 1) = -arma::trans(u.submat(0, j, i - 1, j)); + } + omega_psi.rows(j * (k_dom - 1), (j + 1) * (k_dom - 1) - 1) = arma::solve(omega_i.submat(j * k_dom + 1, 1, (j + 1) * k_dom - 1, k_dom - 1), arma::eye(k_dom - 1, k_dom - 1)); + } + + if (psi_bvs) { + psi_z_bvs = psi_z; + psi_z = psi_z_bvs * psi_lambda; + } + + psi_y = u.rows(1, k_dom - 1); + psi = bvartools::kalman_dk(psi_y, psi_z, omega_psi, psi_sigma_v, psi_b, psi_init, psi_sigma_v); + psi = psi.cols(1, tt); + + if (psi_bvs) { + + // Reorder positions of variable selection + psi_varsel_include_draw = shuffle(psi_varsel_include); + + psi_z = psi_z_bvs; + psi_AG = psi_lambda * psi; + for (int j = 0; j < psi_varsel_n; j++){ + psi_varsel_pos = psi_varsel_include_draw(j); + psi_randu = arma::log(arma::randu(1)); + if (psi_lambda(psi_varsel_pos, psi_varsel_pos) == 1 && psi_randu(0) >= psi_bvs_lprior_1(psi_varsel_pos)){continue;} + if (psi_lambda(psi_varsel_pos, psi_varsel_pos) == 0 && psi_randu(0) >= psi_bvs_lprior_0(psi_varsel_pos)){continue;} + if ((psi_lambda(psi_varsel_pos, psi_varsel_pos) == 1 && psi_randu(0) < psi_bvs_lprior_1(psi_varsel_pos)) || (psi_lambda(psi_varsel_pos, psi_varsel_pos) == 0 && psi_randu(0) < psi_bvs_lprior_0(psi_varsel_pos))){ + // Candidate exclude + psi_theta0 = psi_AG; + psi_theta0.row(psi_varsel_pos) = arma::zeros(1, tt);; + // Candidate include + psi_theta1 = psi_AG; + psi_theta1.row(psi_varsel_pos) = psi.row(psi_varsel_pos); + // Obtain errors + psi_l0 = 0; + psi_l1 = 0; + for (int i = 0; i < tt; i++) { + psi_theta0_res.col(i) = psi_y.col(i) - psi_z.rows(i * (k_dom - 1), (i + 1) * (k_dom - 1) - 1) * psi_theta0.col(i); + psi_theta1_res.col(i) = psi_y.col(i) - psi_z.rows(i * (k_dom - 1), (i + 1) * (k_dom - 1) - 1) * psi_theta1.col(i); + psi_l0 = psi_l0 + arma::as_scalar(arma::trans(psi_theta0_res.col(i)) * omega_psi.rows(i * (k_dom - 1), (i + 1) * (k_dom - 1) - 1) * psi_theta0_res.col(i)); + psi_l1 = psi_l1 + arma::as_scalar(arma::trans(psi_theta1_res.col(i)) * omega_psi.rows(i * (k_dom - 1), (i + 1) * (k_dom - 1) - 1) * psi_theta1_res.col(i)); + } + psi_bayes = psi_l1 - psi_l0; + psi_bayes_rand = arma::as_scalar(arma::log(arma::randu(1))); + if (psi_bayes >= psi_bayes_rand){ + psi_lambda(psi_varsel_pos, psi_varsel_pos) = 1; + } else { + psi_lambda(psi_varsel_pos, psi_varsel_pos) = 0; + } + } + } + psi = psi_lambda * psi; + psi_lambda_vec = psi_lambda.diag(); + } + + psi_lag.col(0) = psi_init; + psi_lag.cols(1, tt - 1) = psi.cols(0, tt - 2); + psi_v = arma::trans(psi - psi_lag); + psi_sigma_v_post_scale = 1 / (psi_sigma_v_prior_rate + arma::vectorise(arma::sum(arma::pow(psi_v, 2))) * 0.5); + for (int i = 0; i < n_psi; i++) { + psi_sigma_v_i(i, i) = arma::randg(arma::distr_param(psi_sigma_v_post_shape(i), psi_sigma_v_post_scale(i))); + psi_sigma_v(i, i) = 1 / psi_sigma_v_i(i, i); + } + + // Draw initial state of a + psi_init_post_v = psi_init_prior_vi + psi_sigma_v_i; + psi_init_post_mu = arma::solve(psi_init_post_v, psi_init_prior_vi * psi_init_prior_mu + psi_sigma_v_i * psi.col(0)); + psi_init = psi_init_post_mu + arma::solve(arma::chol(psi_init_post_v), arma::randn(n_psi)); + + for (int j = 0; j < tt; j ++) { + for (int i = 1; i < k_dom; i++) { + Psi.submat((k_dom * j) + i, 0, (k_dom * j) + i, i - 1) = arma::trans(psi.submat(i * (i - 1) / 2, j, (i + 1) * i / 2 - 1, j)); + } + u.col(j) = Psi.rows(j * k_dom, (j + 1) * k_dom - 1) * u.col(j); + } + } + + /////////////////////////////////////////////////////////////////////// + // Draw error variances + /////////////////////////////////////////////////////////////////////// + if (sv) { + + // Draw variances + for (int i = 0; i < k_dom; i++) { + h.col(i) = bvartools::stoch_vol(u.row(i).t(), h.col(i), sigma_h(i), h_init(i), h_const(i)); + } + + // Update sigma_u + for (int i = 0; i < tt; i++) { + omega_i.rows(i * k_dom, (i + 1) * k_dom - 1) = arma::diagmat(1 / arma::exp(h.row(i))); + if (covar) { + sigma_u_i.rows(i * k_dom, (i + 1) * k_dom - 1) = arma::trans(Psi.rows(i * k_dom, (i + 1) * k_dom - 1)) * omega_i.rows(i * k_dom, (i + 1) * k_dom - 1) * Psi.rows(i * k_dom, (i + 1) * k_dom - 1); + sigma_u.rows(i * k_dom, (i + 1) * k_dom - 1) = arma::solve(sigma_u_i.rows(i * k_dom, (i + 1) * k_dom - 1), diag_k); + } else { + sigma_u.rows(i * k_dom, (i + 1) * k_dom - 1) = arma::diagmat(arma::exp(h.row(i))); + } + } + + // Draw sigma_h + h_lag.row(0) = h_init.t(); + h_lag.rows(1, tt - 1) = h.rows(0, tt - 2); + h_lag = h - h_lag; + sigma_post_scale = 1 / (sigma_prior_rate + arma::trans(arma::sum(arma::pow(h_lag, 2))) * 0.5); + for (int i = 0; i < k_dom; 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_dom)); + + } else { + + // Obtain squared errors + sse = u * u.t(); + + if (use_gamma) { + // Draw from gamma distribution + for (int i = 0; i < k_dom; 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))); + } + + if (covar) { + for (int i = 0; i < tt; i++) { + sigma_u_i.rows(i * k_dom, (i + 1) * k_dom - 1) = arma::trans(Psi.rows(i * k_dom, (i + 1) * k_dom - 1)) * omega_i.rows(i * k_dom, (i + 1) * k_dom - 1) * Psi.rows(i * k_dom, (i + 1) * k_dom - 1); + sigma_u.rows(i * k_dom, (i + 1) * k_dom - 1) = arma::solve(sigma_u_i.rows(i * k_dom, (i + 1) * k_dom - 1), diag_k); + } + } else { + sigma_u_i = omega_i; + sigma_u_temp = arma::solve(sigma_u_i, diag_k); + for (int i = 0; i < tt; i++) { + sigma_u.rows(i * k_dom, (i + 1) * k_dom - 1) = sigma_u_temp; + } + } + + } else { + sigma_u_i = arma::wishrnd(arma::solve(sigma_prior_scale + sse, diag_k), sigma_post_df); + sigma_u_temp = arma::solve(sigma_u_i, diag_k); + for (int i = 0; i < tt; i++) { + sigma_u.rows(i * k_dom, (i + 1) * k_dom - 1) = sigma_u_temp; + } + } + } + + + + // Store draws + if (draw >= burnin) { + + pos_draw = draw - burnin; + + if (sv) { + for (int i = 0; i < tt; i ++) { + draws_sigma_u.submat(i * n_sigma, pos_draw, (i + 1) * n_sigma - 1, pos_draw) = arma::vectorise(sigma_u.rows(i * k_dom, (i + 1) * k_dom - 1)); + } + draws_sigma_sigma.col(pos_draw) = arma::vectorise(arma::diagmat(sigma_h)); + } else { + draws_sigma_u.col(pos_draw) = arma::vectorise(sigma_u.rows(0, k_dom - 1)); + } + + if (psi_bvs) { + draws_a0_lambda.col(pos_draw) = psi_lambda_vec; + } + + if (use_rr) { + draws_alpha.col(pos_draw) = arma::vectorise(gamma.rows(alpha_pos_start, alpha_pos_end)); + for (int i = 0; i < tt; i++) { + beta_t = arma::reshape(beta.col(i), k_w, r).t(); + draws_beta_dom.submat(i * k_dom * r, pos_draw, (i + 1) * r * k_dom - 1, pos_draw) = arma::vectorise(arma::trans(beta_t.cols(0, k_dom - 1))); + if (k_for > 0) { + draws_beta_for.submat(i * k_for * r, pos_draw, (i + 1) * r * k_for - 1, pos_draw) = arma::vectorise(arma::trans(beta_t.cols(k_dom, k_dom + k_for - 1))); + } + if (k_glo > 0) { + draws_beta_glo.submat(i * k_glo * r, pos_draw, (i + 1) * r * k_glo - 1, pos_draw) = arma::vectorise(arma::trans(beta_t.cols(k_dom + k_for, k_dom + k_for + k_glo - 1))); + } + if (k_det_r > 0) { + draws_beta_det.submat(i * k_det_r * r, pos_draw, (i + 1) * r * k_det_r - 1, pos_draw) = arma::vectorise(arma::trans(beta_t.cols(k_dom + k_for + k_glo, k_dom + k_for + k_glo + k_det_r - 1))); + } + } + } + + if (n_dom > 0) { + draws_dom.col(pos_draw) = arma::vectorise(gamma.rows(dom_pos_start, dom_pos_end)); + draws_dom_sigma.col(pos_draw) = 1 / arma::mat(gamma_sigma_v_i.submat(dom_pos_start, dom_pos_start, dom_pos_end, dom_pos_end)).diag(); + if (bvs) { + draws_dom_lambda.col(pos_draw) = gamma_lambda_vec.subvec(dom_pos_start, dom_pos_end); + } + } + + if (n_for > 0) { + draws_for.col(pos_draw) = arma::vectorise(gamma.rows(for_pos_start, for_pos_end)); + draws_for_sigma.col(pos_draw) = 1 / arma::mat(gamma_sigma_v_i.submat(for_pos_start, for_pos_start, for_pos_end, for_pos_end)).diag(); + if (bvs) { + draws_for_lambda.col(pos_draw) = gamma_lambda_vec.subvec(for_pos_start, for_pos_end); + } + } + + if (n_glo > 0) { + draws_glo.col(pos_draw) = arma::vectorise(gamma.rows(glo_pos_start, glo_pos_end)); + draws_glo_sigma.col(pos_draw) = 1 / arma::mat(gamma_sigma_v_i.submat(glo_pos_start, glo_pos_start, glo_pos_end, glo_pos_end)).diag(); + if (bvs) { + draws_glo_lambda.col(pos_draw) = gamma_lambda_vec.subvec(glo_pos_start, glo_pos_end); + } + } + + if (n_det_ur > 0) { + draws_det_ur.col(pos_draw) = arma::vectorise(gamma.rows(c_pos_start, c_pos_end)); + draws_det_ur_sigma.col(pos_draw) = 1 / arma::mat(gamma_sigma_v_i.submat(c_pos_start, c_pos_start, c_pos_end, c_pos_end)).diag(); + if (bvs) { + draws_det_ur_lambda.col(pos_draw) = gamma_lambda_vec.subvec(c_pos_start, c_pos_end); + } + } + + if (structural) { + draws_a0.col(pos_draw) = arma::vectorise(gamma.rows(a0_pos_start, a0_pos_end)); + draws_a0_sigma.col(pos_draw) = 1 / arma::mat(gamma_sigma_v_i.submat(a0_pos_start, a0_pos_start, a0_pos_end, a0_pos_end)).diag(); + if (bvs) { + draws_a0_lambda.col(pos_draw) = gamma_lambda_vec.subvec(a0_pos_start, a0_pos_end); + } + } + } // End posterior storing + } // End Gibbs sampler + + Rcpp::List posteriors = Rcpp::List::create(Rcpp::Named("a0") = R_NilValue, + Rcpp::Named("alpha") = R_NilValue, + Rcpp::Named("beta_dom") = R_NilValue, + Rcpp::Named("beta_for") = R_NilValue, + Rcpp::Named("beta_glo") = R_NilValue, + Rcpp::Named("beta_det") = R_NilValue, + Rcpp::Named("gamma_dom") = R_NilValue, + Rcpp::Named("gamma_for") = R_NilValue, + Rcpp::Named("gamma_glo") = R_NilValue, + Rcpp::Named("gamma_det") = R_NilValue, + Rcpp::Named("sigma") = R_NilValue); + + if (use_rr) { + posteriors["alpha"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_alpha)); + + // Reformat draws + posteriors["beta_dom"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_beta_dom)); + if (k_for > 0) { + posteriors["beta_for"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_beta_for)); + } + if (k_glo > 0) { + posteriors["beta_glo"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_beta_glo)); + } + if (k_det_r > 0) { + posteriors["beta_det"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_beta_det)); + } + } + + if (n_dom > 0) { + if (bvs) { + posteriors["gamma_dom"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_dom, + Rcpp::Named("sigma") = draws_dom_sigma, + Rcpp::Named("lambda") = draws_dom_lambda)); + } else { + posteriors["gamma_dom"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_dom, + Rcpp::Named("sigma") = draws_dom_sigma)); + } + } + + if (n_for > 0) { + if (bvs) { + posteriors["gamma_for"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_for, + Rcpp::Named("sigma") = draws_for_sigma, + Rcpp::Named("lambda") = draws_for_lambda)); + } else { + posteriors["gamma_for"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_for, + Rcpp::Named("sigma") = draws_for_sigma)); + } + } + + if (n_glo > 0) { + if (bvs) { + posteriors["gamma_glo"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_glo, + Rcpp::Named("sigma") = draws_glo_sigma, + Rcpp::Named("lambda") = draws_glo_lambda)); + } else { + posteriors["gamma_glo"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_glo, + Rcpp::Named("sigma") = draws_glo_sigma)); + } + } + + if (n_det_ur > 0) { + if (bvs) { + posteriors["gamma_det"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_det_ur, + Rcpp::Named("sigma") = draws_det_ur_sigma, + Rcpp::Named("lambda") = draws_det_ur_lambda)); + } else { + posteriors["gamma_det"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_det_ur, + Rcpp::Named("sigma") = draws_det_ur_sigma)); + } + } + + if (structural) { + if (bvs) { + posteriors["a0"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_a0, + Rcpp::Named("sigma") = draws_a0_sigma, + Rcpp::Named("lambda") = draws_a0_lambda)); + } else { + posteriors["a0"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_a0, + Rcpp::Named("sigma") = draws_a0_sigma)); + } + } + + if (psi_bvs) { + if (sv) { + posteriors["sigma"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_sigma_u, + Rcpp::Named("sigma") = draws_sigma_sigma, + Rcpp::Named("lambda") = draws_a0_lambda)); + } else { + posteriors["sigma"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_sigma_u, + Rcpp::Named("lambda") = draws_a0_lambda)); + } + } else { + if (sv) { + posteriors["sigma"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_sigma_u, + Rcpp::Named("sigma") = draws_sigma_sigma)); + } else { + posteriors["sigma"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_sigma_u)); + } + } + + return Rcpp::List::create(Rcpp::Named("data") = object["data"], + Rcpp::Named("model") = object["model"], + Rcpp::Named("priors") = object["priors"], + Rcpp::Named("posteriors") = posteriors); + +} + + +/*** R + +#library(bgvars) +library(Matrix) + +set.seed(123456) + +data("gvar2019") + +# Create regions +temp <- create_regions(country_data = gvar2019$country_data, + weight_data = gvar2019$weight_data, + region_weights = gvar2019$region_weights, + regions = list(EA = c("AT", "BE", "DE", "ES", "FI", "FR", "IT", "NL")), period = 3) + +country_data <- temp$country_data +weight_data <- temp$weight_data + +global_data = gvar2019$global_data + +# Create weights +model_weights <- create_weights(weight_data, period = 3, country_data = country_data) + +# Create models +model_specs <- create_specifications(country_data = country_data, + global_data = global_data, + domestic = list(variables = c("y", "Dp", "r", "lr", "eq", "ep"), lags = 1), + foreign = list(variables = c("y", "Dp", "eq", "r", "lr"), lags = 1), + global = list(variables = "poil", lags = 1), + deterministic = list(const = "unrestricted", + trend = "restricted"), + countries = NULL, + tvp = TRUE, sv = TRUE, + type = "VEC", r = 1, + iterations = 100, + burnin = 100) + +country_models <- create_models(country_data = country_data, + weight_data = model_weights, + global_data = global_data, + model_specs = model_specs) + +# Add priors +models_with_priors <- add_priors(country_models, + coef = list(v_i = 1 / 9, v_i_det = 1 / 9, shape = 3, rate = .0001), + sigma = list(shape = 3, rate = .0001, mu = -9, v_i = .1, sigma_h = .05, constant = .0001, covar = TRUE), + bvs = list(inprior = .5, covar = TRUE)) + +temp <- .bgvectvpalg(models_with_priors[["EA"]]) +a <- temp + +*/