Time Varying Parameters and Stochastic Volatility in Error Correction Models

Introduction

The companion vignette on TVP-SV-VAR models describes what the arguments tvp and error do and which priors they require. They mean the same thing for a vector error correction model, and create_bvecmodel takes them in the same way. What changes is the object they are applied to, and that change is not cosmetic: in a VEC model the time varying block includes the cointegration term, so the model allows the long-run relationship itself to move.

The estimated model is

\[\Delta y_t = \Pi_t w_t + \sum_{l=1}^{p-1} \Gamma_{l,t} \Delta y_{t-l} + C_t d_t + u_t, \qquad u_t \sim N(0, \Sigma_t),\]

with $\Pi_t = \alpha_t \beta_t^{\prime}$ of rank $r$. The coefficients of the non-cointegration part and the loadings $\alpha_t$ follow random walks and $\Sigma_t$ is the stochastic volatility specification of the VAR vignette. The cointegration vectors get a state equation of their own,

\[\beta_t = \rho \beta_{t-1} + \eta_t, \qquad \eta_t \sim N(0, I),\]

which is the specification of Koop, León-González and Strachan (2011). A constant coefficient VEC model instead puts a matric-variate prior on the cointegration space, so add_priors asks for different arguments in the two cases: coint$v_i and coint$p_tau_i for a constant model, coint$rho for a time varying one.

Data

Data set at_macrodata contains, among the domestic series of Austria, a short-term interest rate (r), a long-term interest rate (lr) and inflation (Dp). The model uses the three series in percent per quarter from 1979Q2 to 2019Q4 and asks the question this model is for, and which Koop, León-González and Strachan (2011) ask of UK data: whether the long-run relationship between the two rates and inflation is the same at the end of the sample as it was at its beginning.

library(bvartools)

data("at_macrodata")

data <- window(at_macrodata[["domestic"]][, c("r", "lr", "Dp")], end = c(2019, 4)) * 100

plot(data, main = "Austrian interest and inflation rates")

plot of chunk data

Model set-up

model <- create_bvecmodel(data, p = 2, r = 1, const = "unrestricted",
                          tvp = TRUE, error = "sv+covar", varsel = "bvs",
                          iterations = 3000, burnin = 1000)

model[["model"]][["algorithm"]]
#> [1] "VecTvpStochvol"

p is the lag order of the level VAR, so the model contains $p - 1 = 1$ lag of the differenced series. The cointegration rank is fixed at one, which for a term structure is the natural specification.

Models with a prior on the cointegration space are sensitive to the scale of the series in the error correction term, and scale_error_correction can divide them by the standard deviation of their differences. This model is estimated on the series as they are, for a reason explained with the priors below, so nothing has to be rescaled after estimation either.

Priors

model <- add_priors(model,
                    coef = list(v_i = 1, v_i_det = 1 / 10,
                                shape = 3, rate = 0.000001, rate_det = 0.01),
                    coint = list(rho = 0.999),
                    sigma = list(shape = 3, rate = 0.01,
                                 mu = 0, v_i = 0.01,
                                 state_variance = 0.05, offset = 1e-8),
                    varsel = list(inprior = 0.5, exclude_det = TRUE))

Arguments coef and sigma are the ones of the VAR vignette: prior precisions for the initial state of the coefficients, an inverse gamma prior for the variances of their state equation, and the prior of the stochastic volatility block.

coef$rate is two orders of magnitude smaller than in the VAR vignette, and the series in the error correction term are left unscaled. Both choices address the same problem. The coefficient paths can absorb most of the residuals of the short-term interest rate equation, whose residual standard deviation is then reported at a small fraction of the one of a least squares fit of the same regressors and almost flat over forty years – the symptom the VAR vignette describes. In four chains of the length used here, each with a different seed, that happened in every chain with rate = 0.0001, which put the residual standard deviation of that equation at between 4 and 53 percent of the least squares value, and in every chain with rate = 0.000001 when the series were scaled, at between 3 and 10 percent. With rate = 0.000001 and unscaled series, all four chains put it at about 95 percent.

