Bayesian BEKK-GARCH (NumPyro)¶

Multivariate volatility, part 3c: a full-posterior diagonal BEKK¶

Part 3 fit the diagonal BEKK by maximum likelihood. Here we give it the Bayesian treatment — a full posterior over the covariance dynamics — completing the symmetry with the Bayesian DCC (dcc_numpyro.ipynb). We keep to the diagonal BEKK deliberately: the full BEKK is only identified up to the sign of $A$ and orthogonal rotations (only $A'A$ and $B'B$ are identified), which splits MCMC chains across equivalent modes. The diagonal form has none of that ambiguity and samples cleanly.

The model (bivariate), with $A=\text{diag}(a_1,a_2)$, $B=\text{diag}(b_1,b_2)$ and $C$ lower-triangular: $$H_t = CC' + A'\,\varepsilon_{t-1}\varepsilon_{t-1}'\,A + B'\,H_{t-1}B,\qquad H_t[i,j]=(CC')[i,j]+a_ia_j\,\varepsilon_{i,t-1}\varepsilon_{j,t-1}+b_ib_j\,H_{t-1}[i,j].$$ Positive-definiteness of $H_t$ is automatic (BEKK's signature); the news impact per series is $a_i^2$ and the persistence $b_i^2$, with stationarity $a_i^2+b_i^2<1$.

Data: DJIA & NASDAQ daily returns, equity_returns.csv (1993–2024) — the same series as bekk_python.ipynb, so the posterior can be checked against the MLE ($a_i^2\approx0.076/0.067$, $b_i^2\approx0.913/0.925$).

Model specification and priors¶

Loadings of variance. The awkward part of a Bayesian BEKK is the stationarity constraint $a_i^2+b_i^2<1$ together with the near-unit-root geometry ($a_i^2+b_i^2\approx0.99$ in equity data — a boundary funnel). We handle both at once with a simplex prior: for each series draw $$\big(a_i^2,\ b_i^2,\ \text{rest}_i\big)\sim\text{Dirichlet}(1,8,1),\qquad a_i=\sqrt{a_i^2}\ge0,\quad b_i=\sqrt{b_i^2}\ge0,$$ so $a_i^2+b_i^2<1$ holds by construction (the "rest" slack is what keeps it strictly inside), the $\text{Dirichlet}(1,8,1)$ mean puts $\approx0.8$ mass on persistence, and $a_i,b_i$ are taken positive (identifying their sign — only the products $a_ia_j,\,b_ib_j$ enter, so a common sign is free and we fix it positive).

Intercept. $C$ is lower-triangular with a positive diagonal, $c_{11},c_{22}\sim\text{HalfNormal}(0.5)$, $c_{21}\sim N(0,0.5)$; $CC'$ is then p.d. and small (it targets the unconditional level).

