knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
This vignette provides the code to set up and estimate a Bayesian vector error correction (BVEC) model with the bvartools
package. The presented Gibbs sampler is based on the approach of Koop et al. (2010), who propose a prior on the cointegration space. The estimated model has the following form
$$\Delta y_t = \Pi y_{t - 1} + \sum_{l = 1}^{p - 1} \Gamma_l \Delta y_{t - l} + C d_t + u_t,$$
where $\Pi = \alpha \beta^{\prime}$ with cointegration rank $r$, $u_t \sim N(0, \Sigma)$ and $d_t$ contains an intercept and seasonal dummies. For an introduction to vector error correction models see https://www.r-econometrics.com/timeseries/vecintro/.
To illustrate the workflow of the analysis, data set E6 from Lütkepohl (2006) is used, which contains data on German long-term interest rates and inflation from 1972Q2 to 1998Q4.
library(bvartools) data("e6") plot(e6) # Plot the series
The gen_vec
function produces the data matrices Y
, W
and X
for the BVEC estimator, where Y
is the matrix of dependent variables, W
is a matrix of potentially cointegrated regressors, and X
is the matrix of non-cointegration regressors.
data <- gen_vec(e6, p = 4, r = 1, const = "unrestricted", season = "unrestricted", iterations = 5000, burnin = 1000)
Argument p
represents the lag order of the VAR form of the model and r
is the cointegration rank of Pi
. Function gen_vec
requires to specify the inclusion of intercepts, trends and seasonal dummies separately. This allows to decide on whether they enter the cointegration term or the non-cointegration part of the model.
Function add_priors
adds the necessary prior specifications to object data
. We
For the current application we use non-informative priors.
data <- add_priors(data, coint = list(v_i = 0, p_tau_i = 1), coef = list(v_i = 0, v_i_det = 0), sigma = list(df = 0, scale = .0001))
The following code produces posterior draws using the algortihm described in Koop et al. (2010).
# Reset random number generator for reproducibility set.seed(7654321) # Obtain data matrices y <- t(data$data$Y) w <- t(data$data$W) x <- t(data$data$X) r <- data$model$rank # Set rank tt <- ncol(y) # Number of observations k <- nrow(y) # Number of endogenous variables k_w <- nrow(w) # Number of regressors in error correction term k_x <- nrow(x) # Number of differenced regressors and unrestrictec deterministic terms k_gamma <- k * k_x # Total number of non-cointegration coefficients k_alpha <- k * r # Number of elements in alpha k_beta <- k_w * r # Number of elements in beta # Priors a_mu_prior <- data$priors$noncointegration$mu # Prior means a_v_i_prior <- data$priors$noncointegration$v_i # Inverse of the prior covariance matrix v_i <- data$priors$cointegration$v_i p_tau_i <- data$priors$cointegration$p_tau_i sigma_df_prior <- data$priors$sigma$df # Prior degrees of freedom sigma_scale_prior <- data$priors$sigma$scale # Prior covariance matrix sigma_df_post <- tt + sigma_df_prior # Posterior degrees of freedom # Initial values beta <- matrix(0, k_w, r) beta[1:r, 1:r] <- diag(1, r) sigma_i <- diag(1 / .0001, k) g_i <- sigma_i iterations <- data$model$iterations # Number of iterations of the Gibbs sampler burnin <- data$model$burnin # Number of burn-in draws draws <- iterations + burnin # Total number of draws # Data containers draws_alpha <- matrix(NA, k_alpha, iterations) draws_beta <- matrix(NA, k_beta, iterations) draws_pi <- matrix(NA, k * k_w, iterations) draws_gamma <- matrix(NA, k_gamma, iterations) draws_sigma <- matrix(NA, k^2, iterations) # Start Gibbs sampler for (draw in 1:draws) { # Draw conditional mean parameters temp <- post_coint_kls(y = y, beta = beta, w = w, x = x, sigma_i = sigma_i, v_i = v_i, p_tau_i = p_tau_i, g_i = g_i, gamma_mu_prior = a_mu_prior, gamma_v_i_prior = a_v_i_prior) alpha <- temp$alpha beta <- temp$beta Pi <- temp$Pi gamma <- temp$Gamma # Draw variance-covariance matrix u <- y - Pi %*% w - matrix(gamma, k) %*% x sigma_scale_post <- solve(tcrossprod(u) + v_i * alpha %*% tcrossprod(crossprod(beta, p_tau_i) %*% beta, alpha)) sigma_i <- matrix(rWishart(1, sigma_df_post, sigma_scale_post)[,, 1], k) sigma <- solve(sigma_i) # Update g_i g_i <- sigma_i # Store draws if (draw > burnin) { draws_alpha[, draw - burnin] <- alpha draws_beta[, draw - burnin] <- beta draws_pi[, draw - burnin] <- Pi draws_gamma[, draw - burnin] <- gamma draws_sigma[, draw - burnin] <- sigma } }
Obtain point estimates of cointegration variables:
beta <- apply(t(draws_beta) / t(draws_beta)[, 1], 2, mean) # Obtain means for every row beta <- matrix(beta, k_w) # Transform mean vector into a matrix beta <- round(beta, 3) # Round values dimnames(beta) <- list(dimnames(w)[[1]], NULL) # Rename matrix dimensions beta # Print
bvec
objectsThe bvec
function can be used to collect output of the Gibbs sampler in a standardised object, which can be used further for forecasting, impulse response analysis or forecast error variance decomposition.
# Number of non-deterministic coefficients k_nondet <- (k_x - 4) * k # Generate bvec object bvec_est <- bvec(y = data$data$Y, w = data$data$W, x = data$data$X[, 1:6], x_d = data$data$X[, -(1:6)], Pi = draws_pi, r = 1, Gamma = draws_gamma[1:k_nondet,], C = draws_gamma[(k_nondet + 1):nrow(draws_gamma),], Sigma = draws_sigma)
Posterior draws an be inspected with plot
:
plot(bvec_est)
Obtain summaries of posterior draws
summary(bvec_est)
As an alternative to a user-specific algorithm function draw_posterior
can be used estimate such a model as well:
bvec_est <- draw_posterior(data)
Posterior draws can be thinned with function thin
:
bvec_est <- thin(bvec_est, thin = 5)
The function bvec_to_bvar
can be used to transform the VEC model into a VAR in levels:
bvar_form <- bvec_to_bvar(bvec_est)
The output of bvec_to_bvar
is an object of class bvar
for which summary statistics, forecasts, impulse responses and variance decompositions can be obtained in the usual manner using summary
, predict
, irf
and fevd
, respectively.
summary(bvar_form)
Forecasts with credible bands can be obtained with the function predict
. If the model contains deterministic terms, new values can be provided in the argument new_d
. If no values are provided, the function sets them to zero. For the current model, seasonal dummies need to be provided. They are taken from the original series. The number of rows of new_d
must be the same as the argument n.ahead
.
# Generate deterministc terms for function predict new_d <- data$data$X[3 + 1:10, c("const", "season.1", "season.2", "season.3")] # Genrate forecasts bvar_pred <- predict(bvar_form, n.ahead = 10, new_d = new_d) # Plot forecasts plot(bvar_pred)
FEIR <- irf(bvar_form, impulse = "R", response = "Dp", n.ahead = 20) plot(FEIR, main = "Forecast Error Impulse Response", xlab = "Period", ylab = "Response")
bvar_fevd_oir <- fevd(bvar_form, response = "Dp", n.ahead = 20) plot(bvar_fevd_oir, main = "OIR-based FEVD of inflation")
Koop, G., León-González, R., & Strachan R. W. (2010). Efficient posterior simulation for cointegrated models with priors on the cointegration space. Econometric Reviews, 29(2), 224-242. https://doi.org/10.1080/07474930903382208
Lütkepohl, H. (2006). New introduction to multiple time series analysis (2nd ed.). Berlin: Springer.
Pesaran, H. H., & Shin, Y. (1998). Generalized impulse response analysis in linear multivariate models, Economics Letters, 58, 17-29. https://doi.org/10.1016/S0165-1765(97)00214-0
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.