---
title: "Time Varying Parameters and Stochastic Volatility in Error Correction Models"
author: "Franz X. Mohr"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Time Varying Parameters and Stochastic Volatility in Error Correction Models}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
## 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.
``` r
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")
```
## Model set-up
``` r
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
``` r
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.
``` r
set.seed(20260912)
model <- add_initial_values(model)
```
``` r
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.
``` r
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 *
```
``` r
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.
``` r
plot(model)
```
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`.
``` r
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")
```
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.
``` r
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")
```
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.
``` r
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:
``` r
round(eigen(model_ml[["priors"]][["beta"]][["p_tau"]], symmetric = TRUE)$values, 3)
#> [1] 1.000 0.989 0.168
```
``` r
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:
``` r
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")
```
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:
``` r
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:
``` r
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:
``` r
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")
```
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.
``` r
model <- add_forecast_input(model, n_ahead = 8)
model <- add_posterior_forecasts(model)
plot(predict(model, n_ahead = 8), n_pre = 20)
```
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:
``` r
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](https://doi.org/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.
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.
Korobilis, D. (2013). VAR forecasting using Bayesian variable selection. *Journal of Applied Econometrics, 28*(2), 204-230.
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.
