Stochastic volatility from scratch — the state-space view¶

A Kim–Shephard–Chib (1998) mixture sampler, built by hand¶

The other SV notebooks in this arc lean on ready-made engines — sv_numpyro.ipynb (NUTS/JAX), sv_R.ipynb (the stochvol package). This one builds the sampler from scratch to bring out the idea at the centre of modern volatility modelling: stochastic volatility is a state-space model, and the machinery for state-space models — the Kalman filter and a simulation smoother — is exactly what fits it.

The model (returns already de-meaned). Volatility is a latent AR(1) process; returns are noisy readings of it: $$\underbrace{r_t = e^{h_t/2}\,\varepsilon_t,\quad \varepsilon_t\sim N(0,1)}_{\text{measurement equation}}\qquad \underbrace{h_t = \mu + \phi\,(h_{t-1}-\mu) + \sigma_\eta\,\eta_t,\quad \eta_t\sim N(0,1)}_{\text{state (transition) equation}}$$ $h_t$ is the latent log-variance (the state), $\phi$ its persistence, $\sigma_\eta$ the vol-of-vol.

Why "state-space" is the right lens. In the taxonomy of Cox (1981), GARCH is observation-driven — tomorrow's variance is a deterministic function of today's return, so the volatility is known once the data are. SV is parameter-driven — volatility has its own shock $\eta_t$ and is never observed, only inferred. That single change (a latent state with its own innovation) is what makes SV a genuine state-space model, and what makes it harder: we must integrate over the whole hidden path $h_{1:T}$.

The obstacle and the KSC trick. The measurement equation is non-linear in $h_t$. Square and take logs to linearise it: $$y^*_t \equiv \log r_t^2 = h_t + \underbrace{\log\varepsilon_t^2}_{z_t},\qquad z_t\sim\log\chi^2_1\ \ (\text{mean }-1.2704,\ \text{var }\pi^2/2).$$ Now it is linear but the error $z_t$ is non-Gaussian. Kim, Shephard & Chib (1998) approximate $\log\chi^2_1$ by a 7-component Gaussian mixture with indicator $s_t\in\{1,\dots,7\}$. Conditional on the indicators the model is a linear Gaussian state-space model, so the whole path $h_{1:T}$ can be drawn in one block by the Kalman filter + backward sampling. The Gibbs sweep cycles three blocks: $$s\mid h \;\longrightarrow\; h\mid s,\theta\ \ (\text{FFBS})\;\longrightarrow\;\theta=(\mu,\phi,\sigma_\eta)\mid h.$$ All of this is implemented by hand in sv_scratch.py (no PPL, no filtering library).

Step 1 — linearise: the log-χ² mixture¶

The far-left components (with tiny weight) absorb the inliers: days when $r_t\approx 0$ send $\log r_t^2\to-\infty$, and those rare huge-negative values need a home. The mixture reproduces the mean $-1.2704$ and variance $\pi^2/2=4.9348$ of $\log\chi^2_1$ to four decimals — a first internal check that the constants are right.

In [1]:
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import sv_scratch as S

q, m, v2 = S.mixture()                                  # KSC (1998) 7-component mixture for log chi^2_1
mean = np.sum(q * m); var = np.sum(q * (v2 + m**2)) - mean**2
print("KSC mixture moments: mean %.4f (target -1.2704),  var %.4f (target %.4f)" % (mean, var, np.pi**2 / 2))

z = np.linspace(-16, 6, 800)
g_exact = (1/np.sqrt(2*np.pi)) * np.exp(z/2 - np.exp(z)/2)               # exact density of log(chi^2_1)
g_mix = sum(q[i]/np.sqrt(2*np.pi*v2[i])*np.exp(-0.5*(z-m[i])**2/v2[i]) for i in range(7))
fig, ax = plt.subplots(1, 2, figsize=(13, 4.3))
ax[0].plot(z, g_exact, "k", lw=2.2, label=r"exact  $\log\chi^2_1$")
ax[0].plot(z, g_mix, "--", color="firebrick", lw=1.8, label="KSC 7-component mixture")
for i in range(7):
    ax[0].plot(z, q[i]/np.sqrt(2*np.pi*v2[i])*np.exp(-0.5*(z-m[i])**2/v2[i]), color="steelblue", lw=.7, alpha=.6)
