Factor stochastic volatility (from scratch, NumPyro)¶

Multivariate volatility, part 4: Aguilar & West (2000) — the Bayesian capstone¶

MGARCH (CCC/DCC/BEKK) treats the conditional covariance as a deterministic function of past returns. Multivariate stochastic volatility (MSV) instead makes the volatilities latent random processes — their own hidden dynamics — which is why MSV lives natively in the Bayesian world: with latent states, MCMC is not an option but the estimator. The scalable form is the factor SV model of Aguilar & West (2000): a few latent factors, each with its own stochastic volatility, drive the whole cross-section.

The model. With $N$ assets and one common factor, $$r_t = B\,f_t + e_t,\qquad f_t\sim N\!\big(0,e^{h^f_t}\big),\qquad e_{it}\sim N\!\big(0,e^{h^e_{it}}\big),$$ and every log-variance follows its own AR(1) stochastic-volatility process, $$h_t = \mu + \phi\,(h_{t-1}-\mu) + \sigma\,\eta_t.$$ Marginalising the factor gives a covariance that is low-rank-plus-diagonal, hence positive-definite for free and cheap to invert: $$H_t = e^{h^f_t}\,BB' + \operatorname{diag}\!\big(e^{h^e_t}\big).$$ This is the best of both earlier worlds — BEKK's guaranteed positive-definiteness and DCC's scalability — because the factor structure needs only $N$ loadings and $N{+}1$ SV processes, not $O(N^2)$ parameters.

Data. Daily returns on five SPDR sector ETFs — the S&P 500 sub-indices tradeable as funds:

ticker sector
XLF Financials
XLK Technology
XLE Energy
XLV Health Care
XLU Utilities

Downloaded from Yahoo Finance and stored in etf_returns.csv as daily log returns in percent, $r_{it}=100\times\ln(P_{it}/P_{i,t-1})$, over 5,283 common trading days, Jan 2004 – Dec 2024. The five are deliberately chosen to span the spectrum of market sensitivity — cyclical, market-driven sectors (Financials, Technology, Energy) alongside defensive, more idiosyncratic ones (Health Care, Utilities) — so the common market factor and the sector-specific idiosyncrasies are both clearly identified.

Model specification¶

Write $M = N+1$ for the number of latent volatility processes (one factor $+$ $N$ idiosyncratic). The full generative model is hierarchical:

Observation layer. Conditional on the latent log-variances and loadings, $$r_t = B\,f_t + e_t,\qquad f_t \sim N\!\big(0,\,e^{h^f_t}\big),\qquad e_{it}\sim N\!\big(0,\,e^{h^e_{it}}\big),\quad i=1,\dots,N,$$ with $f_t$ (scalar, one factor) and $e_t$ mutually and serially independent given the volatilities. $B=(b_1,\dots,b_N)'$ is the $N\times1$ vector of factor loadings.

Volatility layer. Each of the $M$ log-variances is an independent stationary AR(1) stochastic-volatility process, $$h_{j,t} = \mu_j + \phi_j\,(h_{j,t-1}-\mu_j) + \sigma_j\,\eta_{j,t},\qquad \eta_{j,t}\sim N(0,1),\qquad h_{j,1}\sim N\!\Big(\mu_j,\ \tfrac{\sigma_j^2}{1-\phi_j^2}\Big),$$ where $\mu_j$ is the mean log-variance (the vol level), $\phi_j\in(0,1)$ the persistence, and $\sigma_j$ the volatility-of-volatility; the initial draw is from the AR(1) stationary distribution.

Implied covariance. Because $f_t$ is Gaussian it can be integrated out analytically, giving $r_t\mid h_t \sim N(0,H_t)$ with the diagonal-plus-rank-one conditional covariance $$H_t = e^{h^f_t}\,BB' + \operatorname{diag}\!\big(e^{h^e_{1t}},\dots,e^{h^e_{Nt}}\big),$$ positive-definite by construction and requiring only $N$ loadings $+$ $M$ SV triples $(\mu_j,\phi_j,\sigma_j)$ — $O(N)$ parameters, versus $O(N^2)$ for BEKK.

Priors and identification. A one-factor model has a scale/sign indeterminacy ($f_t\!\to\!cf_t,\,B\!\to\!B/c$); we fix it by setting $b_1\equiv1$ and constraining the remaining loadings positive (a market factor loads positively on every sector — this also removes the sign multimodality that makes factor models hard to sample): $$b_1=1,\quad b_i\sim\text{HalfNormal}(1.5)\ (i\ge2),\qquad \mu_j\sim N(0,2),\quad \phi_j\sim\text{Beta}(20,1.5),\quad \sigma_j\sim\text{HalfNormal}(0.3).$$ The $\text{Beta}(20,1.5)$ prior has mean $\approx0.93$, encoding the strong persistence typical of volatility; $b_i$ and $\sigma_j$ are half-normal to enforce positivity.

