From daebf776309542a540e990a3537b442aca1eb77d Mon Sep 17 00:00:00 2001 From: Franz Mohr Date: Sat, 26 Aug 2023 13:01:47 +0200 Subject: [PATCH] Drastic measure --- DESCRIPTION | 12 - R/RcppExports.R | 12 - R/add_priors.gvecsubmodels.R | 867 --------------------------- R/bvecxpost.R | 359 ----------- R/create_pi_matrices.R | 159 ----- R/ctryvec_to_ctryvar.R | 376 ------------ R/draw_posterior.gvecsubmodels.R | 89 --- R/plot.bgvecest.R | 30 - R/plot.ctryvecest.R | 485 --------------- R/print.summary.bgvecest.R | 15 - R/print.summary.ctryvecest.R | 113 ---- R/summary.bgvecest.R | 45 -- R/summary.ctryvecest.R | 215 ------- R/thin.bgvecest.R | 74 --- src/RcppExports.cpp | 36 -- src/bgvartvpalg.cpp | 737 ----------------------- src/bgvecalg.cpp | 890 ---------------------------- src/bgvectvpalg.cpp | 985 ------------------------------- 18 files changed, 5499 deletions(-) delete mode 100644 R/add_priors.gvecsubmodels.R delete mode 100644 R/bvecxpost.R delete mode 100644 R/create_pi_matrices.R delete mode 100644 R/ctryvec_to_ctryvar.R delete mode 100644 R/draw_posterior.gvecsubmodels.R delete mode 100644 R/plot.bgvecest.R delete mode 100644 R/plot.ctryvecest.R delete mode 100644 R/print.summary.bgvecest.R delete mode 100644 R/print.summary.ctryvecest.R delete mode 100644 R/summary.bgvecest.R delete mode 100644 R/summary.ctryvecest.R delete mode 100644 R/thin.bgvecest.R delete mode 100644 src/bgvartvpalg.cpp delete mode 100644 src/bgvecalg.cpp delete mode 100644 src/bgvectvpalg.cpp diff --git a/DESCRIPTION b/DESCRIPTION index 28dc153..66a62fe 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -24,22 +24,17 @@ 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' - 'create_pi_matrices.R' 'create_regions.R' 'create_specifications.R' 'create_weights.R' - 'ctryvec_to_ctryvar.R' 'data.R' 'diff_variables.R' 'draw_posterior.gvarsubmodels.R' - 'draw_posterior.gvecsubmodels.R' 'gen_varx.R' 'gen_vecx.R' 'get_regressor_names.R' @@ -53,23 +48,16 @@ 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/R/RcppExports.R b/R/RcppExports.R index e4e5716..b0947d8 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -5,18 +5,6 @@ .Call(`_bgvars_bgvaralg`, object) } -.bgvartvpalg <- function(object) { - .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 deleted file mode 100644 index 62c3d3b..0000000 --- a/R/add_priors.gvecsubmodels.R +++ /dev/null @@ -1,867 +0,0 @@ -#' 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 deleted file mode 100644 index 9342fe9..0000000 --- a/R/bvecxpost.R +++ /dev/null @@ -1,359 +0,0 @@ -#' 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/create_pi_matrices.R b/R/create_pi_matrices.R deleted file mode 100644 index 4e84f9d..0000000 --- a/R/create_pi_matrices.R +++ /dev/null @@ -1,159 +0,0 @@ - -# Helper function to obtain Pi matrices from draws of alpha and beta_* - -.create_pi_matrices <- function(object, drop_beta = TRUE) { - - draws <- NULL - specs <- NULL - vars <- c("a0", "alpha", - "beta_domestic", "beta_foreign", "beta_global", "beta_deterministic", - "gamma_domestic", "gamma_foreign", "gamma_global", "gamma_deterministic") - for (i in vars) { - if (is.null(draws)) { - if (!is.null(object[["posteriors"]][[i]])) { - if (is.list(object[["posteriors"]][[i]])) { - draws <- nrow(object[["posteriors"]][[i]][[1]]) - } else { - draws <- nrow(object[["posteriors"]][[i]]) - } - } - } - if (is.null(specs)) { - if (is.list(object[["posteriors"]][[i]])) { - specs <- attr(object[["posteriors"]][[i]][[1]], "mcpar") - } else { - specs <- attr(object[["posteriors"]][[i]], "mcpar") - } - } - } - - k_dom <- NCOL(object[["data"]][["Y"]]) - k_for <- length(object[["model"]][["foreign"]][["variables"]]) - tt <- NROW(object[["data"]][["Y"]]) - tvp <- object[["model"]][["tvp"]] - p <- object[["model"]][["domestic"]][["lags"]] - r <- object[["model"]][["rank"]] - - ## Domestic ---- - if (is.null(object[["posteriors"]][["pi_domestic"]])) { - if (!is.null(object[["posteriors"]][["alpha"]]) & !is.null(object[["posteriors"]][["beta_domestic"]])){ - if (tvp) { - object[["posteriors"]][["pi_domestic"]] <- list() - for (j in 1:tt) { - object[["posteriors"]][["pi_domestic"]][[j]] <- matrix(NA_real_, draws, k_dom * k_dom) - for (i in 1:draws) { - object[["posteriors"]][["pi_domestic"]][[j]][i,] <- matrix(object[["posteriors"]][["alpha"]][[j]][i,], k_dom) %*% t(matrix(object[["posteriors"]][["beta_domestic"]][[j]][i,], k_dom)) - } - object[["posteriors"]][["pi_domestic"]][[j]] <- coda::mcmc(object[["posteriors"]][["pi_domestic"]][[j]]) - attr(object[["posteriors"]][["pi_domestic"]][[j]], "mcpar") <- specs - } - } else { - object[["posteriors"]][["pi_domestic"]] <- matrix(NA_real_, draws, k_dom * k_dom) - for (i in 1:draws) { - object[["posteriors"]][["pi_domestic"]][i,] <- matrix(object[["posteriors"]][["alpha"]][i,], k_dom) %*% t(matrix(object[["posteriors"]][["beta_domestic"]][i,], k_dom)) - } - object[["posteriors"]][["pi_domestic"]] <- coda::mcmc(object[["posteriors"]][["pi_domestic"]]) - attr(object[["posteriors"]][["pi_domestic"]], "mcpar") <- specs - } - } - if (drop_beta) { - object[["posteriors"]][["beta_domestic"]] <- NULL - } - } - - ## Foreign ---- - - if (is.null(object[["posteriors"]][["pi_foreign"]])) { - if (!is.null(object[["posteriors"]][["alpha"]]) & !is.null(object[["posteriors"]][["beta_foreign"]])){ - if (tvp) { - object[["posteriors"]][["pi_foreign"]] <- list() - for (j in 1:tt) { - object[["posteriors"]][["pi_foreign"]][[j]] <- matrix(NA_real_, draws, k_dom * k_for) - for (i in 1:draws) { - object[["posteriors"]][["pi_foreign"]][[j]][i,] <- matrix(object[["posteriors"]][["alpha"]][[j]][i,], k_dom) %*% t(matrix(object[["posteriors"]][["beta_foreign"]][[j]][i,], k_for)) - } - object[["posteriors"]][["pi_foreign"]][[j]] <- coda::mcmc(object[["posteriors"]][["pi_foreign"]][[j]]) - attr(object[["posteriors"]][["pi_foreign"]][[j]], "mcpar") <- specs - } - } else { - object[["posteriors"]][["pi_foreign"]] <- matrix(NA_real_, draws, k_dom * k_for) - for (i in 1:draws) { - object[["posteriors"]][["pi_foreign"]][i,] <- matrix(object[["posteriors"]][["alpha"]][i,], k_dom) %*% t(matrix(object[["posteriors"]][["beta_foreign"]][i,], k_for)) - } - object[["posteriors"]][["pi_foreign"]] <- coda::mcmc(object[["posteriors"]][["pi_foreign"]]) - attr(object[["posteriors"]][["pi_foreign"]], "mcpar") <- specs - } - } - if (drop_beta) { - object[["posteriors"]][["beta_foreign"]] <- NULL - } - } - - ## Global ---- - global <- !is.null(object[["model"]][["global"]]) - if (global) { - k_glo <- length(object[["model"]][["global"]][["variables"]]) - if (is.null(object[["posteriors"]][["pi_global"]])) { - if (!is.null(object[["posteriors"]][["alpha"]]) & !is.null(object[["posteriors"]][["beta_global"]])){ - if (tvp) { - object[["posteriors"]][["pi_global"]] <- list() - for (j in 1:tt) { - object[["posteriors"]][["pi_global"]][[j]] <- matrix(NA_real_, draws, k_dom * k_glo) - for (i in 1:draws) { - object[["posteriors"]][["pi_global"]][[j]][i,] <- matrix(object[["posteriors"]][["alpha"]][[j]][i,], k_dom) %*% t(matrix(object[["posteriors"]][["beta_global"]][[j]][i,], k_glo)) - } - object[["posteriors"]][["pi_global"]][[j]] <- coda::mcmc(object[["posteriors"]][["pi_global"]][[j]]) - attr(object[["posteriors"]][["pi_global"]][[j]], "mcpar") <- specs - } - } else { - object[["posteriors"]][["pi_global"]] <- matrix(NA_real_, draws, k_dom * k_glo) - for (i in 1:draws) { - object[["posteriors"]][["pi_global"]][i,] <- matrix(object[["posteriors"]][["alpha"]][i,], k_dom) %*% t(matrix(object[["posteriors"]][["beta_global"]][i,], k_glo)) - } - object[["posteriors"]][["pi_global"]] <- coda::mcmc(object[["posteriors"]][["pi_global"]]) - attr(object[["posteriors"]][["pi_global"]], "mcpar") <- specs - } - } - if (drop_beta) { - object[["posteriors"]][["beta_global"]] <- NULL - } - } - } - - ## Deterministic ---- - k_det_r <- 0 - if (!is.null(object[["model"]][["deterministic"]][["restricted"]])) { - k_det_r <- length(object[["model"]][["deterministic"]][["restricted"]]) - if (is.null(object[["posteriors"]][["pi_deterministic"]])) { - if (!is.null(object[["posteriors"]][["alpha"]]) & !is.null(object[["posteriors"]][["beta_deterministic"]])){ - if (tvp) { - object[["posteriors"]][["pi_deterministic"]] <- list() - for (j in 1:tt) { - object[["posteriors"]][["pi_deterministic"]][[j]] <- matrix(NA_real_, draws, k_dom * k_det_r) - for (i in 1:draws) { - object[["posteriors"]][["pi_deterministic"]][[j]][i,] <- matrix(object[["posteriors"]][["alpha"]][[j]][i,], k_dom) %*% t(matrix(object[["posteriors"]][["beta_deterministic"]][[j]][i,], k_det_r)) - } - object[["posteriors"]][["pi_deterministic"]][[j]] <- coda::mcmc(object[["posteriors"]][["pi_deterministic"]][[j]]) - attr(object[["posteriors"]][["pi_deterministic"]][[j]], "mcpar") <- specs - } - } else { - object[["posteriors"]][["pi_deterministic"]] <- matrix(NA_real_, draws, k_dom * k_det_r) - for (i in 1:draws) { - object[["posteriors"]][["pi_deterministic"]][i,] <- matrix(object[["posteriors"]][["alpha"]][i,], k_dom) %*% t(matrix(object[["posteriors"]][["beta_deterministic"]][i,], k_det_r)) - } - object[["posteriors"]][["pi_deterministic"]] <- coda::mcmc(object[["posteriors"]][["pi_deterministic"]]) - attr(object[["posteriors"]][["pi_deterministic"]], "mcpar") <- specs - } - } - if (drop_beta) { - object[["posteriors"]][["beta_deterministic"]] <- NULL - } - } - } - - if (drop_beta) { - object[["posteriors"]][["alpha"]] <- NULL - } - - return(object) -} \ No newline at end of file diff --git a/R/ctryvec_to_ctryvar.R b/R/ctryvec_to_ctryvar.R deleted file mode 100644 index 0baa18d..0000000 --- a/R/ctryvec_to_ctryvar.R +++ /dev/null @@ -1,376 +0,0 @@ -#' Transform the BVEC Submodels of a BGVAR model to a BVAR Submodel in Levels -#' -#' An object of class \code{"bgvecest"} is transformed to a VAR in level representation. -#' -#' @param object an object of class \code{"bgvecest"}. -#' @param ctry character. Name of the element in argument \code{object}, for which -#' the transformation should be carried out. If \code{NULL} (default), all submodels -#' are transformed. -#' -#' @return An object of class \code{"bgvarest"}. -#' -#' @references -#' -#' Lütkepohl, H. (2006). \emph{New introduction to multiple time series analysis} (2nd ed.). Berlin: Springer. -#' -#' @export -ctryvec_to_ctryvar <- function(object, ctry = NULL) { - - if (!any(c("bgvecest", "ctryvecest") %in% class(object))) { - stop("Argument 'object' must be of class 'bgvecest' or 'ctryvecest'.") - } - - if ("bgvecest" %in% class(object)) { - if (!is.null(ctry)) { - # Determine position of countries in object - pos <- which(names(object) %in% ctry) - if (length(pos) == 0) { - stop("Specified countries not available.") - } - # Perform transformations for selected submodels - result <- NULL - for (i in pos) { - temp <- .transform_ctryvec(object[[i]]) - result <- c(result, list(temp)) - rm(temp) - } - names(result) <- names(object)[pos] - object <- result - rm(result) - } else { - object <- lapply(object, .transform_ctryvec) - } - class(object) <- c("bgvarest", "list") - } - - if ("ctryvecest" %in% class(object)) { - if (!is.null(ctry)) { - warning("Only one model provided. Argument 'ctry' will be ignored.") - } - object <- .transform_ctryvec(object) - class(object) <- c("ctryvarest", "list") - } - - return(object) -} - -.transform_ctryvec <- function(object) { - - # Skip tests if posterior simulation was not successful - cond <- is.null(object[["error"]]) - if (!cond) { - if (object[["error"]]) { - cond <- FALSE - } - } - - if (cond) { - # Get number of draws and draw information - draws <- NULL - specs <- NULL - vars <- c("a0", "alpha", - "beta_domestic", "beta_foreign", "beta_global", "beta_deterministic", - "gamma_domestic", "gamma_foreign", "gamma_global", "gamma_deterministic") - for (i in vars) { - if (is.null(draws)) { - if (!is.null(object[["posteriors"]][[i]])) { - if (is.list(object[["posteriors"]][[i]])) { - draws <- nrow(object[["posteriors"]][[i]][[1]]) - } else { - draws <- nrow(object[["posteriors"]][[i]]) - } - } - } - if (is.null(specs)) { - if (is.list(object[["posteriors"]][[i]])) { - specs <- attr(object[["posteriors"]][[i]][[1]], "mcpar") - } else { - specs <- attr(object[["posteriors"]][[i]], "mcpar") - } - } - } - - # Model specs - k <- NCOL(object[["data"]][["Y"]]) - tt <- NROW(object[["data"]][["Y"]]) - tvp <- object[["model"]][["tvp"]] - p <- object[["model"]][["domestic"]][["lags"]] - r <- object[["model"]][["rank"]] - - # Calculate Pi matrices ---- - # Obtain Pi matrices - object <- .create_pi_matrices(object) - - # Produce VAR matrices ---- - - ## Domestic ---- - A <- NULL - if (!is.null(object[["posteriors"]][["gamma_domestic"]])) { - - W <- diag(-1, k * p) - W[1:k, 1:k] <- diag(1, k) - W[-(1:k), -(k * (p - 1) + 1:k)] <- W[-(1:k),-(k * (p - 1) + 1:k)] + diag(k * (p - 1)) - J <- matrix(0, k, k * p) - J[1:k, 1:k] <- diag(1, k) - - n_gamma <- k * k * p - - if (tvp) { - - A <- list() - for (i in 1:tt) { - pi_temp <- matrix(0, k, k) - A[[i]] <- matrix(NA, n_gamma, draws) - for (draw in 1:draws) { - if (!is.null(object[["posteriors"]][["pi_domestic"]])) { - if (tvp) { - pi_temp <- matrix(object[["posteriors"]][["pi_domestic"]][[i]][draw, ], k) - } else { - pi_temp <- matrix(object[["posteriors"]][["pi_domestic"]][draw, ], k) - } - } - if (tvp) { - gamma_temp <- matrix(object[["posteriors"]][["gamma_domestic"]][[i]][draw, ], k) - } else { - gamma_temp <- matrix(object[["posteriors"]][["gamma_domestic"]][draw, ], k) - } - A[[i]][, draw] <- cbind(pi_temp, gamma_temp) %*% W + J - } - A[[i]] <- t(A[[i]]) - } - - } else { - - A <- matrix(NA, n_gamma, draws) - for (draw in 1:draws) { - if (is.null(object[["posteriors"]][["pi_domestic"]])) { - A[, draw] <- cbind(matrix(0, k, k), matrix(object[["posteriors"]][["gamma_domestic"]][draw, ], k)) %*% W + J - } else { - A[, draw] <- cbind(matrix(object[["posteriors"]][["pi_domestic"]][draw, ], k), matrix(object[["posteriors"]][["gamma_domestic"]][draw, ], k)) %*% W + J - } - } - A <- t(A) - } - - } else { - - if (!is.null(object[["posteriors"]][["pi_domestic"]])) { - n_a <- k * k - if (tvp) { - A <- list() - for (i in 1:tt) { - A[[i]] <- matrix(NA, n_a, draws) - for (draw in 1:draws) { - A[[i]][, draw] <- matrix(object[["posteriors"]][["pi_domestic"]][[i]][draw, ], k) + matrix(diag(1, k), k) - } - A[[i]] <- t(A[[i]]) - } - } else { - A <- matrix(NA, n_a, draws) - for (draw in 1:draws) { - A[, draw] <- matrix(object[["posteriors"]][["pi_domestic"]][draw, ], k) + matrix(diag(1, k), k) - } - A <- t(A) - } - } else { - A <- matrix(0, k * k, draws) - } - } - - object[["posteriors"]][["domestic"]] <- A - rm(A) - object[["posteriors"]][["pi_domestic"]] <- NULL - object[["posteriors"]][["gamma_domestic"]] <- NULL - object[["posteriors"]][["gamma_domestic_lambda"]] <- NULL - object[["posteriors"]][["gamma_domestic_sigma"]] <- NULL - - ## Foreign ---- - B <- NULL - if (!is.null(object[["posteriors"]][["gamma_foreign"]])) { - - k_for <- length(object[["model"]][["foreign"]][["variables"]]) - p_for <- object[["model"]][["foreign"]][["lags"]] - - W <- diag(-1, k_for * (p_for + 1)) - W[1:k_for, 1:k_for] <- 0 - W[1:k_for, k_for + 1:k_for] <- diag(1, k_for) - W[-(1:k_for), 1:(k_for * p_for)] <- W[-(1:k_for), 1:(k_for * p_for)] + diag(1, k_for * p_for) - - n_b <- k * k_for * (p_for + 1) - - if (tvp) { - B <- list() - for (i in 1:tt) { - B[[i]] <- matrix(NA, n_b, draws) - for (draw in 1:draws){ - pix_temp <- matrix(0, k, k_for) - if (!is.null(object[["posteriors"]][["pi_foreign"]])) { - if (tvp) { - pix_temp <- matrix(object[["posteriors"]][["pi_foreign"]][[i]][draw, ], k) - } else { - pix_temp <- matrix(object[["posteriors"]][["pi_foreign"]][draw, ], k) - } - } - if (tvp) { - ups_temp <- matrix(object[["posteriors"]][["gamma_foreign"]][[i]][draw, ], k) - } else { - ups_temp <- matrix(object[["posteriors"]][["gamma_foreign"]][draw, ], k) - } - B[[i]][, draw] <- cbind(pix_temp, ups_temp) %*% W - } - B[[i]] <- t(B[[i]]) - } - - } else { - - B <- matrix(NA, n_b, draws) - for (draw in 1:draws){ - if (!is.null(object[["posteriors"]][["pi_foreign"]])) { - B[, draw] <- cbind(matrix(object[["posteriors"]][["pi_foreign"]][draw, ], k), matrix(object[["posteriors"]][["gamma_foreign"]][draw, ], k)) %*% W - } else { - B[, draw] <- cbind(matrix(0, k, k_for), matrix(object[["posteriors"]][["gamma_foreign"]][draw, ], k)) %*% W - } - } - B <- t(B) - } - } - - object[["posteriors"]][["foreign"]] <- B - rm(B) - object[["posteriors"]][["pi_foreign"]] <- NULL - object[["posteriors"]][["gamma_foreign"]] <- NULL - object[["posteriors"]][["gamma_foreign_lambda"]] <- NULL - object[["posteriors"]][["gamma_foreign_sigma"]] <- NULL - - ## Global ---- - B <- NULL - if (!is.null(object[["posteriors"]][["gamma_global"]])) { - - k_glo <- length(object[["model"]][["global"]][["variables"]]) - p_glo <- object[["model"]][["global"]][["lags"]] - - W <- diag(-1, k_glo * (p_glo + 1)) - W[1:k_glo, 1:k_glo] <- 0 - W[1:k_glo, k_glo + 1:k_glo] <- diag(1, k_glo) - W[-(1:k_glo), 1:(k_glo * p_glo)] <- W[-(1:k_glo), 1:(k_glo * p_glo)] + diag(1, k_glo * p_glo) - - n_b <- k * k_glo * (p_glo + 1) - - if (tvp) { - - B <- list() - for (i in 1:tt) { - B[[i]] <- matrix(NA, n_b, draws) - for (draw in 1:draws){ - pix_temp <- matrix(0, k, k_glo) - if (!is.null(object[["posteriors"]][["pi_global"]])) { - if (tvp) { - pix_temp <- matrix(object[["posteriors"]][["pi_global"]][[i]][draw, ], k) - } else { - pix_temp <- matrix(object[["posteriors"]][["pi_global"]][draw, ], k) - } - } - if (tvp) { - ups_temp <- matrix(object[["posteriors"]][["gamma_global"]][[i]][draw, ], k) - } else { - ups_temp <- matrix(object[["posteriors"]][["gamma_global"]][draw, ], k) - } - B[[i]][, draw] <- cbind(pix_temp, ups_temp) %*% W - } - B[[i]] <- t(B[[i]]) - } - - } else { - - B <- matrix(NA, n_b, draws) - for (draw in 1:draws){ - if (!is.null(object[["posteriors"]][["pi_global"]])) { - B[, draw] <- cbind(matrix(object[["posteriors"]][["pi_global"]][draw, ], k), matrix(object[["posteriors"]][["gamma_global"]][draw, ], k)) %*% W - } else { - B[, draw] <- cbind(matrix(0, k, k_glo), matrix(object[["posteriors"]][["gamma_global"]][draw, ], k)) %*% W - } - } - B <- t(B) - } - object[["posteriors"]][["global"]] <- B - } - rm(B) - object[["posteriors"]][["pi_global"]] <- NULL - object[["posteriors"]][["gamma_global"]] <- NULL - object[["posteriors"]][["gamma_global_lambda"]] <- NULL - object[["posteriors"]][["gamma_global_sigma"]] <- NULL - - ## Deterministic ---- - k_det_r <- 0 - if (r > 0 & !is.null(object[["model"]][["deterministic"]][["restricted"]])) { - k_det_r <- length(object[["model"]][["deterministic"]][["restricted"]]) - } - k_det_ur <- 0 - if (!is.null(object[["model"]][["deterministic"]][["unrestricted"]])) { - k_det_ur <- length(object[["model"]][["deterministic"]][["unrestricted"]]) - } - - if (k_det_r + k_det_ur > 0) { - if (tvp) { - object[["posteriors"]][["deterministic"]] <- list() - for (i in 1:tt) { - object[["posteriors"]][["deterministic"]][[i]] <- matrix(NA, draws, (k_det_r + k_det_ur) * k) - if (k_det_r > 0) { - object[["posteriors"]][["deterministic"]][[i]][, 1:(k_det_r * k)] <- object[["posteriors"]][["pi_deterministic"]][[i]] - } - if (k_det_ur > 0) { - object[["posteriors"]][["deterministic"]][[i]][, (k_det_r * k) + 1:(k_det_ur * k)] <- object[["posteriors"]][["gamma_deterministic"]][[i]] - } - } - object[["posteriors"]][["beta_deterministic"]] <- NULL - object[["posteriors"]][["gamma_deterministic"]] <- NULL - - } else { - object[["posteriors"]][["deterministic"]] <- matrix(NA_real_, draws, (k_det_r + k_det_ur) * k) - if (k_det_r > 0) { - object[["posteriors"]][["deterministic"]][, 1:(k_det_r * k)] <- object[["posteriors"]][["pi_deterministic"]] - object[["posteriors"]][["pi_deterministic"]] <- NULL - } - if (k_det_ur > 0) { - object[["posteriors"]][["deterministic"]][, (k_det_r * k) + 1:(k_det_ur * k)] <- object[["posteriors"]][["gamma_deterministic"]] - object[["posteriors"]][["gamma_deterministic"]] <- NULL - } - } - } - - for (i in c("a0", "domestic", "foreign", "global", "deterministic")) { - if (!is.null(object[["posteriors"]][[i]])) { - if (tvp) { - for (j in 1:tt) { - object[["posteriors"]][[i]][[j]] <- coda::mcmc(object[["posteriors"]][[i]][[j]]) - attr(object[["posteriors"]][[i]][[j]], "mcpar") <- specs - } - } else { - object[["posteriors"]][[i]] <- coda::mcmc(object[["posteriors"]][[i]]) - attr(object[["posteriors"]][[i]], "mcpar") <- specs - } - } - } - - # Update raw data ---- - - # Prepare deterministic terms for .gen_varx - if (k_det_r + k_det_ur > 0) { - object[["data"]][["deterministic"]] <- cbind(object[["data"]][["deterministic"]][["restricted"]], - object[["data"]][["deterministic"]][["unrestricted"]]) - } - temp <- .gen_varx(object) - - object[["data"]][["Y"]] <- temp[["Y"]] - object[["data"]][["Z"]] <- temp[["Z"]] - object[["data"]][["SUR"]] <- temp[["SUR"]] - object[["data"]][["W"]] <- NULL - object[["data"]][["X"]] <- NULL - - object[["model"]][["type"]] <- "VAR" - - class(object) <- c("ctryvarest", "list") - } - - return(object) -} diff --git a/R/draw_posterior.gvecsubmodels.R b/R/draw_posterior.gvecsubmodels.R deleted file mode 100644 index 5b6509a..0000000 --- a/R/draw_posterior.gvecsubmodels.R +++ /dev/null @@ -1,89 +0,0 @@ -#' 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 deleted file mode 100644 index 2472264..0000000 --- a/R/plot.bgvecest.R +++ /dev/null @@ -1,30 +0,0 @@ -#' 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 deleted file mode 100644 index ea4bb87..0000000 --- a/R/plot.ctryvecest.R +++ /dev/null @@ -1,485 +0,0 @@ -#' 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 deleted file mode 100644 index b9c0d69..0000000 --- a/R/print.summary.bgvecest.R +++ /dev/null @@ -1,15 +0,0 @@ -#' @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 deleted file mode 100644 index b54d28e..0000000 --- a/R/print.summary.ctryvecest.R +++ /dev/null @@ -1,113 +0,0 @@ -#' @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 deleted file mode 100644 index 7ffcb49..0000000 --- a/R/summary.bgvecest.R +++ /dev/null @@ -1,45 +0,0 @@ -#' 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 deleted file mode 100644 index 9ace33d..0000000 --- a/R/summary.ctryvecest.R +++ /dev/null @@ -1,215 +0,0 @@ -#' 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 deleted file mode 100644 index 55dc8c7..0000000 --- a/R/thin.bgvecest.R +++ /dev/null @@ -1,74 +0,0 @@ -#' 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/src/RcppExports.cpp b/src/RcppExports.cpp index 8d541ad..7f506e8 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -22,39 +22,6 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// bgvartvpalg -Rcpp::List bgvartvpalg(Rcpp::List object); -RcppExport SEXP _bgvars_bgvartvpalg(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(bgvartvpalg(object)); - 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) { @@ -118,9 +85,6 @@ 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/bgvartvpalg.cpp b/src/bgvartvpalg.cpp deleted file mode 100644 index 147fd57..0000000 --- a/src/bgvartvpalg.cpp +++ /dev/null @@ -1,737 +0,0 @@ -#include -#include -// [[Rcpp::depends(RcppArmadillo)]] -// [[Rcpp::depends(bvartools)]] - -// [[Rcpp::export(.bgvartvpalg)]] -Rcpp::List bgvartvpalg(Rcpp::List object) { - - // Initialise variables - Rcpp::List data = object["data"]; - arma::mat y = arma::trans(Rcpp::as(data["Y"])); - const arma::mat yvec = arma::vectorise(y); - arma::mat z = Rcpp::as(data["SUR"]); - const int n_tot = z.n_cols; - - // Model information - Rcpp::List model = object["model"]; - Rcpp::CharacterVector model_names = model.names(); - Rcpp::List domestic = model["domestic"]; - Rcpp::CharacterVector endo_names = Rcpp::as(domestic["variables"]); - - // Define useful variables - const int tt = y.n_cols; - const int k_dom = y.n_rows; - const int p_dom = Rcpp::as(domestic["lags"]); - const int n_dom = k_dom * k_dom * p_dom; - 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 = 0; - int n_det = 0; - const int n_sigma = k_dom * k_dom; - const arma::mat diag_k = arma::eye(k_dom, k_dom); - const arma::sp_mat diag_tt = arma::speye(tt, tt); - const arma::vec vec_tt = arma::ones(tt); // T vector - const bool sv = Rcpp::as(model["sv"]); - const bool structural = Rcpp::as(model["structural"]); - int n_a0 = 0; - int n_psi = 0; - 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 + 1); - - // 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 + 1); - } - - // Deterministic terms - Rcpp::CharacterVector det_names; - if (std::find(model_names.begin(), model_names.end(), "deterministic") != model_names.end()) { - det_names = Rcpp::as(model["deterministic"]); - k_det = det_names.length(); - n_det = k_dom * k_det; - } - - // Priors & initial values ---- - Rcpp::List priors = object["priors"]; - Rcpp::CharacterVector priors_names = priors.names(); - Rcpp::List initial = object["initial"]; - - // Priors - Coefficients - Rcpp::List priors_coefficients = priors["coefficients"]; - Rcpp::CharacterVector prcoeff_names = priors_coefficients.names(); - Rcpp::List init_coeffs = initial["coefficients"]; - - Rcpp::List a_prior_varsel; - arma::vec a_init, a_init_post_mu, a_post_mu, a_sigma_v_post_shape, a_sigma_v_prior_rate, a_sigma_v_post_scale; - arma::vec a_bvs_lprior_0, a_bvs_lprior_1, a_bvs_l0_res, a_bvs_l1_res, a_randu, a_prior_incl, a_varsel_include, a_varsel_include_draw; - arma::mat a, a_AG, a_init_prior_mu, a_init_prior_vi, a_init_post_v, a_lag, a_mat, a_post_v, a_theta0, a_theta1, a_v, zz_bvs; - arma::mat a_lambda, a_sigma_v, a_sigma_v_i; - int a_varsel_n, a_varsel_pos; - double a_l0, a_l1, a_bayes, a_bayes_rand; - arma::sp_mat z_large = arma::zeros(tt * k_dom, tt * n_tot); // Final data matrix - for (int i = 0; i < tt; i++) { - z_large.submat(i * k_dom, i * n_tot, (i + 1) * k_dom - 1, (i + 1) * n_tot - 1) = z.rows(i * k_dom, (i + 1) * k_dom - 1); - } - - // Measurement - a_init_prior_mu = Rcpp::as(priors_coefficients["mu"]); - a_init_prior_vi = Rcpp::as(priors_coefficients["v_i"]); - a_init_post_mu = a_init_prior_mu * 0; - // State - a_sigma_v_post_shape = Rcpp::as(priors_coefficients["shape"]) + 0.5 * tt; - a_sigma_v_prior_rate = Rcpp::as(priors_coefficients["rate"]); - - const arma::mat a_b = arma::eye(n_tot, n_tot); - a_post_mu = arma::zeros(n_tot * tt); - a_post_v = arma::zeros(tt * n_tot, tt * n_tot); - a = arma::zeros(n_tot, tt); - a_lag = a; - a_sigma_v = arma::eye(n_tot, n_tot); - a_sigma_v.diag() = a_sigma_v_prior_rate; - a_sigma_v_i = arma::eye(n_tot, n_tot); - a_sigma_v_i.diag() = 1 / a_sigma_v_prior_rate; - - a_init = Rcpp::as(init_coeffs["draw"]); - - // Priors - Coefficients - BVS - if (std::find(prcoeff_names.begin(), prcoeff_names.end(), "bvs") != prcoeff_names.end()) { - bvs = true; - a_prior_varsel = priors_coefficients["bvs"]; - a_bvs_lprior_0 = arma::log(1 - Rcpp::as(a_prior_varsel["inprior"])); - a_bvs_lprior_1 = arma::log(Rcpp::as(a_prior_varsel["inprior"])); - a_varsel_include = Rcpp::as(a_prior_varsel["include"]) - 1; - a_varsel_n = size(a_varsel_include)(0); - a_lambda = arma::eye(n_tot, n_tot); - a_l0 = 0; - a_l1 = 0; - a_bayes = 0; - a_bayes_rand = 0; - zz_bvs = z; - } - - // Priors - Covar coefficients - 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::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, psi_init, psi_init_post_mu, psi_lag, 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_Psi, psi_AG, psi_init_post_v, psi_mat, psi_post_v, psi_theta0, psi_theta1, psi_v, psi_z_bvs; - arma::sp_mat diag_covar_omega_i, Psi, psi_hh, psi_h_sigmav_i_h, psi_lambda, psi_sigma_v_i, psi_z, psi_zzss_i; - if (covar) { - n_psi = k_dom * (k_dom - 1) / 2; - Psi = arma::speye(k_dom * tt, k_dom * tt); - Rcpp::List init_psi = initial["psi"]; - psi_init = Rcpp::as(init_psi["draw"]); - psi_lag = arma::zeros(n_psi * tt); - psi_z = arma::zeros((k_dom - 1) * tt, n_psi * tt); - psi_hh = arma::speye(n_psi * tt, n_psi * tt); - psi_hh.diag(-n_psi) = -arma::ones(n_psi * (tt - 1)); - diag_covar_omega_i = arma::zeros(tt * (k_dom - 1), tt * (k_dom - 1)); - psi_sigma_v_i = arma::speye(n_psi, n_psi); - psi_sigma_v_i.diag() = 1 / psi_sigma_v_prior_rate; - if (psi_bvs) { - psi_varsel_include = Rcpp::as(psi_prior_varsel["include"]) - 1; - psi_varsel_n = size(psi_varsel_include)(0); - psi_lambda = arma::speye(n_psi, n_psi); - psi_l0 = 0; - psi_l1 = 0; - 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, sse, sigma_u; - arma::mat u = y * 0; - arma::mat sigma_u_i, omega_i; - arma::sp_mat diag_sigma_u_i = arma::zeros(k_dom * tt, k_dom * tt); - arma::sp_mat diag_omega_i = arma::zeros(k_dom * tt, k_dom * tt); - if (covar || sv) { - sigma_u = 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)); - sigma_u_i = arma::diagmat(1 / exp(h_init)); - for (int i = 0; i < tt; i++) { - sigma_u.rows(i * k_dom, (i + 1) * k_dom - 1) = arma::solve(sigma_u_i, diag_k); - } - h_const = Rcpp::as(init_sigma["constant"]); - } else { - if (use_gamma) { - omega_i = arma::zeros(k_dom, k_dom); - omega_i.diag() = Rcpp::as(init_sigma["sigma_i"]).diag(); - diag_omega_i.diag() = arma::repmat(omega_i.diag(), tt, 1); - if (covar) { - for (int i = 0; i < tt; i++) { - sigma_u.rows(i * k_dom, (i + 1) * k_dom - 1) = arma::solve(omega_i, diag_k); - } - } - } else { - omega_i = Rcpp::as(init_sigma["sigma_i"]); - sigma_u = arma::solve(omega_i, diag_k); - } - sigma_u_i = omega_i; - } - diag_sigma_u_i.diag() = arma::repmat(sigma_u_i.diag(), tt, 1); - if (covar || sv) { - diag_omega_i = diag_sigma_u_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 = 0; - const int dom_pos_start = 0; - const int dom_pos_end = n_dom - 1; - const int for_pos_start = n_dom; - const int for_pos_end = n_dom + n_for - 1; - const int glo_pos_start = n_dom + n_for; - const int glo_pos_end = n_dom + n_for + n_glo - 1; - const int det_pos_start = n_dom + n_for + n_glo; - const int det_pos_end = n_dom + n_for + n_glo + n_det - 1; - const int a0_pos_start = n_dom + n_for + n_glo + n_det; - const int a0_pos_end = n_dom + n_for + n_glo + n_det + n_a0 - 1; - - arma::mat draws_a0 = arma::zeros(n_a0 * tt, iter); - arma::mat draws_sigma_a0 = arma::zeros(n_a0, iter); - arma::mat draws_dom = arma::zeros(n_dom * tt, iter); - arma::mat draws_sigma_dom = arma::zeros(n_dom, iter); - arma::mat draws_for = arma::zeros(n_for * tt, iter); - arma::mat draws_sigma_for = arma::zeros(n_for, iter); - arma::mat draws_glo = arma::zeros(n_glo * tt, iter); - arma::mat draws_sigma_glo = arma::zeros(n_glo, iter); - arma::mat draws_det = arma::zeros(n_det * tt, iter); - arma::mat draws_sigma_det = arma::zeros(n_det, iter); - arma::mat draws_sigma_u, draws_sigma_sigma; - if (sv || covar) { - draws_sigma_u = arma::zeros(k_dom * k_dom * tt, iter); - if (sv) { - draws_sigma_sigma = arma::zeros(k_dom * k_dom, iter); - } - } else { - draws_sigma_u = 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_det; - if (bvs) { - 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_det = arma::zeros(n_det, iter); - } - if (covar && psi_bvs) { - 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 interruption every now and then - Rcpp::checkUserInterrupt(); - } - - // Draw a - if (bvs) { - z = zz_bvs * a_lambda; - } - a = bvartools::kalman_dk(y, z, sigma_u, a_sigma_v, a_b, a_init, a_sigma_v); - a = a.cols(1, tt); - - // Draw sigma_v_i - a_lag.col(0) = a_init; - a_lag.cols(1, tt - 1) = a.cols(0, tt - 2); - a_v = arma::trans(a - a_lag); - a_sigma_v_post_scale = 1 / (a_sigma_v_prior_rate + arma::vectorise(arma::sum(arma::pow(a_v, 2))) * 0.5); - for (int i = 0; i < n_tot; i++) { - a_sigma_v_i(i, i) = arma::randg(arma::distr_param(a_sigma_v_post_shape(i), a_sigma_v_post_scale(i))); - } - - // Draw initial state of a - a_init_post_v = a_init_prior_vi + a_sigma_v_i; - a_init_post_mu = arma::solve(a_init_post_v, a_init_prior_vi * a_init_prior_mu + a_sigma_v_i * a.col(0)); - a_init = a_init_post_mu + arma::solve(arma::chol(a_init_post_v), arma::randn(n_tot)); - - // BVS - if (bvs) { - z = zz_bvs; - a_mat = a; - a_AG = a_lambda * a_mat; - a_varsel_include_draw = shuffle(a_varsel_include); // Reorder positions of variable selection - 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) = arma::zeros(1, tt); - a_theta1.row(a_varsel_pos) = a_mat.row(a_varsel_pos); - a_bvs_l0_res = yvec - z_large * arma::vectorise(a_theta0); - a_bvs_l1_res = yvec - z_large * arma::vectorise(a_theta1); - a_l0 = -arma::as_scalar(trans(a_bvs_l0_res) * diag_sigma_u_i * a_bvs_l0_res) * 0.5 + arma::as_scalar(a_bvs_lprior_0(a_varsel_pos)); - a_l1 = -arma::as_scalar(trans(a_bvs_l1_res) * diag_sigma_u_i * a_bvs_l1_res) * 0.5 + 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; - } - } - } - a = a_lambda * a_mat; - lambda_vec = a_lambda.diag(); - } - - for (int i = 0; i < tt; i++) { - u.col(i) = y.col(i) - z.rows(i * k_dom, (i + 1) * k_dom - 1) * a.col(i); - } - - // Covariances - if (covar) { - - // Prepare data - u_vec = arma::vectorise(u); - 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, - j * n_psi + i * (i - 1) / 2, - j * (k_dom - 1) + i - 1, - j * n_psi + (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 * arma::kron(diag_tt, psi_lambda); - } - - psi_zzss_i = arma::trans(psi_z) * diag_covar_omega_i; - psi_h_sigmav_i_h = arma::trans(psi_hh) * arma::kron(diag_tt, psi_sigma_v_i) * psi_hh; - psi_post_v = psi_h_sigmav_i_h + psi_zzss_i * psi_z; - psi_post_mu = arma::solve(psi_post_v, psi_h_sigmav_i_h * arma::kron(vec_tt, psi_init) + psi_zzss_i * psi_y); - psi = psi_post_mu + arma::solve(arma::chol(psi_post_v), arma::randn(n_psi * tt)); - - if (psi_bvs) { - - // Reorder positions of variable selection - psi_varsel_include_draw = shuffle(psi_varsel_include); - - psi_z = psi_z_bvs; - psi_mat = arma::reshape(psi, n_psi, tt); - psi_AG = psi_lambda * psi_mat; - 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 = arma::vectorise(psi_lambda * psi_mat); - psi_lambda_vec = psi_lambda.diag(); - } - - for (int j = 0; j < tt; j ++) { - for (int i = 1; i < k_dom; i++) { - Psi.submat((k_dom * j) + i, k_dom * j, (k_dom * j) + i, (k_dom * j) + i - 1) = arma::trans(psi.subvec((n_psi * j) + i * (i - 1) / 2, (n_psi * j) + (i + 1) * i / 2 - 1)); - } - } - - u = arma::reshape(Psi * u_vec, k_dom, tt); - - // Draw sigma_v_i - psi_lag.subvec(0, n_psi - 1) = psi_init; - psi_lag.subvec(n_psi, tt * n_psi - 1) = psi.subvec(0, (tt - 1) * n_psi - 1); - psi_v = arma::trans(arma::reshape(psi - psi_lag, n_psi, tt)); - 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))); - } - - // 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.subvec(0, n_psi - 1)); - psi_init = psi_init_post_mu + arma::solve(arma::chol(psi_init_post_v), arma::randn(n_psi)); - } - - /////////////////////////////////////////////////////////////////////// - // 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 diag_omega_i - diag_omega_i.diag() = 1 / exp(arma::vectorise(h.t())); - // Update diag_sigma_u_i - if (covar) { - diag_sigma_u_i = arma::trans(Psi) * diag_omega_i * Psi; - } else { - diag_sigma_u_i = diag_omega_i; - } - // Update sigma_u - for (int i = 0; i < tt; i++) { - sigma_u.rows(i * k_dom, (i + 1) * k_dom - 1) = arma::solve(arma::mat(diag_sigma_u_i.submat(i * k_dom, i * k_dom, (i + 1) * k_dom - 1, (i + 1) * k_dom - 1)), diag_k); - } - - // 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))); - } - // Repeat output for large (sparse) diagonal matrix - diag_omega_i.diag() = arma::repmat(omega_i.diag(), tt, 1); - - if (covar) { - diag_sigma_u_i = arma::trans(Psi) * diag_omega_i * Psi; - for (int i = 0; i < tt; i++) { - sigma_u.rows(i * k_dom, (i + 1) * k_dom - 1) = arma::solve(arma::mat(diag_sigma_u_i.submat(i * k_dom, i * k_dom, (i + 1) * k_dom - 1, (i + 1) * k_dom - 1)), diag_k); - } - } else { - sigma_u_i = omega_i; - diag_sigma_u_i = diag_omega_i; - sigma_u = arma::solve(sigma_u_i, diag_k); - } - - } else { - sigma_u_i = arma::wishrnd(arma::solve(sigma_prior_scale + sse, diag_k), sigma_post_df); - sigma_u = arma::solve(sigma_u_i, diag_k); - if (bvs) { - diag_sigma_u_i = arma::kron(diag_tt, arma::sp_mat(sigma_u_i)); - } - } - } - - // Store draws - if (draw >= burnin) { - - pos_draw = draw - burnin; - - if (sv || covar) { - 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(arma::solve(arma::mat(diag_sigma_u_i.submat(i * k_dom, i * k_dom, (i + 1) * k_dom - 1, (i + 1) * k_dom - 1)), diag_k)); - } - if (sv) { - draws_sigma_sigma.col(pos_draw) = arma::vectorise(arma::diagmat(sigma_h)); - } - } else { - draws_sigma_u.col(pos_draw) = arma::vectorise(sigma_u); - } - - if (psi_bvs) { - draws_lambda_a0.col(pos_draw) = psi_lambda_vec; - } - - a_mat = a; - - if (n_dom > 0) { - draws_dom.col(pos_draw) = arma::vectorise(a_mat.rows(dom_pos_start, dom_pos_end)); - draws_sigma_dom.col(pos_draw) = 1 / arma::mat(a_sigma_v_i.submat(dom_pos_start, dom_pos_start, dom_pos_end, dom_pos_end)).diag(); - if (bvs) { - 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(a_mat.rows(for_pos_start, for_pos_end)); - draws_sigma_for.col(pos_draw) = 1 / arma::mat(a_sigma_v_i.submat(for_pos_start, for_pos_start, for_pos_end, for_pos_end)).diag(); - if (bvs) { - 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(a_mat.rows(glo_pos_start, glo_pos_end)); - draws_sigma_glo.col(pos_draw) = 1 / arma::mat(a_sigma_v_i.submat(glo_pos_start, glo_pos_start, glo_pos_end, glo_pos_end)).diag(); - if (bvs) { - draws_lambda_glo.col(pos_draw) = lambda_vec.subvec(glo_pos_start, glo_pos_end); - } - } - if (n_det > 0) { - draws_det.col(pos_draw) = arma::vectorise(a_mat.rows(det_pos_start, det_pos_end)); - draws_sigma_det.col(pos_draw) = 1 / arma::mat(a_sigma_v_i.submat(det_pos_start, det_pos_start, det_pos_end, det_pos_end)).diag(); - if (bvs) { - draws_lambda_det.col(pos_draw) = lambda_vec.subvec(det_pos_start, det_pos_end); - } - } - if (structural) { - draws_a0.col(pos_draw) = arma::vectorise(a_mat.rows(a0_pos_start, a0_pos_end)); - draws_sigma_a0.col(pos_draw) = 1 / arma::mat(a_sigma_v_i.submat(a0_pos_start, a0_pos_start, a0_pos_end, a0_pos_end)).diag(); - if (bvs) { - draws_lambda_a0.col(pos_draw) = lambda_vec.subvec(a0_pos_start, a0_pos_end); - } - } - } // End storage condition - } // End loop - - Rcpp::List posteriors = Rcpp::List::create(Rcpp::Named("a0") = R_NilValue, - Rcpp::Named("domestic") = R_NilValue, - Rcpp::Named("foreign") = R_NilValue, - Rcpp::Named("global") = R_NilValue, - Rcpp::Named("deterministic") = R_NilValue, - Rcpp::Named("sigma") = R_NilValue); - - if (n_dom > 0) { - if (bvs) { - posteriors["domestic"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_dom, - Rcpp::Named("sigma") = draws_sigma_dom, - Rcpp::Named("lambda") = draws_lambda_dom)); - } else { - posteriors["domestic"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_dom, - Rcpp::Named("sigma") = draws_sigma_dom)); - } - } - - if (n_for > 0) { - if (bvs) { - posteriors["foreign"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_for, - Rcpp::Named("sigma") = draws_sigma_for, - Rcpp::Named("lambda") = draws_lambda_for)); - } else { - posteriors["foreign"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_for, - Rcpp::Named("sigma") = draws_sigma_for)); - } - } - - if (n_glo > 0) { - if (bvs) { - posteriors["global"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_glo, - Rcpp::Named("sigma") = draws_sigma_glo, - Rcpp::Named("lambda") = draws_lambda_glo)); - } else { - posteriors["global"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_glo, - Rcpp::Named("sigma") = draws_sigma_glo)); - } - } - - if (n_det > 0) { - if (bvs) { - posteriors["deterministic"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_det, - Rcpp::Named("sigma") = draws_sigma_det, - Rcpp::Named("lambda") = draws_lambda_det)); - } else { - posteriors["deterministic"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_det, - Rcpp::Named("sigma") = draws_sigma_det)); - } - } - - if (structural) { - if (bvs) { - posteriors["a0"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_a0, - Rcpp::Named("sigma") = draws_sigma_a0, - Rcpp::Named("lambda") = draws_lambda_a0)); - } else { - posteriors["a0"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_a0, - Rcpp::Named("sigma") = draws_sigma_a0)); - } - } - - 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_lambda_a0)); - } else { - posteriors["sigma"] = Rcpp::wrap(Rcpp::List::create(Rcpp::Named("coeffs") = draws_sigma_u, - Rcpp::Named("lambda") = draws_lambda_a0)); - } - } 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); - - - // return Rcpp::List::create(Rcpp::Named("test") = sigma_i); -} - - -/*** R - -#library(bgvars) -library(Matrix) - -# 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 - -# Make series stationary -country_data <- diff_variables(country_data, variables = c("y", "Dp", "r"), multi = 100) -global_data <- diff_variables(global_data, multi = 100) - -# 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" = TRUE), - countries = c("EA", "US", "GB", "CA", "JP", "IN"), - type = "VAR", - tvp = TRUE, sv = TRUE, - iterations = 10, - burnin = 10) - -# Country-specific specifications -model_specs[["US"]][["domestic"]][["variables"]] <- c("poil", "r", "Dp", "y") -model_specs[["US"]][["foreign"]][["variables"]] <- c("y", "Dp") - -# 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, shape = 3, rate = .0001, rate_det = .01), - #bvs = list(inprior = .5), - sigma = list(shape = 3, rate = .0001, mu = 0, v_i = .1, hinit = 0.05, constant = .0001, covar = TRUE)) - -.bgvartvpalg(temp_model[[1]]) - -***/ \ No newline at end of file diff --git a/src/bgvecalg.cpp b/src/bgvecalg.cpp deleted file mode 100644 index 9895b60..0000000 --- a/src/bgvecalg.cpp +++ /dev/null @@ -1,890 +0,0 @@ -#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 deleted file mode 100644 index b58e2a1..0000000 --- a/src/bgvectvpalg.cpp +++ /dev/null @@ -1,985 +0,0 @@ -#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 - -*/