ax[0].set_title("Linearising the measurement equation\n"
                r"$\log\varepsilon_t^2$ approximated by a 7-normal mixture (KSC 1998)")
ax[0].set_xlabel(r"$z=\log\varepsilon^2$"); ax[0].set_ylabel("density"); ax[0].legend(fontsize=8)
ax[0].text(-15.2, .01, "far-left components\ncapture inliers ($r_t\\approx0$)", fontsize=7, color="dimgray")
rr, hh = S.simulate_sv(400, -0.5, 0.98, 0.15, seed=7)
ax[1].plot(rr, color="0.6", lw=.7, label="returns  $r_t$")
ax[1].plot(np.exp(hh/2), color="firebrick", lw=1.6, label=r"latent vol  $e^{h_t/2}$")
ax[1].plot(-np.exp(hh/2), color="firebrick", lw=1.6)
ax[1].set_title("SV is a state-space model:\nreturns are noisy readings of a latent volatility state")
ax[1].set_xlabel("t"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()

KSC mixture moments: mean -1.2704 (target -1.2704),  var 4.9349 (target 4.9348)
No description has been provided for this image

Step 2 — draw the whole volatility path: Kalman filter + FFBS¶

Conditional on the mixture indicators $s_{1:T}$, the system is linear Gaussian: $$y^*_t = h_t + \big(m_{s_t}-1.2704\big) + \sqrt{v_{s_t}}\,u_t,\qquad h_t = (1-\phi)\mu + \phi h_{t-1} + \sigma_\eta\eta_t .$$ The Kalman filter runs forward computing the filtered mean/variance of $h_t\mid y^*_{1:t}$; the backward pass then samples $h_T, h_{T-1},\dots,h_1$ from the exact joint smoothing distribution — the forward-filter backward-sample (FFBS) simulation smoother of Carter & Kohn (1994) and Frühwirth-Schnatter (1994). Because it draws the entire path at once (a "multi-move" step) it mixes far better than updating each $h_t$ singly, which is the whole reason KSC beat the earlier single-move sampler of Jacquier, Polson & Rossi (1994). The parameter block is standard: $\sigma_\eta^2$ inverse-gamma, $\mu$ Gaussian, and $\phi$ by an independence-Metropolis step with a $\text{Beta}$ prior on $(\phi+1)/2$ (keeping $|\phi|<1$).

Simulation check. Before touching real data, simulate $T=2000$ points with known $(\mu,\phi,\sigma_\eta)=(-0.5,0.98,0.15)$ and confirm the sampler recovers them and the hidden path.

In [2]:
mu_t, phi_t, sig_t = -0.5, 0.98, 0.15
r_sim, h_true = S.simulate_sv(2000, mu_t, phi_t, sig_t, seed=1)
Ds = S.sv_gibbs(r_sim, niter=4000, burn=1000, seed=2)
for nm, tru in [("mu", mu_t), ("phi", phi_t), ("sigma", sig_t)]:
    x = Ds[nm]
    print("%-6s posterior %7.4f  [%.4f, %.4f]   true %.3f" % (nm, x.mean(), np.quantile(x, .025), np.quantile(x, .975), tru))
print("corr(posterior mean h, true h) = %.3f" % np.corrcoef(Ds['h_mean'], h_true)[0, 1])

fig, ax = plt.subplots(1, 3, figsize=(14, 3.9))
for a, nm, tru, lab in [(ax[0], "phi", phi_t, r"$\phi$"), (ax[1], "sigma", sig_t, r"$\sigma_\eta$")]:
    a.hist(Ds[nm], bins=45, color="steelblue", alpha=.8, density=True)
    a.axvline(tru, color="firebrick", lw=2, ls="--", label="true = %.2f" % tru)
    a.set_title("posterior of %s" % lab); a.legend(fontsize=8); a.set_yticks([])
w = slice(600, 900)
ax[2].plot(np.arange(2000)[w], h_true[w], color="firebrick", lw=1.8, label="true $h_t$")
ax[2].plot(np.arange(2000)[w], Ds['h_mean'][w], color="steelblue", lw=1.3, label="posterior mean")
ax[2].fill_between(np.arange(2000)[w], (Ds['h_mean']-2*Ds['h_sd'])[w], (Ds['h_mean']+2*Ds['h_sd'])[w],
                   color="steelblue", alpha=.2)
ax[2].set_title("latent log-variance recovered\n(corr = %.2f)" % np.corrcoef(Ds['h_mean'], h_true)[0, 1])
ax[2].set_xlabel("t"); ax[2].legend(fontsize=8)
plt.suptitle("Simulation check: the from-scratch KSC sampler recovers known parameters and the latent path", y=1.03)
plt.tight_layout(); plt.show()

mu     posterior -0.5623  [-1.0203, -0.0381]   true -0.500
phi    posterior  0.9852  [0.9754, 0.9938]   true 0.980
sigma  posterior  0.1510  [0.1215, 0.1834]   true 0.150
corr(posterior mean h, true h) = 0.923
No description has been provided for this image

Step 3 — fit the S&P 500 and cross-check the engines¶

Same sp500ret.csv used throughout the arc (5,523 daily returns, 1987–2009). The from-scratch posterior should land on the same values as the NUTS fit (sv_numpyro) and the stochvol package (sv_R) — a three-engine validation.

In [3]:
df = pd.read_csv("sp500ret.csv")
ret = df["ret"].values.astype(float); dates = pd.to_datetime(df["date"])
D = S.sv_gibbs(ret, niter=6000, burn=1000, seed=3)                       # ~90s (pure-numpy Kalman + FFBS)
print("From-scratch KSC posterior (S&P 500):")
for nm in ["mu", "phi", "sigma"]:
    x = D[nm]; print("  %-6s %7.4f   [%.4f, %.4f]" % (nm, x.mean(), np.quantile(x, .025), np.quantile(x, .975)))
print("cross-check  ->  NumPyro/NUTS: phi 0.987, sigma 0.155, mu -0.23   |   stochvol: phi 0.987, sigma 0.152, mu -0.231")

vol = np.exp(D['h_mean'] / 2)
lo = np.exp(np.quantile(D['h_paths'], .05, axis=0) / 2)
hi = np.exp(np.quantile(D['h_paths'], .95, axis=0) / 2)
fig = plt.figure(figsize=(13, 7.4))
gs = fig.add_gridspec(2, 2, height_ratios=[2, 1], hspace=.34, wspace=.18)
axt = fig.add_subplot(gs[0, :])
axt.plot(dates, np.abs(D['r']), color="0.75", lw=.4, label="|return|")
axt.plot(dates, vol, color="firebrick", lw=1.0, label="SV posterior vol  $e^{h_t/2}$")
axt.fill_between(dates, lo, hi, color="firebrick", alpha=.25, label="90% credible band")
axt.set_ylim(0, 8); axt.set_title("From-scratch KSC stochastic-volatility fit — S&P 500 daily returns 1987–2009")
axt.set_ylabel("daily volatility (%)"); axt.legend(fontsize=8, loc="upper left")
ax1 = fig.add_subplot(gs[1, 0])
ax1.hist(D['phi'], bins=50, color="steelblue", alpha=.8, density=True)
ax1.axvline(0.9865, color="firebrick", ls="--", lw=1.6, label="NumPyro 0.987 /\nstochvol 0.987")
ax1.set_yticks([]); ax1.legend(fontsize=8); ax1.set_xlabel(r"$\phi$  (persistence)"); ax1.set_title(r"posterior of $\phi$")
ax2 = fig.add_subplot(gs[1, 1])
ax2.hist(D['sigma'], bins=50, color="seagreen", alpha=.75, density=True)
ax2.axvline(0.1533, color="firebrick", ls="--", lw=1.6, label="NumPyro 0.155 /\nstochvol 0.152")
ax2.set_yticks([]); ax2.legend(fontsize=8); ax2.set_xlabel(r"$\sigma_\eta$  (vol-of-vol)"); ax2.set_title(r"posterior of $\sigma_\eta$")
plt.show()

From-scratch KSC posterior (S&P 500):
  mu     -0.2353   [-0.5766, 0.1131]
  phi     0.9879   [0.9814, 0.9931]
  sigma   0.1468   [0.1226, 0.1742]
cross-check  ->  NumPyro/NUTS: phi 0.987, sigma 0.155, mu -0.23   |   stochvol: phi 0.987, sigma 0.152, mu -0.231
No description has been provided for this image

Results¶

  • The from-scratch sampler is correct. It recovers known parameters on simulated data (Step 2) and, on the S&P 500, lands on $\phi=0.988$, $\sigma_\eta=0.147$, $\mu=-0.235$ — matching the NUTS fit (sv_numpyro: 0.987 / 0.155 / −0.23) and the stochvol package (sv_R: 0.987 / 0.152 / −0.231) to well within the credible intervals. Three independent engines, one posterior — the from-scratch $\sigma_\eta$ sits ~5% below the exact samplers (0.147 vs 0.152–0.155), the small downward bias expected from the KSC 7-component approximation to the $\log\chi^2$ errors.
  • Volatility is persistent but its own random process. $\phi\approx0.99$ (highly persistent, but stationary — below the near-unit-root $\alpha+\beta\approx0.997$ that GARCH reports on the same data), and $\sigma_\eta\approx0.15$ is a genuine volatility-of-volatility: unlike GARCH, the fitted vol carries a posterior credible band (widening through the 1987 crash and 2008), because $h_t$ is a latent state we infer rather than a deterministic function of past returns.
  • The pay-off of the state-space view. Writing SV as a linear Gaussian state-space model (conditional on the mixture indicators) turns an intractable non-linear filtering problem into a Kalman-filter + FFBS block that draws the entire 5,523-point volatility path at once. That is the KSC insight, and doing it by hand shows there is no black box — just a filter, a smoother, and three Gibbs blocks.

References¶

  • Taylor, S. J. (1982, 1986). The original stochastic-volatility model.
  • Cox, D. R. (1981). Statistical analysis of time series. Scand. J. Statist. — the observation- vs parameter-driven distinction.
  • Jacquier, E., Polson, N. G. & Rossi, P. E. (1994). Bayesian analysis of stochastic volatility models. JBES — the single-move Gibbs sampler.
  • Carter, C. K. & Kohn, R. (1994). On Gibbs sampling for state space models. Biometrika; Frühwirth-Schnatter, S. (1994). Data augmentation and dynamic linear models. JTSA — the FFBS simulation smoother.
  • Kim, S., Shephard, N. & Chib, S. (1998). Stochastic volatility: likelihood inference and comparison with ARCH models. Review of Economic Studies — the mixture / multi-move sampler implemented here.
  • Durbin, J. & Koopman, S. J. (2012). Time Series Analysis by State Space Methods (2nd ed.).
  • Omori, Y., Chib, S., Shephard, N. & Nakajima, J. (2007). Stochastic volatility with leverage. J. Econometrics — the 10-component refinement with leverage.
  • Kastner, G. & Frühwirth-Schnatter, S. (2014). Ancillarity-sufficiency interweaving (ASIS). CSDA — the algorithm behind the stochvol package used in sv_R.

Next in the SV arc: sv_ext_numpyro adds Student-$t$ innovations and leverage (the SV analogue of GJR); a head-to-head sv_vs_garch compares rolling predictive scores and VaR/ES.