Estimation algorithm¶

The model carries $M\times T \approx 32{,}000$ latent log-volatilities, so inference is a genuine stress-test for MCMC. Four ingredients make it both correct and well-mixed:

1. Marginalise the factors. The $f_t$ are integrated out in closed form (above), so the sampler works only with the log-vol states and the static parameters — never the factors themselves.

2. Woodbury likelihood — no $N\times N$ Cholesky. Evaluating $N(0,H_t)$ naively needs a Cholesky of each $H_t$; JAX's LAPACK Cholesky crashes on this Windows build (the same BLAS-DLL failure that felled rmgarch). But $H_t=D_t+e^{h^f_t}BB'$ with $D_t=\operatorname{diag}(e^{h^e_t})$ is rank-one-plus-diagonal, so the Sherman–Morrison / matrix-determinant lemma reduces everything to scalars. Writing $D_t^{-1}=\operatorname{diag}(e^{-h^e_t})$ and $$m_t = e^{-h^f_t} + \textstyle\sum_i b_i^2\,e^{-h^e_{it}},\qquad c_t = \textstyle\sum_i b_i\,r_{it}\,e^{-h^e_{it}},$$ the quadratic form and log-determinant needed for the likelihood are $$r_t'H_t^{-1}r_t = \sum_i r_{it}^2\,e^{-h^e_{it}} - \frac{c_t^2}{m_t},\qquad \log|H_t| = \sum_i h^e_{it} + h^f_t + \log m_t.$$ These are exactly the lines in the code cell below — pure elementwise arithmetic, crash-proof and much faster than a dense solve.

3. Non-centred parameterisation. Sampling the log-vol states directly creates the classic SV funnel (state variance depends on $\sigma_j$). Instead we sample standardised innovations $\varepsilon_{j,t}\sim N(0,1)$ and construct $h_{j,t}$ from them via the AR(1) recursion (ar1 below), which decouples the states from $\sigma_j$ and lets NUTS explore efficiently.

