# 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`:

    plot(bvec_est)

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