Neither choice reaches all of the drift, though. The cointegration vectors take steps of unit variance whatever coef$rate is, and how far those steps move $\Pi_t$ depends on the prior precision of the loadings relative to the scale of $\Pi$ and on the levels of the series in the error correction term, as section ‘Prior on the cointegration space’ of ?add_priors.bvecmodel explains. Small rates alone therefore do not guarantee that a model passes the comparison with least squares, which is why that comparison is the check to make, whichever block would fit the residuals.

Argument coint is where a VEC model differs. rho is the autocorrelation of the state equation of $\beta$ and must be smaller than one. The prior on the state before the sample is that equation’s own stationary distribution, $N(0, I / (1 - \rho^2))$, which is what makes the prior proper, so a value close to one is intended and a value far from it draws a warning. Only the product $\alpha_t \beta_t^{\prime}$ is identified, and the prior variance of the loadings is shrunk by $1 - \rho^2$ to compensate, which leaves the product on the scale that coef$v_i asks for.

The loadings are never subject to variable selection, so the cointegration term is always in the model. Variable selection applies to the coefficients of the differenced regressors, and exclude_det = TRUE keeps the unrestricted constant out of it as well.

Initial values and posterior draws

Initial values of the state variances and volatilities are drawn from their priors, so the seed is set before them.

set.seed(20260912)
model <- add_initial_values(model)

model <- add_posterior_coefficients(model)

Evaluation

Summary statistics for a period

As for a TVP-SV-VAR, a summary table is the summary of a period, and the last one is used by default. The first block of columns is the cointegration matrix $\Pi_t$ rather than $\alpha_t$ and $\beta_t$ separately, since the factors are identified only up to an invertible $r \times r$ transformation while their product is not.

