Stochastic Volatility — extensions: fat tails & leverage (NumPyro)¶

SV-t and the asymmetric (leverage) SV — the SV analogues of GARCH-t and GJR¶

Two ingredients on top of the plain SV of sv_numpyro.ipynb, fit as a $2\times2$ family:

  • Student-t innovations (SV-t): $\varepsilon_t\sim t_\nu$ (standardized so $\mathrm{Var}=e^{h_t}$). Even with stochastic vol, the conditional return distribution can be fat-tailed — the SV analogue of GARCH-t.
  • Leverage (asymmetric SV): correlate the return and log-vol shocks, $\mathrm{corr}(\varepsilon_t,\eta_t)=\rho$. A negative return then pushes next-period volatility up — $\rho<0$ is the SV analogue of the GJR $\gamma>0$. In the non-centered scan the log-vol innovation gets a $\rho\,\varepsilon_{t-1}$ term, with $\varepsilon_{t-1}=r_{t-1}e^{-h_{t-1}/2}$.

$$r_t=e^{h_t/2}\varepsilon_t,\quad h_t=\mu+\phi(h_{t-1}-\mu)+\sigma_\eta\big(\rho\,\varepsilon_{t-1}+\sqrt{1-\rho^2}\,z_t\big).$$

We fit all four (SV, SV-t, SV-leverage, SV-t-leverage) and rank them by WAIC.

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, log_likelihood
from jax import lax
numpyro.set_host_device_count(1)
d = pd.read_csv("sp500ret.csv"); y = d["ret"].values - d["ret"].values.mean()

def sv_model(y, student=False, leverage=False):
    T = y.shape[0]
    mu    = numpyro.sample("mu", dist.Normal(0., 5.))
    phi   = numpyro.deterministic("phi", 2*numpyro.sample("phi01", dist.Beta(20., 1.5)) - 1)
    sigma = numpyro.sample("sigma", dist.HalfNormal(0.3))
    rho = numpyro.sample("rho", dist.Uniform(-0.999, 0.999)) if leverage else 0.0
    nu  = numpyro.sample("nu", dist.Uniform(2.5, 40.)) if student else 0.0
    z   = numpyro.sample("z", dist.Normal(0., 1.).expand([T]))
    h0  = mu + sigma/jnp.sqrt(1 - phi**2)*z[0]
    if leverage:
        def step(hp, inp):
            zt, rp = inp
            eps = rp*jnp.exp(-hp/2)                                   # standardized past return
            return (mu + phi*(hp - mu) + sigma*(rho*eps + jnp.sqrt(1 - rho**2)*zt),)*2
        _, hr = lax.scan(step, h0, (z[1:], y[:-1]))
    else:
        def step(hp, zt):
            return (mu + phi*(hp - mu) + sigma*zt,)*2
        _, hr = lax.scan(step, h0, z[1:])
    h = numpyro.deterministic("h", jnp.concatenate([h0[None], hr]))
    scale = jnp.exp(h/2)
    if student:
        numpyro.sample("obs", dist.StudentT(nu, 0., scale*jnp.sqrt((nu - 2)/nu)), obs=y)
    else:
        numpyro.sample("obs", dist.Normal(0., scale), obs=y)

def fit_waic(student, leverage, key, draws=800):
    m = MCMC(NUTS(sv_model, target_accept_prob=0.9), num_warmup=draws, num_samples=draws, num_chains=2, progress_bar=False)
    m.run(jax.random.PRNGKey(key), y, student=student, leverage=leverage)
    s = m.get_samples()
    ll = np.asarray(log_likelihood(sv_model, s, y, student=student, leverage=leverage)["obs"])   # (draws, T)
    mx = ll.max(0); lppd = (mx + np.log(np.exp(ll - mx).mean(0))).sum(); p_waic = ll.var(0).sum()
    return s, -2*(lppd - p_waic)
In [2]:
variants = [("SV", False, False), ("SV-t", True, False), ("SV-leverage", False, True), ("SV-t-leverage", True, True)]
R = {}
for i, (tag, st, lv) in enumerate(variants):
    s, w = fit_waic(st, lv, i + 1)
    R[tag] = dict(s=s, waic=w, phi=s["phi"].mean(), sigma=s["sigma"].mean(),
                  rho=(s["rho"].mean() if lv else np.nan), nu=(s["nu"].mean() if st else np.nan))
    print("%-15s WAIC %8.0f   phi %.3f  sigma %.3f  rho %6s  nu %5s" %
          (tag, w, R[tag]["phi"], R[tag]["sigma"],
           ("%.3f" % R[tag]["rho"]) if lv else "-", ("%.1f" % R[tag]["nu"]) if st else "-"))