4. NUTS with careful tuning. We run the No-U-Turn sampler (JAX/NumPyro — PyMC's scan is unusable here, as in the Bayesian-DCC notebook) with init_to_median, target_accept=0.95 (short steps for the curved geometry), 1000 warmup $+$ 500 draws, 2 chains. Convergence is checked with $\hat R$. Even so, the model only mixes with the identification choices above: free-sign loadings split the chains across rotation modes ($\hat R\!\sim\!10^{15}$), and a second factor mixes poorly even when identified (see below); one positive-loading factor gives $\hat R\approx1.02$.

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

df = pd.read_csv("etf_returns.csv", parse_dates=["date"], index_col="date")
names = list(df.columns); r = jnp.asarray(df.values - df.values.mean(0)); T, N = r.shape

def ar1(mu, phi, sigma, eps):                              # non-centred AR(1) log-vol
    def step(hp, e): h = mu + phi * (hp - mu) + sigma * e; return h, h
    h0 = mu + sigma / jnp.sqrt(1 - phi ** 2) * eps[0]
    _, hr = jax.lax.scan(step, h0, eps[1:])
    return jnp.concatenate([h0[None], hr], 0)

def model(r):
    M = 1 + N
    mu = numpyro.sample("mu", dist.Normal(0.0, 2.0).expand([M]).to_event(1))
    phi = numpyro.sample("phi", dist.Beta(20.0, 1.5).expand([M]).to_event(1))
    sigma = numpyro.sample("sigma", dist.HalfNormal(0.3).expand([M]).to_event(1))
    eps = numpyro.sample("eps", dist.Normal(0, 1).expand([T, M]).to_event(2))
    h = ar1(mu, phi, sigma, eps)
    df_ = jnp.exp(h[:, 0]); de = jnp.exp(h[:, 1:])         # factor variance, idiosyncratic variances
    bl = numpyro.sample("bl", dist.HalfNormal(1.5).expand([N - 1]).to_event(1))   # positive loadings
    B = jnp.concatenate([jnp.array([1.0]), bl])
    invde = 1.0 / de                                       # --- Woodbury log-likelihood (only a scalar inverse) ---
    Mw = (invde * B ** 2).sum(1) + 1.0 / df_
    br = (invde * B * r).sum(1)
    quad = (r ** 2 * invde).sum(1) - br ** 2 / Mw
    logdetH = jnp.log(de).sum(1) + jnp.log(df_) + jnp.log(Mw)
    numpyro.factor("ll", -0.5 * jnp.sum(N * jnp.log(2 * jnp.pi) + logdetH + quad))

mcmc = MCMC(NUTS(model, target_accept_prob=0.95, init_strategy=numpyro.infer.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(); sg = diag.summary(mcmc.get_samples(group_by_chain=True), prob=0.9)
B_mean = np.concatenate([[1.0], np.array(S["bl"]).mean(0)])
print("Factor SV (1 common factor + %d idiosyncratic), NumPyro NUTS, T=%d:" % (N, T))
print("  max r_hat:", round(max(np.nanmax(sg[p]["r_hat"]) for p in ["mu","phi","sigma","bl"]), 3))
print("  loadings B:", np.round(B_mean, 3), "(", ", ".join(names), ")")
print("  persistence phi: factor %.3f | idio %s" % (np.array(S["phi"]).mean(0)[0], np.round(np.array(S["phi"]).mean(0)[1:], 3)))

Factor SV (1 common factor + 5 idiosyncratic), NumPyro NUTS, T=5283:
  max r_hat: 1.022
  loadings B: [1.    0.869 0.995 0.68  0.538] ( XLF, XLK, XLE, XLV, XLU )
  persistence phi: factor 0.972 | idio [0.99  0.983 0.994 0.963 0.977]
In [2]:
# reconstruct the posterior-mean covariance path H_t = e^{h_f} B B' + diag(e^{h_e}), and read off summaries
@jax.jit
def draw_H(mu, phi, sigma, eps, B):
    h = ar1(mu, phi, sigma, eps); dff = jnp.exp(h[:, 0]); de = jnp.exp(h[:, 1:])
    return dff[:, None, None] * jnp.outer(B, B)[None] + de[:, :, None] * jnp.eye(N)[None], dff, de

ND = len(S["mu"]); Hsum = np.zeros((T, N, N)); dfsum = np.zeros(T)
Bd = np.concatenate([np.ones((ND, 1)), np.array(S["bl"])], 1)
for i in range(ND):
    H, dff, _ = draw_H(jnp.asarray(S["mu"][i]), jnp.asarray(S["phi"][i]), jnp.asarray(S["sigma"][i]),
                       jnp.asarray(S["eps"][i]), jnp.asarray(Bd[i]))
    Hsum += np.array(H); dfsum += np.array(dff)
Hbar = Hsum / ND; np.save("fsv_Hbar.npy", Hbar)            # covariance path -> portfolio horse-race
fac_vol = np.sqrt(252 * dfsum / ND)                        # common-factor vol (annualized %)
di = np.sqrt(np.diagonal(Hbar, axis1=1, axis2=2)); corr = Hbar / (di[:, :, None] * di[:, None, :])
iu = np.triu_indices(N, 1); avg_corr = corr[:, iu[0], iu[1]].mean(1)
common = (B_mean ** 2)[None, :] * (dfsum / ND)[:, None] / np.diagonal(Hbar, axis1=1, axis2=2)
common = common.mean(0)

import matplotlib.pyplot as plt
dates = df.index; crises = {"GFC": "2008-09-15", "COVID": "2020-03-16"}
fig, ax = plt.subplots(1, 3, figsize=(16, 4.3), gridspec_kw=dict(width_ratios=[1.4, 1.4, 1]))
ax[0].plot(dates, fac_vol, color="firebrick", lw=.5)
for lab, dt in crises.items(): ax[0].axvline(pd.Timestamp(dt), color="0.6", lw=.7, ls=":")
ax[0].set_ylabel("annualized volatility (%)"); ax[0].set_title("The common market-vol factor\n(latent SV, spikes in every crisis)")
ax[1].plot(dates, avg_corr, color="steelblue", lw=.3, alpha=.4)
ax[1].plot(dates, pd.Series(avg_corr, index=dates).rolling(21, center=True).mean(), color="navy", lw=1.1)
for lab, dt in crises.items(): ax[1].axvline(pd.Timestamp(dt), color="0.6", lw=.7, ls=":")
ax[1].set_ylim(0, 1); ax[1].set_ylabel("avg pairwise correlation")
ax[1].set_title("Correlation is DRIVEN by the factor\n(rises in crises = diversification meltdown)")
o = np.argsort(-common)
ax[2].bar(range(N), np.array(common)[o], color="firebrick", label="common factor")
ax[2].bar(range(N), 1 - np.array(common)[o], bottom=np.array(common)[o], color="0.8", label="idiosyncratic")
ax[2].set_xticks(range(N)); ax[2].set_xticklabels([names[i] for i in o]); ax[2].set_ylim(0, 1)
ax[2].set_ylabel("share of variance"); ax[2].legend(fontsize=7); ax[2].set_title("Variance decomposition\n(how 'market' each sector is)")
for a in (ax[0], ax[1]):
    for l in a.get_xticklabels(): l.set_rotation(20)
plt.tight_layout(); plt.show()
No description has been provided for this image

Reading the figure¶

  • Left — a latent common-volatility factor. Unlike GARCH, this market-vol process is not a deterministic function of past returns but a hidden state inferred by MCMC. It spikes at every crisis — 2008 (~110% annualised), 2020 (~120%), the 2011 euro crisis, 2018, 2022 — and its persistence is $\phi\approx0.97$.
  • Middle — the factor is the correlation. The model-implied average pairwise correlation swings from ~0.1 to ~0.9, high exactly when the common factor's volatility is high. This is the mechanistic story behind the "diversification meltdown": when market-wide volatility dominates the idiosyncratic terms, $BB'$ dominates $\operatorname{diag}(e^{h^e})$ and everything correlates. DCC described this; factor SV explains it.
  • Right — how "market" each sector is. The variance decomposition (share of each ETF's variance from the common factor) ranges from XLF 68% — financials move with the market — down to XLU 32%, utilities being the most idiosyncratic, defensive sector. The loadings ($b_{\text{XLF}}=1$, $b_{\text{XLE}}=0.99$, $b_{\text{XLU}}=0.54$) are exactly the market betas.

How many factors? One here — and the limits of vanilla NUTS¶

Why a single factor? We also tested a two-factor version (positive-lower-triangular loadings to break the rotation/sign symmetry, same Woodbury likelihood). A modest second factor is plausible — XLU's communality is only 0.32, so one factor leaves real idiosyncratic structure — and factorstochvol's specialised ASIS interweaving sampler is built to handle the multi-factor geometry that defeats plain NUTS.

But plain NUTS will not mix a two-factor model — exploratory runs gave $\hat R$ in the tens to low hundreds, versus $\approx1.02$ for one factor. Multi-factor SV is a notoriously hard HMC target: near-symmetric posterior modes, compounded here by a weakly identified second factor (these five sector ETFs are dominated by one common factor). That contrast is the lesson, not a defect — for a single common factor off-the-shelf NUTS is perfect, but beyond that the model needs a purpose-built sampler (ASIS/interweaving), which is exactly why factorstochvol exists. So we keep the from-scratch model — and the factorsv_R.ipynb cross-check — at one factor (clean and fully identified across both engines); a purpose-built sampler like factorstochvol's is the route to more.

Results — the Bayesian, scalable resolution¶

  • Both p.d. and scalable. Factor SV delivers a guaranteed positive-definite $H_t$ (like BEKK) with only $O(N)$ parameters (like DCC), by routing the co-movement through a low-rank factor. On $N=5$ that is modest; at $N=100$ it is the difference between feasible and hopeless.
  • Natively Bayesian. With latent volatility states there is no two-step QMLE shortcut — MCMC estimates the loadings, the SV parameters, and the entire latent covariance path jointly, returning a full posterior over $H_t$. Convergence is clean ($\hat R\approx1.02$) despite ~32,000 latent states, thanks to marginalisation + Woodbury + the identification/geometry fixes.
  • What it costs. This was the heaviest fit in the project (~44 min): the price of latent-state inference. That is the MGARCH-vs-MSV trade-off in one line — MGARCH is fast and deterministic; MSV is expensive but gives volatility its own honest, uncertainty-quantified dynamics.
  • Next: does it pay off where it matters? The whole point of a covariance model is decisions. The final notebook feeds this posterior $H_t$ into a minimum-variance portfolio and races it, out of sample, against the frequentist DCC — weights and their uncertainty.

References¶

  • Aguilar, O. & West, M. (2000). Bayesian dynamic factor models and portfolio allocation. J. Business & Economic Statistics 18, 338–357.
  • Kastner, G., Frühwirth-Schnatter, S. & Lopes, H. F. (2017). Efficient Bayesian inference for multivariate factor stochastic volatility models. J. Computational & Graphical Statistics 26, 905–917.

Next: factorsv_R.ipynb (the factorstochvol cross-check) and portfolio_horserace.ipynb (the min-variance showdown).