summary(model)
#> 
#> Bayesian TVP-SV-VEC model with p = 2 
#> 
#> Variable selection algorithm: Bayesian variable selection (Korobilis, 2013)
#> 
#> Endogenous variables: r, lr, Dp
#> 
#> Period: 161 
#> 
#> Variable: r 
#> 
#>                  Mean         SD     Naive SD Time-series SD        2.5%
#> l.r      -0.002498654 0.02045598 0.0003734734   0.0009380368 -0.05054164
#> l.lr     -0.002775249 0.02110438 0.0003853115   0.0006605850 -0.05339489
#> l.Dp      0.009239946 0.04651781 0.0008492951   0.0019721790 -0.08252893
#> d.r.l01   0.003389352 0.03545343 0.0006472881   0.0017520944 -0.05356195
#> d.lr.l01  0.202578621 0.07661494 0.0013987910   0.0036999522  0.00000000
#> d.Dp.l01 -0.002548845 0.01088882 0.0001988017   0.0012879152 -0.04381623
#> const     0.012968736 0.02375067 0.0004336258   0.0011102930 -0.04037986
#>                    50%      97.5% Incl. prob. 
#> l.r      -0.0004160284 0.03767029   1.0000000 
#> l.lr     -0.0006078286 0.04158929   1.0000000 
#> l.Dp      0.0077881922 0.11013957   1.0000000 
#> d.r.l01   0.0000000000 0.11304485   0.1020000 
#> d.lr.l01  0.2065857662 0.34351616   0.9693333 
#> d.Dp.l01  0.0000000000 0.00000000   0.0720000 
#> const     0.0149747008 0.05598934   1.0000000 
#> 
#> Variable: lr 
#> 
#>                   Mean          SD     Naive SD Time-series SD        2.5%
#> l.r      -7.314892e-04 0.025952633 4.738281e-04   5.721513e-04 -0.06071701
#> l.lr     -8.463268e-03 0.030936850 5.648270e-04   1.643052e-03 -0.08791141
#> l.Dp      1.722778e-02 0.062409379 1.139434e-03   1.648010e-03 -0.10883931
#> d.r.l01   2.286392e-03 0.017083117 3.118936e-04   6.509040e-04  0.00000000
#> d.lr.l01  2.543608e-02 0.068193781 1.245042e-03   4.566388e-03  0.00000000
#> d.Dp.l01 -9.506516e-05 0.002314801 4.226228e-05   7.264354e-05  0.00000000
#> const    -2.154395e-02 0.043378138 7.919728e-04   1.067755e-03 -0.10600798
#>                    50%      97.5% Incl. prob. 
#> l.r      -5.783027e-05 0.05376886   1.0000000 
#> l.lr     -2.379598e-03 0.04318360   1.0000000 
#> l.Dp      1.555078e-02 0.14475878   1.0000000 
#> d.r.l01   0.000000e+00 0.04804089   0.0780000 
#> d.lr.l01  0.000000e+00 0.24977720   0.1843333 
#> d.Dp.l01  0.000000e+00 0.00000000   0.0100000 
#> const    -2.086259e-02 0.06176590   1.0000000 
#> 
#> Variable: Dp 
#> 
#>                 Mean         SD     Naive SD Time-series SD        2.5%
#> l.r       0.01551050 0.35092524 0.0064069890    0.022891321 -0.69245441
#> l.lr      0.21359443 0.30835903 0.0056298398    0.027875870 -0.37945527
#> l.Dp     -0.82804085 0.23206010 0.0042368184    0.012514884 -1.29892873
#> d.r.l01   0.48389659 0.32607282 0.0059532480    0.016788728  0.00000000
#> d.lr.l01  0.01357624 0.17896592 0.0032674556    0.004513047 -0.39130970
#> d.Dp.l01  0.01254822 0.04552313 0.0008311348    0.002353052  0.00000000
#> const     0.33472474 0.14041238 0.0025635675    0.008738602  0.06544187
#>                  50%      97.5% Incl. prob.  
#> l.r       0.02792628  0.6652067   1.0000000  
#> l.lr      0.21122375  0.8318072   1.0000000  
#> l.Dp     -0.81862910 -0.4091824   1.0000000 *
#> d.r.l01   0.52046963  1.0621197   0.8116667  
#> d.lr.l01  0.00000000  0.4803338   0.2746667  
#> d.Dp.l01  0.00000000  0.1732222   0.1430000  
#> const     0.32998103  0.6064128   1.0000000 *
#> 
#> Variance-covariance matrix:
#> 
#>               Mean           SD     Naive SD Time-series SD         2.5%
#> r_r   0.0004006287 0.0006519833 1.190353e-05   4.643144e-05 3.908238e-07
#> r_lr  0.0000927324 0.0001510870 2.758459e-06   1.076273e-05 7.384997e-08
#> r_Dp  0.0004687951 0.0007442306 1.358773e-05   5.572635e-05 4.244548e-07
#> lr_lr 0.0026949518 0.0011332436 2.069010e-05   9.061756e-05 9.806436e-04
#> lr_Dp 0.0029027804 0.0017870810 3.262749e-05   9.990530e-05 2.311158e-04
#> Dp_Dp 0.0511828476 0.0184742867 3.372928e-04   1.009279e-03 2.414547e-02
#>                50%        97.5%  
#> r_r   1.884831e-04 0.0020620257 *
#> r_lr  4.002893e-05 0.0004874886 *
#> r_Dp  2.101437e-04 0.0024302265 *
#> lr_lr 2.525613e-03 0.0053977892 *
#> lr_Dp 2.604870e-03 0.0071282121 *
#> Dp_Dp 4.816964e-02 0.0961273539 *

period_of <- function(object, year, quarter) {
  which(stats::time(object[["data"]][["train"]][["y"]]) == year + (quarter - 1) / 4)
}

