Stochastic Volatility — Bayesian (NumPyro)¶

The volatility is a latent process, not a deterministic recursion¶

Phase 3 of the volatility arc. In GARCH the conditional variance $\sigma_t^2$ is a deterministic function of past returns; in a stochastic-volatility (SV) model the log-variance is its own latent AR(1) driven by an independent shock:

$$r_t = e^{h_t/2}\,\varepsilon_t,\qquad h_t = \mu + \phi\,(h_{t-1}-\mu) + \sigma_\eta\,\eta_t,\qquad \varepsilon_t,\eta_t\sim N(0,1).$$

  • $\phi$ — persistence of log-volatility (≈ the GARCH $\alpha+\beta$ analogue), $\sigma_\eta$ — vol-of-vol, $\mu$ — mean log-variance.
  • Because $h_t$ has its own innovation $\eta_t$, the volatility path carries posterior uncertainty — a credible band — which a GARCH filter (vol known given the past) cannot produce. That is the headline difference.

Engine: NumPyro/NUTS with the non-centered parameterization (sampling standardized innovations $z_t$ and rebuilding $h_t$ via a JAX scan) — essential to defuse the SV funnel. On this machine PyMC is slow on scan models, so NumPyro is the Python SV engine; sv_R.ipynb will validate with Kastner's stochvol. Data: daily S&P 500 returns, 1987–2009.

In [1]:
import os
_lib = r"C:\Users\user\anaconda3\envs\pymc-env\Library\bin"
if _lib not in os.environ.get("PATH", ""): os.environ["PATH"] = _lib + os.pathsep + os.environ.get("PATH", "")
import numpy as np, pandas as pd, jax
jax.config.update("jax_enable_x64", True)
import jax.numpy as jnp, numpyro, numpyro.distributions as dist, matplotlib.pyplot as plt
from numpyro.infer import MCMC, NUTS
from jax import lax
numpyro.set_host_device_count(1)

d = pd.read_csv("sp500ret.csv", parse_dates=["date"]); dates = d["date"].values
y = d["ret"].values - d["ret"].values.mean()                # demeaned daily % returns
print("%d days %s..%s   sample vol %.2f%%/day" % (len(y), str(dates[0])[:10], str(dates[-1])[:10], np.std(y)))

5523 days 1987-03-10..2009-01-30   sample vol 1.19%/day
In [2]:
# --- Bayesian SV: latent log-variance AR(1), non-centered (defuses the funnel), sampled by NUTS ---
def sv_model(y):
    T     = y.shape[0]
    mu    = numpyro.sample("mu", dist.Normal(0., 5.))                                     # mean log-variance
    phi   = numpyro.deterministic("phi", 2*numpyro.sample("phi01", dist.Beta(20., 1.5)) - 1)   # persistence in (-1,1)
    sigma = numpyro.sample("sigma", dist.HalfNormal(0.3))                                 # vol-of-vol
    z     = numpyro.sample("z", dist.Normal(0., 1.).expand([T]))                          # standardized innovations
    h0    = mu + sigma/jnp.sqrt(1 - phi**2)*z[0]                                          # stationary initial log-var
    def step(h_prev, z_t):
        h = mu + phi*(h_prev - mu) + sigma*z_t
        return h, h
    _, h_rest = lax.scan(step, h0, z[1:])
    h = numpyro.deterministic("h", jnp.concatenate([h0[None], h_rest]))                   # latent log-variance path
    numpyro.sample("obs", dist.Normal(0., jnp.exp(h/2)), obs=y)                           # r_t = exp(h_t/2) * eps_t

mcmc = MCMC(NUTS(sv_model, target_accept_prob=0.9), num_warmup=1000, num_samples=1000, num_chains=2, progress_bar=False)
mcmc.run(jax.random.PRNGKey(1), y)
ps = mcmc.get_samples(group_by_chain=False)
psc = mcmc.get_samples(group_by_chain=True)
rhat = lambda k: float(numpyro.diagnostics.gelman_rubin(psc[k]))
print("posterior means:  phi %.4f   sigma %.4f   mu %.4f" % (ps["phi"].mean(), ps["sigma"].mean(), ps["mu"].mean()))
print("r_hat:            phi %.3f   sigma %.3f   mu %.3f" % (rhat("phi"), rhat("sigma"), rhat("mu")))
print("95%% CrI  phi [%.4f, %.4f]   sigma [%.3f, %.3f]" %
      (*np.percentile(ps["phi"], [2.5, 97.5]), *np.percentile(ps["sigma"], [2.5, 97.5])))

C:\Users\user\AppData\Local\Temp\ipykernel_21484\3821468155.py:16: UserWarning: There are not enough devices to run parallel chains: expected 2 but got 1. Chains will be drawn sequentially. If you are running MCMC in CPU, consider using `numpyro.set_host_device_count(2)` at the beginning of your program. You can double-check how many devices are available in your system using `jax.local_device_count()`.
  mcmc = MCMC(NUTS(sv_model, target_accept_prob=0.9), num_warmup=1000, num_samples=1000, num_chains=2, progress_bar=False)