Likelihood. For the bivariate case everything is a $2\times2$ closed form — no Cholesky is needed: $$\det H_t = h_{11,t}h_{22,t}-h_{12,t}^2,\qquad r_t'H_t^{-1}r_t=\frac{h_{22,t}r_{1t}^2-2h_{12,t}r_{1t}r_{2t}+h_{11,t}r_{2t}^2}{\det H_t},$$ $$\log p(r_t\mid\cdot)=-\tfrac12\big[2\log2\pi+\log\det H_t + r_t'H_t^{-1}r_t\big].$$ (Using the closed form rather than MultivariateNormal also keeps the notebook robust to environments where JAX's LAPACK Cholesky misbehaves.)

Estimation algorithm¶

The conditional covariance $H_t$ is a deterministic recursion in the parameters and past returns, written as a jax.lax.scan that carries the three unique entries $(h_{11},h_{12},h_{22})$ forward — the exact BEKK filter of bekk.py, differentiable for HMC. There are no latent states (unlike factor SV): only the 7 static parameters $(c_{11},c_{21},c_{22},a_1,a_2,b_1,b_2)$, so this is a small, fast target. We sample with NUTS (JAX/NumPyro — PyMC's scan is unusable on the build that generated the sibling notebooks), 2 chains, init_to_median, and check $\hat R$. The Dirichlet parameterisation keeps the sampler off the $a_i^2+b_i^2=1$ boundary, so no special funnel handling is needed beyond that.

In [3]:
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import jax, jax.numpy as jnp, numpyro, numpyro.distributions as dist
from numpyro.infer import MCMC, NUTS, init_to_median
import numpyro.diagnostics as diag
numpyro.enable_x64(); numpyro.set_host_device_count(2)

df = pd.read_csv("equity_returns.csv", parse_dates=["date"], index_col="date")
dates = df.index; rf = df.values - df.values.mean(0); r = jnp.asarray(rf); T, N = r.shape
Sc = np.cov(rf, rowvar=False); H0 = jnp.array([Sc[0, 0], Sc[0, 1], Sc[1, 1]])   # initialise H_1 at sample cov

def bekk_paths(cdiag, c21, a, b):                          # deterministic H_t = (h11,h12,h22), (T,3)
    CC00 = cdiag[0] ** 2; CC10 = c21 * cdiag[0]; CC11 = c21 ** 2 + cdiag[1] ** 2
    aa = jnp.outer(a, a); bb = jnp.outer(b, b)
    def step(Hp, ep):
        h11 = CC00 + aa[0, 0] * ep[0] * ep[0] + bb[0, 0] * Hp[0]
        h12 = CC10 + aa[0, 1] * ep[0] * ep[1] + bb[0, 1] * Hp[1]
        h22 = CC11 + aa[1, 1] * ep[1] * ep[1] + bb[1, 1] * Hp[2]
        Hn = jnp.array([h11, h12, h22]); return Hn, Hn
    _, Hr = jax.lax.scan(step, H0, r[:-1])
    return jnp.concatenate([H0[None], Hr], 0)

def model(r):
    cdiag = numpyro.sample("cdiag", dist.HalfNormal(0.5).expand([2]).to_event(1))
    c21 = numpyro.sample("c21", dist.Normal(0., 0.5))
    w = numpyro.sample("w", dist.Dirichlet(jnp.array([1., 8., 1.])).expand([2]).to_event(1))  # [a^2, b^2, rest]
    a = jnp.sqrt(w[:, 0]); b = jnp.sqrt(w[:, 1])
    numpyro.deterministic("news", w[:, 0]); numpyro.deterministic("persist", w[:, 1])
    H = bekk_paths(cdiag, c21, a, b)
    det = H[:, 0] * H[:, 2] - H[:, 1] ** 2
    quad = (H[:, 2] * r[:, 0] ** 2 - 2 * H[:, 1] * r[:, 0] * r[:, 1] + H[:, 0] * r[:, 1] ** 2) / det
    numpyro.factor("ll", -0.5 * jnp.sum(2 * jnp.log(2 * jnp.pi) + jnp.log(det) + quad))

mcmc = MCMC(NUTS(model, target_accept_prob=0.9, init_strategy=init_to_median),
            num_warmup=1000, num_samples=500, num_chains=2, chain_method="sequential", progress_bar=False)
mcmc.run(jax.random.PRNGKey(0), r=r)
S = mcmc.get_samples(); news = np.array(S["news"]); persist = np.array(S["persist"])
sg = diag.summary(mcmc.get_samples(group_by_chain=True), prob=0.9)
print("Bayesian diagonal BEKK (NumPyro NUTS), DJIA-NASDAQ:")
print("  news a^2:    mean %s   (MLE 0.076 / 0.067)" % np.round(news.mean(0), 3))
print("  persist b^2: mean %s   (MLE 0.913 / 0.925)" % np.round(persist.mean(0), 3))
print("  a^2+b^2:     mean %s" % np.round((news + persist).mean(0), 3))
print("  r_hat w max %.3f" % float(np.nanmax(sg["w"]["r_hat"])))

Bayesian diagonal BEKK (NumPyro NUTS), DJIA-NASDAQ:
  news a^2:    mean [0.077 0.068]   (MLE 0.076 / 0.067)
  persist b^2: mean [0.913 0.924]   (MLE 0.913 / 0.925)
  a^2+b^2:     mean [0.99  0.993]
  r_hat w max 1.007
In [4]:
# posterior conditional-correlation band: recompute the H_t path for a batch of draws in one vmap (closed form, no Cholesky)
idx = np.linspace(0, len(news) - 1, 200).astype(int)
Hp = jax.vmap(bekk_paths, in_axes=(0, 0, 0, 0))(
    jnp.asarray(S["cdiag"][idx]), jnp.asarray(S["c21"][idx]),
    jnp.sqrt(jnp.asarray(S["w"][idx][:, :, 0])), jnp.sqrt(jnp.asarray(S["w"][idx][:, :, 1])))
H = np.array(Hp)                                           # (ndraw, T, 3)
rho = H[:, :, 1] / np.sqrt(H[:, :, 0] * H[:, :, 2])
rlo, rmd, rhi = np.percentile(rho, [5, 50, 95], axis=0)
mle_rho = np.load("bekk_mle_corr.npy")

fig, ax = plt.subplots(1, 2, figsize=(15, 4.6), gridspec_kw=dict(width_ratios=[2, 1]))
ax[0].fill_between(dates, rlo, rhi, color="darkgreen", alpha=.22, lw=0, label="90% posterior band")
ax[0].plot(dates, rmd, color="darkgreen", lw=.6, label="posterior median")
ax[0].plot(dates, mle_rho, color="black", lw=.4, alpha=.6, label="MLE (bekk_python)")
for dt in ("2000-03-10", "2008-09-15", "2020-03-16"): ax[0].axvline(pd.Timestamp(dt), color="0.6", lw=.7, ls=":")
ax[0].set_ylim(0.2, 1.0); ax[0].set_ylabel("conditional correlation"); ax[0].legend(fontsize=8, loc="lower right")
ax[0].set_title("Bayesian BEKK: posterior of the DJIA-NASDAQ correlation path")
for l in ax[0].get_xticklabels(): l.set_rotation(20)
ax[1].scatter(news[:, 0] + persist[:, 0], news[:, 1] + persist[:, 1], s=3, alpha=.15, color="darkgreen")
ax[1].set_xlabel("$a_1^2+b_1^2$ (DJIA)"); ax[1].set_ylabel("$a_2^2+b_2^2$ (NASDAQ)")
ax[1].set_title("Posterior persistence per series\n(near the stationarity boundary)")
plt.tight_layout(); plt.show()
No description has been provided for this image

Results¶

  • The posterior concentrates on the MLE. The news $a_i^2$ and persistence $b_i^2$ posteriors sit around the maximum-likelihood values ($a_i^2\approx0.07$, $b_i^2\approx0.92$), now with full uncertainty and the stationarity constraint respected exactly. Both series' $a_i^2+b_i^2$ posteriors crowd just below 1 — the near-integrated persistence, quantified rather than point-estimated.
  • A posterior over the whole covariance path. The conditional-correlation band tracks the MLE path (low in the late-1990s tech divergence, high in every crisis) and, like the Bayesian DCC, hands a decision-maker a distribution over $H_t$ rather than a single trajectory.
  • Why diagonal, and where this sits. Keeping $A,B$ diagonal sidesteps the full-BEKK identification problem, so NUTS mixes without the multi-modality that plagues the full form. This completes the Bayesian trio of the covariance-parameterising models — DCC, and now BEKK — alongside the natively-Bayesian factor SV.

References¶

  • Engle, R. F. & Kroner, K. F. (1995). Multivariate simultaneous generalized ARCH. Econometric Theory 11, 122–150.
  • Vrontos, I. D., Dellaportas, P. & Politis, D. N. (2003). A full-factor multivariate GARCH model. Econometrics Journal 6, 312–334.