summary(model, period = period_of(model, 1985, 1))
#> 
#> Bayesian TVP-SV-VEC model with p = 2 
#> 
#> Variable selection algorithm: Bayesian variable selection (Korobilis, 2013)
#> 
#> Endogenous variables: r, lr, Dp
#> 
#> Period: 22 
#> 
#> Variable: r 
#> 
#>                  Mean         SD     Naive SD Time-series SD         2.5%
#> l.r      -0.024727489 0.03063389 0.0005592958    0.003479470 -0.101048536
#> l.lr     -0.030201496 0.03017514 0.0005509201    0.003084038 -0.104114477
#> l.Dp      0.115201420 0.06266083 0.0011440250    0.002585367 -0.001300008
#> d.r.l01   0.003408164 0.03553712 0.0006488161    0.001757510 -0.048960191
#> d.lr.l01  0.202567497 0.07671751 0.0014006637    0.003703166  0.000000000
#> d.Dp.l01 -0.002615366 0.01113142 0.0002032310    0.001388843 -0.044910971
#> const    -0.021736091 0.06560893 0.0011978497    0.005818083 -0.149932518
#>                  50%      97.5% Incl. prob. 
#> l.r      -0.01831887 0.02067567   1.0000000 
#> l.lr     -0.02464811 0.01464793   1.0000000 
#> l.Dp      0.11307566 0.24326960   1.0000000 
#> d.r.l01   0.00000000 0.11362155   0.1020000 
#> d.lr.l01  0.20650861 0.34236350   0.9693333 
#> d.Dp.l01  0.00000000 0.00000000   0.0720000 
#> const    -0.02189295 0.10821392   1.0000000 
#> 
#> Variable: lr 
#> 
#>                   Mean          SD     Naive SD Time-series SD        2.5%
#> l.r      -7.022402e-03 0.012573656 2.295625e-04   7.716274e-04 -0.03656391
#> l.lr     -1.036862e-02 0.014801274 2.702330e-04   1.294448e-03 -0.04649749
#> l.Dp      3.745882e-02 0.036910212 6.738852e-04   7.855834e-04 -0.03290696
#> d.r.l01   2.377112e-03 0.016936882 3.092237e-04   6.487778e-04  0.00000000
#> d.lr.l01  2.553780e-02 0.068357597 1.248033e-03   4.595080e-03  0.00000000
#> d.Dp.l01 -9.806079e-05 0.002361541 4.311564e-05   6.406874e-05  0.00000000
#> const    -9.201111e-03 0.036213469 6.611645e-04   1.367843e-03 -0.08113406
#>                   50%      97.5% Incl. prob. 
#> l.r      -0.004546849 0.01435386   1.0000000 
#> l.lr     -0.006911344 0.01117245   1.0000000 
#> l.Dp      0.037284501 0.10924475   1.0000000 
#> d.r.l01   0.000000000 0.05234066   0.0780000 
#> d.lr.l01  0.000000000 0.25437243   0.1843333 
#> d.Dp.l01  0.000000000 0.00000000   0.0100000 
#> const    -0.009138791 0.06249934   1.0000000 
#> 
#> Variable: Dp 
#> 
#>                 Mean         SD     Naive SD Time-series SD        2.5%
#> l.r       0.19829265 0.19696738 0.0035961158    0.023079824 -0.19371389
#> l.lr      0.27701483 0.21458603 0.0039177870    0.037041415 -0.09939393
#> l.Dp     -1.02426330 0.17632114 0.0032191688    0.008967970 -1.40841504
#> d.r.l01   0.48364279 0.32604041 0.0059526562    0.016782942  0.00000000
#> d.lr.l01  0.01357996 0.17894439 0.0032670627    0.004513039 -0.39392692
#> d.Dp.l01  0.01261751 0.04563327 0.0008331458    0.002300732  0.00000000
#> const     0.01258847 0.29199929 0.0053311533    0.049141454 -0.57819240
#>                  50%      97.5% Incl. prob.  
#> l.r       0.20618031  0.5638109   1.0000000  
#> l.lr      0.26190785  0.7212998   1.0000000  
#> l.Dp     -1.01661609 -0.6978181   1.0000000 *
#> d.r.l01   0.52073128  1.0673605   0.8116667  
#> d.lr.l01  0.00000000  0.4815487   0.2746667  
#> d.Dp.l01  0.00000000  0.1736515   0.1430000  
#> const     0.03404646  0.5515325   1.0000000  
#> 
#> Variance-covariance matrix:
#> 
#>              Mean          SD     Naive SD Time-series SD        2.5%
#> r_r   0.027415701 0.017852707 3.259444e-04   0.0006748980 0.007492173
#> r_lr  0.006194967 0.004368520 7.975789e-05   0.0001599025 0.001376332
#> r_Dp  0.032525923 0.024160024 4.410997e-04   0.0011223228 0.007088188
#> lr_lr 0.003870373 0.001551034 2.831788e-05   0.0001294902 0.001787518
#> lr_Dp 0.009884457 0.006179528 1.128222e-04   0.0002391376 0.002994233
#> Dp_Dp 0.164531632 0.053639847 9.793251e-04   0.0042006867 0.087772035
#>               50%       97.5%  
#> r_r   0.023089770 0.073534158 *
#> r_lr  0.005170804 0.017456977 *
#> r_Dp  0.026819258 0.092845246 *
#> lr_lr 0.003550990 0.007686343 *
#> lr_Dp 0.008373436 0.025435736 *
#> Dp_Dp 0.156124730 0.288342413 *