posterior means:  phi 0.9865   sigma 0.1546   mu -0.2334
r_hat:            phi 1.000   sigma 1.000   mu 1.003
95% CrI  phi [0.9804, 0.9922]   sigma [0.132, 0.179]
In [3]:
# GARCH(1,1)-t on the same series, for the deterministic-vol contrast
from arch import arch_model
gm = arch_model(y, mean="Zero", vol="GARCH", p=1, q=1, dist="t").fit(disp="off")
garch_vol = np.asarray(gm.conditional_volatility)
gp = gm.params
print("GARCH(1,1)-t: omega %.4f  alpha %.3f  beta %.3f  persistence %.4f  nu %.1f" %
      (gp["omega"], gp["alpha[1]"], gp["beta[1]"], gp["alpha[1]"]+gp["beta[1]"], gp["nu"]))
print("SV persistence phi %.4f  (vs GARCH alpha+beta %.4f)" % (ps["phi"].mean(), gp["alpha[1]"]+gp["beta[1]"]))

GARCH(1,1)-t: omega 0.0061  alpha 0.061  beta 0.936  persistence 0.9969  nu 6.2
SV persistence phi 0.9865  (vs GARCH alpha+beta 0.9969)
In [4]:
# HEADLINE FIGURE: SV latent volatility (with a 90% credible BAND) vs GARCH's single deterministic line
vol = np.exp(ps["h"]/2)                                        # (draws, T) posterior volatility paths
vlo, vmed, vhi = np.percentile(vol, [5, 50, 95], axis=0)
dt = pd.to_datetime(dates)
fig, ax = plt.subplots(2, 1, figsize=(12, 7.2))
ax[0].plot(dt, np.abs(y), color="0.8", lw=.35, label="|return|")
ax[0].fill_between(dt, vlo, vhi, color="steelblue", alpha=.30, label="SV vol — 90% credible band")
ax[0].plot(dt, vmed, color="steelblue", lw=.9, label="SV vol (posterior median)")
ax[0].plot(dt, garch_vol, color="firebrick", lw=.7, label="GARCH(1,1)-t vol (deterministic)")
ax[0].set_ylabel("daily volatility (%)"); ax[0].legend(fontsize=8, ncol=2, loc="upper left")
ax[0].set_title("Stochastic volatility carries UNCERTAINTY (a band); GARCH gives a single line")
m = (dt >= "2008-06-01") & (dt <= "2009-06-30")
ax[1].fill_between(dt[m], vlo[m], vhi[m], color="steelblue", alpha=.30, label="SV vol — 90% band")
ax[1].plot(dt[m], vmed[m], color="steelblue", lw=1.4, label="SV vol (median)")
ax[1].plot(dt[m], garch_vol[m], color="firebrick", lw=1.1, label="GARCH vol")
ax[1].set_ylabel("vol (%)"); ax[1].set_title("2008 crisis zoom — the SV band widens exactly when vol is hardest to pin down")
ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig("sv_numpyro_vol.png", dpi=120, bbox_inches="tight"); plt.show()
No description has been provided for this image
In [5]:
# SV parameter posteriors: persistence and vol-of-vol
fig, ax = plt.subplots(1, 2, figsize=(11, 4))
ax[0].hist(ps["phi"], bins=45, color="steelblue", density=True, alpha=.8)
ax[0].axvline(np.median(ps["phi"]), color="k", lw=1)
ax[0].axvline(gp["alpha[1]"]+gp["beta[1]"], color="firebrick", ls="--", lw=1.5, label="GARCH $\\alpha+\\beta$")
ax[0].set_xlabel("$\\phi$  (log-vol persistence)"); ax[0].set_ylabel("posterior density")
ax[0].set_title("Persistence $\\phi$ = %.3f" % ps["phi"].mean()); ax[0].legend(fontsize=8)
ax[1].hist(ps["sigma"], bins=45, color="seagreen", density=True, alpha=.8)
ax[1].axvline(np.median(ps["sigma"]), color="k", lw=1)
ax[1].set_xlabel("$\\sigma_\\eta$  (vol-of-vol)")
ax[1].set_title("Vol-of-vol $\\sigma_\\eta$ = %.3f  (>0 => genuinely stochastic vol)" % ps["sigma"].mean())
plt.tight_layout(); plt.savefig("sv_numpyro_params.png", dpi=120, bbox_inches="tight"); plt.show()
No description has been provided for this image

Results¶

  • Persistence $\phi \approx 0.987$ — log-volatility is highly persistent but stationary, a touch below the GARCH $\alpha+\beta \approx 0.997$ (the two aren't the same parameter, but both say "shocks to vol die out very slowly").
  • Vol-of-vol $\sigma_\eta \approx 0.15$, clearly bounded away from 0 — the volatility is genuinely stochastic, not a deterministic function of returns. That is the qualitative break from GARCH.
  • The band is the point. SV's volatility path comes with a 90% credible band that widens exactly in turbulent stretches (2008): the model is honest that vol is estimated, not observed. GARCH returns a single line with no such uncertainty.
  • Median and crisis-peak vol track GARCH closely, so the two agree on the level of volatility while disagreeing on its epistemics.

Next in the arc: sv_R.ipynb (Kastner stochvol, two-engine validation), then SV-t and leverage SV ($\text{corr}(\varepsilon,\eta)<0$, the SV analogue of GJR), and the honest sv_vs_garch head-to-head (rolling one-step log-predictive score + VaR/ES backtest).