best = min(R, key=lambda t: R[t]["waic"])
print("\nbest by WAIC:", best)

C:\Users\user\AppData\Local\Temp\ipykernel_53100\2216044468.py:39: 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()`.
  m = MCMC(NUTS(sv_model, target_accept_prob=0.9), num_warmup=draws, num_samples=draws, num_chains=2, progress_bar=False)

SV              WAIC    14592   phi 0.986  sigma 0.154  rho      -  nu     -

SV-t            WAIC    14560   phi 0.994  sigma 0.105  rho      -  nu   9.1

SV-leverage     WAIC    14513   phi 0.980  sigma 0.183  rho -0.608  nu     -

SV-t-leverage   WAIC    14460   phi 0.987  sigma 0.142  rho -0.727  nu   9.2

best by WAIC: SV-t-leverage
In [3]:
fig, ax = plt.subplots(1, 3, figsize=(15, 4.3))
# (A) leverage rho << 0
for tag, col in [("SV-leverage", "steelblue"), ("SV-t-leverage", "firebrick")]:
    ax[0].hist(R[tag]["s"]["rho"], bins=45, density=True, alpha=.55, color=col, label="%s (rho %.2f)" % (tag, R[tag]["rho"]))
ax[0].axvline(0, color="k", lw=1)
ax[0].set_xlabel(r"$\rho$  (return / vol-shock correlation)"); ax[0].set_ylabel("posterior density")
ax[0].set_title("Leverage: $\\rho \\ll 0$\n(bad news raises future vol = SV analogue of GJR)"); ax[0].legend(fontsize=8)
# (B) fat tails: nu
for tag, col in [("SV-t", "seagreen"), ("SV-t-leverage", "firebrick")]:
    ax[1].hist(R[tag]["s"]["nu"], bins=45, density=True, alpha=.55, color=col, label="%s (nu %.1f)" % (tag, R[tag]["nu"]))
ax[1].set_xlabel(r"$\nu$  (Student-t d.o.f.)"); ax[1].set_title("Fat tails: finite $\\nu$\n(extra kurtosis beyond stochastic vol)"); ax[1].legend(fontsize=8)
# (C) WAIC comparison
tags = list(R); waics = [R[t]["waic"] for t in tags]; bi = int(np.argmin(waics))
ax[2].barh(range(len(tags)), waics, color=["firebrick" if i == bi else "steelblue" for i in range(len(tags))])
ax[2].set_yticks(range(len(tags))); ax[2].set_yticklabels(tags); ax[2].invert_yaxis()
ax[2].set_xlim(min(waics) - 40, max(waics) + 40)
for i, w in enumerate(waics): ax[2].text(w, i, " %.0f" % w, va="center", fontsize=8)
ax[2].set_xlabel("WAIC (lower = better)"); ax[2].set_title("Model comparison — WAIC")
plt.tight_layout(); plt.savefig("sv_ext_numpyro.png", dpi=120, bbox_inches="tight"); plt.show()
No description has been provided for this image

Results¶

  • Leverage is strong and unambiguous: $\rho \approx -0.6$ to $-0.7$, entirely below 0 — a negative return today drives volatility up tomorrow. This is the stochastic-volatility counterpart of the GJR $\gamma>0$ leverage term, and the single most important extension.
  • Fat tails persist even with stochastic vol: $\nu \approx 8$–$12$ — the conditional return is still leptokurtic on top of the vol dynamics (as in GARCH-t, $\nu\approx6$).
  • WAIC ranking: adding leverage buys the largest improvement; t adds a further, smaller gain, so SV-t-leverage is preferred — the full asymmetric, fat-tailed SV. (The exact gaps are printed above.)

Cross-model picture so far: GARCH said persistence $\alpha+\beta\approx0.997$, $\nu\approx6$, GJR $\gamma>0$; SV says $\phi\approx0.98$, $\nu\approx9$, $\rho<0$ — the same three stylized facts (persistence, fat tails, leverage) seen through a latent-volatility lens.

Next: sv_ext_R.ipynb — stochvol::svtsample / svlsample to validate $\nu$ and $\rho$; then the honest sv_vs_garch head-to-head (rolling one-step log-predictive score + VaR/ES backtest), which also becomes the explorer's #1 panel.