The inclusion probabilities of the $\Pi$ columns are one by construction, as explained above, and the inclusion probabilities of the remaining columns do not carry a period: variable selection decides on a regressor for the whole sample.

Coefficient and volatility paths

plot draws one figure per block of coefficients: the cointegration matrix $\Pi_t$, the lagged differences, the deterministic terms and the covariance matrix of the error term, each on its own.

plot(model)

plot of chunk paths

In the figure of $\Pi_t$, reading down a column gives the response of each differenced series to the level of one variable, and reading across a row gives the error correction of one equation. A panel that sits at zero among the $\Gamma$ coefficients is one that variable selection has switched off.

A single element of $\Pi_t$ is easier to read on its own. The following computes it from the draws of the loadings, which are the first $Kr$ coefficients of each period of posterior$a$coeffs, and the draws of the cointegration vectors in posterior$beta$coeffs.

pi_path <- function(object, equation, variable, ci = 0.9) {
  k <- object[["model"]][["k"]]
  r <- object[["model"]][["rank"]]
  m <- ncol(object[["data"]][["train"]][["z"]])
  k_ect <- ncol(object[["data"]][["train"]][["w"]])
  tt <- nrow(object[["data"]][["train"]][["y"]])
  draws <- nrow(object[["posterior"]][["a"]][["coeffs"]])

  i <- which(object[["model"]][["endogen"]] == equation)
  j <- which(dimnames(object[["data"]][["train"]][["w"]])[[2]] == variable)

  result <- matrix(NA_real_, draws, tt)
  for (period in 1:tt) {
    alpha <- object[["posterior"]][["a"]][["coeffs"]][, (period - 1) * m + 1:(k * r)]
    beta <- object[["posterior"]][["beta"]][["coeffs"]][, (period - 1) * k_ect * r + 1:(k_ect * r)]
    result[, period] <- rowSums(matrix(alpha, draws, k * r)[, k * 0:(r - 1) + i, drop = FALSE] *
                                  matrix(beta, draws, k_ect * r)[, k_ect * 0:(r - 1) + j, drop = FALSE])
  }

  bands <- t(apply(result, 2, stats::quantile,
                   probs = c((1 - ci) / 2, .5, 1 - (1 - ci) / 2)))

  stats::ts(bands, start = stats::start(object[["data"]][["train"]][["y"]]),
            frequency = stats::frequency(object[["data"]][["train"]][["y"]]))
}

