This vignette provides the code to set up and estimate a basic 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.
To illustrate the estimation process, the dataset E6 from Lütkepohl (2007) is used, which contains data on German long-term interest rates and inflation from 1972Q2 to 1998Q4.
The gen_vec function produces the inputs 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.
# Reset random number generator for reproducibility
set.seed(7654321)
iter <- 10000 # Number of iterations of the Gibbs sampler
burnin <- 5000 # Number of burn-in draws
store <- iter - burnin
r <- 1 # 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_alpha <- k * r # Number of elements in alpha
k_beta <- k_w * r # Number of elements in beta
k_gamma <- k * k_x
# Set uninformative priors
a_mu_prior <- matrix(0, k_x * k) # Vector of prior parameter means
a_v_i_prior <- diag(0, k_x * k) # Inverse of the prior covariance matrix
v_i <- 0
p_tau_i <- matrix(0, k_w, k_w)
p_tau_i[1:r, 1:r] <- diag(1, r)
sigma_df_prior <- k + r # Prior degrees of freedom
sigma_scale_prior <- diag(0, k) # 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(.00001, k)
sigma <- solve(sigma_i)
g_i <- sigma_i
# Data containers
draws_alpha <- matrix(NA, k_alpha, store)
draws_beta <- matrix(NA, k_beta, store)
draws_pi <- matrix(NA, k * k_w, store)
draws_gamma <- matrix(NA, k_gamma, store)
draws_sigma <- matrix(NA, k^2, store)
# Start Gibbs sampler
for (draw in 1:iter) {
# 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
#> [,1]
#> l.R 1.000
#> l.Dp -3.966bvec 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 = y, w = w, x = x,
Pi = draws_pi,
Gamma = draws_gamma[1:k_nondet,],
C = draws_gamma[(k_nondet + 1):nrow(draws_gamma),],
Sigma = draws_sigma)Obtain summaries of posterior draws
summary(bvec_est)
#>
#> Model:
#>
#> y ~ l.R + l.Dp + d.R.1 + d.Dp.1 + d.R.2 + d.Dp.2 + d.R.3 + d.Dp.3 + const + season.1 + season.2 + season.3
#>
#> Variable: d.R
#>
#> Mean SD Naive SD Time-series SD 2.5% 50%
#> l.R -0.1027368 0.055975 7.916e-04 2.787e-03 -0.215454 -0.1013532
#> l.Dp 0.3777113 0.189482 2.680e-03 5.073e-03 0.001146 0.3772081
#> d.R.1 0.2670722 0.106362 1.504e-03 1.955e-03 0.057599 0.2664760
#> d.Dp.1 -0.1876540 0.158783 2.246e-03 3.780e-03 -0.495507 -0.1876838
#> d.R.2 -0.0168650 0.106392 1.505e-03 1.655e-03 -0.229040 -0.0159563
#> d.Dp.2 -0.2073960 0.128465 1.817e-03 2.886e-03 -0.460733 -0.2058733
#> d.R.3 0.2245038 0.105663 1.494e-03 1.950e-03 0.015177 0.2243238
#> d.Dp.3 -0.0991798 0.085253 1.206e-03 1.571e-03 -0.266135 -0.1001086
#> const 0.0017789 0.004308 6.092e-05 1.791e-04 -0.006395 0.0016503
#> season.1 0.0014277 0.005070 7.171e-05 7.171e-05 -0.008566 0.0014706
#> season.2 0.0088265 0.005208 7.365e-05 7.627e-05 -0.001403 0.0088886
#> season.3 -0.0004358 0.005001 7.072e-05 7.072e-05 -0.010238 -0.0004517
#> 97.5%
#> l.R -0.000283
#> l.Dp 0.747755
#> d.R.1 0.476899
#> d.Dp.1 0.130833
#> d.R.2 0.190115
#> d.Dp.2 0.039490
#> d.R.3 0.429270
#> d.Dp.3 0.070469
#> const 0.010544
#> season.1 0.011484
#> season.2 0.019000
#> season.3 0.009347
#>
#> Variable: d.Dp
#>
#> Mean SD Naive SD Time-series SD 2.5% 50%
#> l.R 0.1456078 0.047596 6.731e-04 1.249e-03 0.048922 0.14561
#> l.Dp -0.5670470 0.197330 2.791e-03 8.417e-03 -0.956057 -0.57178
#> d.R.1 0.0742445 0.100122 1.416e-03 1.416e-03 -0.123126 0.07283
#> d.Dp.1 -0.3818973 0.164995 2.333e-03 6.542e-03 -0.707916 -0.37996
#> d.R.2 -0.0005796 0.102168 1.445e-03 1.445e-03 -0.201152 -0.00139
#> d.Dp.2 -0.4202953 0.128535 1.818e-03 4.496e-03 -0.672229 -0.42055
#> d.R.3 0.0248960 0.097099 1.373e-03 1.548e-03 -0.161884 0.02355
#> d.Dp.3 -0.3604562 0.082801 1.171e-03 2.127e-03 -0.521991 -0.36060
#> const 0.0106490 0.003971 5.616e-05 1.233e-04 0.002763 0.01072
#> season.1 -0.0342656 0.004838 6.842e-05 6.842e-05 -0.043825 -0.03427
#> season.2 -0.0179348 0.005053 7.145e-05 7.145e-05 -0.027908 -0.01792
#> season.3 -0.0165758 0.004761 6.734e-05 6.558e-05 -0.025843 -0.01652
#> 97.5%
#> l.R 0.238231
#> l.Dp -0.179161
#> d.R.1 0.276389
#> d.Dp.1 -0.058483
#> d.R.2 0.199223
#> d.Dp.2 -0.176216
#> d.R.3 0.214603
#> d.Dp.3 -0.200261
#> const 0.018226
#> season.1 -0.024735
#> season.2 -0.008060
#> season.3 -0.007116
#>
#> Variance-covariance matrix:
#>
#> Mean SD Naive SD Time-series SD 2.5% 50%
#> d.R_d.R 2.878e-05 4.272e-06 6.041e-08 7.578e-08 2.166e-05 2.838e-05
#> d.R_d.Dp -1.874e-06 2.842e-06 4.019e-08 4.334e-08 -7.497e-06 -1.860e-06
#> d.Dp_d.R -1.874e-06 2.842e-06 4.019e-08 4.334e-08 -7.497e-06 -1.860e-06
#> d.Dp_d.Dp 2.602e-05 3.923e-06 5.548e-08 6.670e-08 1.931e-05 2.561e-05
#> 97.5%
#> d.R_d.R 3.844e-05
#> d.R_d.Dp 3.505e-06
#> d.Dp_d.R 3.505e-06
#> d.Dp_d.Dp 3.476e-05Posterior draws can be thinned with function thin:
The function bvec_to_bvar can be used to transform the VEC model into a VAR in levels:
Impulse responses for VECs can be constructed from their VAR respresentations.
IR <- irf(bvar_form, impulse = "R", response = "Dp", n.ahead = 20)
plot(IR, main = "Forecast Error Impulse Response", xlab = "Year", ylab = "Response")
OIR <- irf(bvar_form, impulse = "R", response = "Dp", n.ahead = 20, type = "oir")
plot(OIR, main = "Orthogonalised Impulse Response", xlab = "Year", ylab = "Response")
GIR <- irf(bvar_form, impulse = "R", response = "Dp", n.ahead = 20, type = "gir")
plot(GIR, main = "Generalised Impulse Response", xlab = "Year", ylab = "Response")
bvec_fevd <- fevd(bvar_form, response = "Dp", n.ahead = 20)
plot(bvec_fevd, main = "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. (2007). 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