The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.
Bayesian Error Correction Models with
Priors on the Cointegration Space
Franz X. Mohr
2024-01-08
Introduction
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/.
Data
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.
Priors
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))
Estimation
User-specific algorithm
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
#> [,1]
#> l.R 1.000
#> l.Dp -3.946
bvec objects
The 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:
Obtain summaries of posterior draws
summary(bvec_est)
#>
#> Bayesian VEC model with p = 4
#>
#> Model:
#>
#> d.y ~ l.R + l.Dp + d.R.l01 + d.Dp.l01 + d.R.l02 + d.Dp.l02 + d.R.l03 + d.Dp.l03 + const + season.1 + season.2 + season.3
#>
#> Variable: d.R
#>
#> Mean SD Naive SD Time-series SD 2.5% 50%
#> l.R -0.1058257 0.057128 8.079e-04 2.934e-03 -0.2218739 -0.104748
#> l.Dp 0.3816631 0.190012 2.687e-03 5.035e-03 -0.0008244 0.381684
#> d.R.l01 0.2686749 0.108454 1.534e-03 1.818e-03 0.0503097 0.268921
#> d.Dp.l01 -0.1902379 0.159116 2.250e-03 3.841e-03 -0.5031594 -0.193049
#> d.R.l02 -0.0166691 0.109222 1.545e-03 1.798e-03 -0.2292546 -0.016324
#> d.Dp.l02 -0.2098870 0.130365 1.844e-03 2.559e-03 -0.4693754 -0.209137
#> d.R.l03 0.2259926 0.105552 1.493e-03 1.924e-03 0.0200153 0.228738
#> d.Dp.l03 -0.1024420 0.085916 1.215e-03 1.215e-03 -0.2734722 -0.102478
#> const 0.0018718 0.004447 6.290e-05 1.847e-04 -0.0065256 0.001778
#> season.1 0.0015565 0.005165 7.304e-05 7.059e-05 -0.0082831 0.001627
#> season.2 0.0089724 0.005277 7.463e-05 7.124e-05 -0.0012790 0.008950
#> season.3 -0.0003219 0.005168 7.308e-05 7.308e-05 -0.0103220 -0.000263
#> 97.5%
#> l.R 0.0002007
#> l.Dp 0.7544752
#> d.R.l01 0.4799672
#> d.Dp.l01 0.1260587
#> d.R.l02 0.2016127
#> d.Dp.l02 0.0485518
#> d.R.l03 0.4331961
#> d.Dp.l03 0.0645323
#> const 0.0106365
#> season.1 0.0116696
#> season.2 0.0191547
#> season.3 0.0097738
#>
#> Variable: d.Dp
#>
#> Mean SD Naive SD Time-series SD 2.5% 50%
#> l.R 0.145090 0.049654 7.022e-04 1.514e-03 0.046217 0.1449742
#> l.Dp -0.556494 0.199952 2.828e-03 8.912e-03 -0.942351 -0.5581037
#> d.R.l01 0.076859 0.104437 1.477e-03 1.765e-03 -0.124850 0.0749406
#> d.Dp.l01 -0.390568 0.166755 2.358e-03 6.892e-03 -0.720068 -0.3885522
#> d.R.l02 0.001176 0.104539 1.478e-03 1.473e-03 -0.203920 0.0007347
#> d.Dp.l02 -0.423783 0.129145 1.826e-03 4.618e-03 -0.675218 -0.4239844
#> d.R.l03 0.024011 0.099281 1.404e-03 1.404e-03 -0.168761 0.0232362
#> d.Dp.l03 -0.362624 0.084262 1.192e-03 2.395e-03 -0.528992 -0.3620497
#> const 0.010553 0.004022 5.688e-05 1.308e-04 0.002849 0.0105357
#> season.1 -0.034204 0.004885 6.909e-05 6.781e-05 -0.043993 -0.0341443
#> season.2 -0.018008 0.004999 7.070e-05 7.070e-05 -0.027921 -0.0179772
#> season.3 -0.016431 0.004864 6.878e-05 6.878e-05 -0.025821 -0.0164508
#> 97.5%
#> l.R 0.242481
#> l.Dp -0.154874
#> d.R.l01 0.289202
#> d.Dp.l01 -0.062096
#> d.R.l02 0.210974
#> d.Dp.l02 -0.173093
#> d.R.l03 0.220960
#> d.Dp.l03 -0.198335
#> const 0.018401
#> season.1 -0.024807
#> season.2 -0.008508
#> season.3 -0.006929
#>
#> Variance-covariance matrix:
#>
#> Mean SD Naive SD Time-series SD 2.5% 50%
#> d.R_d.R 2.948e-05 4.562e-06 6.452e-08 7.135e-08 2.178e-05 2.905e-05
#> d.R_d.Dp -1.868e-06 3.017e-06 4.266e-08 4.780e-08 -7.811e-06 -1.856e-06
#> d.Dp_d.R -1.868e-06 3.017e-06 4.266e-08 4.780e-08 -7.811e-06 -1.856e-06
#> d.Dp_d.Dp 2.671e-05 4.065e-06 5.749e-08 7.304e-08 1.980e-05 2.640e-05
#> 97.5%
#> d.R_d.R 3.963e-05
#> d.R_d.Dp 3.891e-06
#> d.Dp_d.R 3.891e-06
#> d.Dp_d.Dp 3.572e-05
Default algorithm
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)
Evaluation
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)
#>
#> Bayesian VAR model with p = 4
#>
#> Model:
#>
#> y ~ R.l01 + Dp.l01 + R.l02 + Dp.l02 + R.l03 + Dp.l03 + R.l04 + Dp.l04 + const + season.1 + season.2 + season.3
#>
#> Variable: R
#>
#> Mean SD Naive SD Time-series SD 2.5% 50%
#> R.l01 1.1618580 0.107931 0.0034131 0.0034131 0.942884 1.165e+00
#> Dp.l01 0.1931743 0.089310 0.0028242 0.0028242 0.019021 1.917e-01
#> R.l02 -0.2822331 0.161264 0.0050996 0.0054410 -0.600631 -2.804e-01
#> Dp.l02 -0.0227417 0.085735 0.0027112 0.0027112 -0.187809 -2.309e-02
#> R.l03 0.2431203 0.156611 0.0049525 0.0049525 -0.070666 2.480e-01
#> Dp.l03 0.1052096 0.089048 0.0028159 0.0028159 -0.055730 1.042e-01
#> R.l04 -0.2276779 0.106858 0.0033792 0.0040310 -0.439024 -2.288e-01
#> Dp.l04 0.1040068 0.083735 0.0026479 0.0026479 -0.057224 1.027e-01
#> const 0.0016440 0.004341 0.0001373 0.0001813 -0.006820 1.449e-03
#> season.1 0.0017077 0.005072 0.0001604 0.0001604 -0.007675 1.595e-03
#> season.2 0.0093040 0.005296 0.0001675 0.0001460 -0.001209 9.225e-03
#> season.3 -0.0001351 0.005040 0.0001594 0.0001522 -0.010260 -9.785e-05
#> 97.5%
#> R.l01 1.373889
#> Dp.l01 0.372050
#> R.l02 0.021192
#> Dp.l02 0.145084
#> R.l03 0.545996
#> Dp.l03 0.280562
#> R.l04 -0.015193
#> Dp.l04 0.279507
#> const 0.010739
#> season.1 0.012468
#> season.2 0.019146
#> season.3 0.009544
#>
#> Variable: Dp
#>
#> Mean SD Naive SD Time-series SD 2.5% 50%
#> R.l01 0.22116 0.100069 0.0031645 0.0031645 0.027522 0.22139
#> Dp.l01 0.04984 0.089053 0.0028161 0.0027241 -0.121803 0.04897
#> R.l02 -0.07106 0.154687 0.0048916 0.0046610 -0.369066 -0.07157
#> Dp.l02 -0.03476 0.089019 0.0028150 0.0030828 -0.197268 -0.03491
#> R.l03 0.01886 0.154435 0.0048837 0.0048837 -0.298118 0.02106
#> Dp.l03 0.05868 0.088388 0.0027951 0.0033079 -0.109995 0.05923
#> R.l04 -0.02168 0.102035 0.0032266 0.0035704 -0.217541 -0.02035
#> Dp.l04 0.36242 0.084178 0.0026619 0.0032900 0.197902 0.35994
#> const 0.01042 0.004094 0.0001295 0.0001748 0.002681 0.01035
#> season.1 -0.03410 0.004928 0.0001559 0.0001559 -0.044104 -0.03407
#> season.2 -0.01801 0.004995 0.0001580 0.0001580 -0.028383 -0.01811
#> season.3 -0.01639 0.004928 0.0001558 0.0001558 -0.026135 -0.01628
#> 97.5%
#> R.l01 0.419654
#> Dp.l01 0.228124
#> R.l02 0.230664
#> Dp.l02 0.134750
#> R.l03 0.317470
#> Dp.l03 0.229644
#> R.l04 0.182453
#> Dp.l04 0.533865
#> const 0.018695
#> season.1 -0.024122
#> season.2 -0.008187
#> season.3 -0.007100
#>
#> Variance-covariance matrix:
#>
#> Mean SD Naive SD Time-series SD 2.5% 50%
#> R_R 2.931e-05 4.400e-06 1.391e-07 1.391e-07 2.190e-05 2.888e-05
#> R_Dp -1.903e-06 3.047e-06 9.634e-08 8.959e-08 -8.123e-06 -1.940e-06
#> Dp_R -1.903e-06 3.047e-06 9.634e-08 8.959e-08 -8.123e-06 -1.940e-06
#> Dp_Dp 2.682e-05 4.350e-06 1.376e-07 1.471e-07 1.985e-05 2.654e-05
#> 97.5%
#> R_R 3.858e-05
#> R_Dp 3.695e-06
#> Dp_R 3.695e-06
#> Dp_Dp 3.676e-05
Forecasts
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)
Forecast error impulse response
FEIR <- irf(bvar_form, impulse = "R", response = "Dp", n.ahead = 20)
plot(FEIR, main = "Forecast Error Impulse Response", xlab = "Period", ylab = "Response")
Forecast error variance decomposition
bvar_fevd_oir <- fevd(bvar_form, response = "Dp", n.ahead = 20)
plot(bvar_fevd_oir, main = "OIR-based FEVD of inflation")
References
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
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.