path <- pi_path(model, equation = "r", variable = "l.r")

stats::ts.plot(path, lty = c(2, 1, 2),
               main = "Short-term rate: adjustment to its own level",
               ylab = "Element of Pi")
abline(h = 0, lty = "dotted")

plot of chunk pi-path

The band of this element is widest in the early 1980s and narrows over the sample, while its median rises from about -0.04 to zero: whatever adjustment of the short-term rate to its own level there is fades out over the four decades.

Note what is not plotted here. The path of $\beta_t$ on its own is not a quantity to read: a cointegration vector is identified only up to normalisation, and normalising it on a variable whose coefficient wanders near zero produces a path that swings between large values without anything having happened to the model. $\Pi_t$ carries the same information and is identified.

The volatilities are extracted from the draws of the error precision exactly as in the VAR case.

volatility <- function(object) {
  k <- object[["model"]][["k"]]
  kk <- k * k
  tt <- nrow(object[["data"]][["train"]][["y"]])
  draws <- object[["posterior"]][["u_sigma_inv"]][["coeffs"]]

  result <- matrix(NA_real_, tt, k)
  for (i in 1:tt) {
    result[i, ] <- rowMeans(apply(draws[, (i - 1) * kk + 1:kk], 1,
                                  function(x) {sqrt(diag(solve(matrix(x, k))))}))
  }

  stats::ts(result, start = stats::start(object[["data"]][["train"]][["y"]]),
            frequency = stats::frequency(object[["data"]][["train"]][["y"]]),
            names = object[["model"]][["endogen"]])
}

plot(volatility(model),
     main = "Posterior mean of the residual standard deviations")

plot of chunk volatility

For these data the volatility block does most of its work in the equation of the short-term interest rate. Its residual standard deviation falls from about a quarter of a percentage point in the early 1980s to a tenth of that once the rate approached its lower bound, with a brief return during the financial crisis. A model that averaged over the sample would report credible bands for the 2010s that are far too wide. The volatilities of the other two equations move much less: the one of inflation falls by about half until 2000, the one of the long-term rate by about a quarter, most of it in the first half of the 1980s.

A prior centred on the maximum likelihood estimate

The prior used above says nothing about where the cointegration space lies: under $\beta_t = \rho \beta_{t-1} + \eta_t$ the marginal prior of $sp(\beta_t)$ is uniform in every period. The working paper version of Koop, León-González and Strachan (2011) shows how to centre it on a given space $sp(H)$ instead, by putting $P_\tau = H H^{\prime} + H_\perp T H_\perp^{\prime}$ into the transition of the state equation,

\[\beta_t = \rho (I_r \otimes P_\tau) \beta_{t-1} + \eta_t, \qquad \eta_t \sim N(0, I).\]

The part of $\beta_t$ along $sp(H)$ keeps the autocorrelation $\rho$ and the part off it decays faster, so the mode of the marginal distribution of $sp(\beta_t)$ is $sp(H)$ in every period, while the space of one period is still centred near the space of the period before.

With coint = list(rho = 0.999, p_tau_i = "ml", weight = 1), add_priors takes $H$ from Johansen’s (1995) maximum likelihood estimate of the space and chooses $T$ so that the prior spread of $sp(\beta_t)$ around it in each period matches the sampling variance of that estimate, with weight counting how many samples’ worth of information the prior adds. The state before the sample gets the stationary distribution of the new state equation. The estimate is computed from the error correction term as it is stored in the model, here on the scale of the data.

model_ml <- create_bvecmodel(data, p = 2, r = 1, const = "unrestricted",
                             tvp = TRUE, error = "sv+covar", varsel = "bvs",
                             iterations = 3000, burnin = 1000)

model_ml <- add_priors(model_ml,
                       coef = list(v_i = 1, v_i_det = 1 / 10,
                                   shape = 3, rate = 0.000001, rate_det = 0.01),
                       coint = list(rho = 0.999, p_tau_i = "ml", weight = 1),
                       sigma = list(shape = 3, rate = 0.01,
                                    mu = 0, v_i = 0.01,
                                    state_variance = 0.05, offset = 1e-8),
                       varsel = list(inprior = 0.5, exclude_det = TRUE))

How informative this prior can be is bounded by $\rho$: even $T = 0$ leaves the space $1 - \rho^2$ of room to move in each period. Where the sampling variance of the estimate in some direction is smaller than that, add_priors sets $T$ to zero in that direction and warns. On these data both directions off $sp(H)$ stay above the bound, so the prior matches the sampling variance without a warning. The eigenvalues of the transition that was stored show the outcome – one for the direction of $sp(H)$ itself, and one for each direction off it:

round(eigen(model_ml[["priors"]][["beta"]][["p_tau"]], symmetric = TRUE)$values, 3)
#> [1] 1.000 0.989 0.168

set.seed(20260913)
model_ml <- add_initial_values(model_ml)
model_ml <- add_posterior_coefficients(model_ml)

The same element of $\Pi_t$ as above, with the bands under the noninformative prior in grey:

path_ml <- pi_path(model_ml, equation = "r", variable = "l.r")

stats::ts.plot(path, path_ml, lty = c(2, 1, 2, 2, 1, 2),
               col = rep(c("grey60", "black"), each = 3),
               main = "Short-term rate: adjustment to its own level",
               ylab = "Element of Pi")
abline(h = 0, lty = "dotted")
legend("bottomleft", legend = c("Noninformative", "Centred on ML"),
       col = c("grey60", "black"), lty = 1, bty = "n")

plot of chunk ml-prior-path

Averaged over the sample, the width of the 90% bands of every element of $\Pi_t$ under the prior centred on the maximum likelihood estimate, relative to the noninformative one – rows are the variables in the error correction term, columns the equations:

band_width <- function(object) {
  sapply(object[["model"]][["endogen"]], function(equation) {
    sapply(dimnames(object[["data"]][["train"]][["w"]])[[2]], function(variable) {
      bands <- pi_path(object, equation = equation, variable = variable)
      mean(bands[, 3] - bands[, 1])
    })
  })
}

round(band_width(model_ml) / band_width(model), 2)
#>         r   lr   Dp
#> l.r  0.72 0.71 0.85
#> l.lr 0.81 0.70 0.97
#> l.Dp 0.81 0.82 1.07

A ratio below one means that the prior centred on the maximum likelihood estimate narrows the band of that element. Here it does so for all elements but one, most for the adjustment to the lagged long-term rate, while the band of the adjustment of inflation to its own level is about five percent wider than under the noninformative prior. It need not narrow every band, and that is a property of how the prior is built rather than a failure of it: a larger weight makes the prior tighter by lowering $T$, but $T$ is also how much of the position of the space carries over from one period to the next. In a direction where $T$ is zero, the position does not carry over at all, each period is informed by its own data alone, and the posterior in that direction can be wider than under the noninformative prior – see ?add_priors.bvecmodel. How this plays out for the elements of $\Pi_t$ depends on how the directions of the space line up with the variables.

Two further things to keep in mind. First, the centre is estimated from the same sample the model is estimated on, so the data enter twice, and where the bands do narrow they overstate how precisely the space is known. Second, and most important for these data, $H$ is the estimate of a constant space. Centring every period on it pulls the path towards the average of the sample, which works against the question the model is estimated to answer, whether the long-run relationship moved. A prior like this is best suited to a model whose cointegration space is expected to move around a stable centre, rather than to drift away from one.

Impulse responses and forecasts

Impulse responses and variance decompositions of a VEC model are obtained from its VAR representation in levels, which vec_to_var computes draw by draw. A time varying model is a VEC model per period, so the transformation is applied period by period, and the result is a VAR model with time varying parameters:

bvar_form <- vec_to_var(model)

bvar_form[["model"]][["algorithm"]]
#> [1] "VarTvpStochvol"

The draws of the variances of the state equations of the VEC coefficients and of $\rho$ are not carried over: they describe how the coefficients of the VEC model drift and have no counterpart among the coefficients in levels. The paths of the coefficients and the draws of the error term are, and they are all the application functions need.

Impulse responses and forecast error variance decompositions are those of one period, which argument period selects, the last period being the default. The response of inflation to a shock to the short-term interest rate in 1985Q1:

feir <- irf(bvar_form, impulse = "r", response = "Dp", n_ahead = 20,
            period = period_of(model, 1985, 1))

plot(feir, main = "Forecast error impulse response, 1985Q1",
     xlab = "Period", ylab = "Response")

plot of chunk feir-1985

Forecasts are made from the VEC model itself. add_forecast_input prepares the regressors of the forecast periods, in levels, and add_posterior_forecasts simulates each draw forward from the last period of the sample: the loadings, the short-run coefficients and the volatilities take a step of their random walks every quarter, the cointegration vectors a step of their state equation, and the VAR in levels is rebuilt from them in every period. The bands therefore reflect both the uncertainty about where the coefficients are at the end of the sample and the drift the model allows after it.

model <- add_forecast_input(model, n_ahead = 8)
model <- add_posterior_forecasts(model)

plot(predict(model, n_ahead = 8), n_pre = 20)

plot of chunk forecasts

The VAR representation above has no state equation for its coefficients in levels, so it can only hold them at their values in the last period. forecast_states = "hold" gives that forecast from the VEC model as well. The width of the 90% intervals eight quarters ahead shows what the drift adds:

held <- add_posterior_forecasts(model, forecast_states = "hold")

interval_width <- function(object) {
  fcst <- predict(object, n_ahead = 8)[["fcst"]]
  apply(fcst[8, , ], 1, function(x) diff(stats::quantile(x, c(0.05, 0.95))))
}

round(rbind(simulated = interval_width(model), held = interval_width(held)), 2)
#>              r   lr   Dp
#> simulated 1.83 2.44 1.67
#> held      0.60 1.19 1.10

Eight quarters ahead, carrying the drift forward about doubles the interval of the long-term rate, widens that of inflation by half and makes that of the short-term rate between two and three times as wide. A forecast that holds the coefficients at the end of the sample understates the uncertainty this model implies.

The constant coefficient model with stochastic volatility, estimated with tvp = FALSE, remains the model a TVP-SV specification has to beat. The two can be compared out of sample by estimating both over an expanding window with use_expanding_window and forecasting the models of its windows directly: add_forecast_input and add_posterior_forecasts take VEC models, forecast them in levels and, for the time varying one, simulate its cointegration vectors and volatility forward over the forecast horizon.

Citing bvartools

If you use bvartools in published work, please cite it. citation("bvartools") prints the reference, and the package has the DOI 10.5281/zenodo.22736604, which always resolves to the latest archived version.

References

Johansen, S. (1995). Likelihood-based inference in cointegrated vector autoregressive models. Oxford: Oxford University Press.

Koop, G., León-González, R., & Strachan R. W. (2011). Bayesian inference in a time varying cointegration model. Journal of Econometrics, 165(2), 210-220. https://doi.org/10.1016/j.jeconom.2011.07.007

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

Korobilis, D. (2013). VAR forecasting using Bayesian variable selection. Journal of Applied Econometrics, 28(2), 204-230. https://doi.org/10.1002/jae.1271

Lütkepohl, H. (2006). New introduction to multiple time series analysis (2nd ed.). Berlin: Springer.

Primiceri, G. E. (2005). Time varying structural vector autoregressions and monetary policy. Review of Economic Studies, 72(3), 821-852. https://doi.org/10.1111/j.1467-937X.2005.00353.x