SV vs GARCH — which volatility model forecasts better?¶
Out-of-sample 1-step density forecasts + Value-at-Risk backtest¶
The volatility arc has fit two paradigms on the same S&P 500 series:
- GARCH (observation-driven): tomorrow's variance is a deterministic function of today's return, $\sigma_t^2=\omega+\alpha r_{t-1}^2+\beta\sigma_{t-1}^2$ — known once the data are.
- Stochastic volatility (parameter-driven): volatility is a latent AR(1) state with its own shock, inferred rather than observed.
This notebook settles the practical question — does the extra flexibility of a latent-state volatility actually forecast better? — with a genuine out-of-sample horse race.
Design (fixed-parameter rolling filter). Estimate each model once on a training window (1987–2001), then roll the volatility filter forward across the test set (2001–2009, including the 2008 crisis) producing true 1-step-ahead predictive densities $p(r_t\mid r_{1:t-1})$. Two scores:
- Log predictive score (LPS) — the mean 1-step log predictive density, a proper scoring rule (higher = better).
- Value-at-Risk backtest — Kupiec (1995) unconditional coverage and Christoffersen (1998) conditional coverage, plus expected shortfall, at the 1% and 5% levels.
A fair fight. To avoid comparing a fully-loaded GARCH against a bare-bones SV, each paradigm gets a plain and a fully-featured member:
- GARCH-N (Gaussian) and GARCH-t and GJR-t (fat tails + leverage) — from
garch_scratch. - SV (Gaussian) and SV-t-leverage (fat tails + leverage) — the SV analogue of GJR-t.
A structural difference shows up immediately. GARCH's 1-step predictive is closed form. SV's is not — because SV is a non-linear, non-Gaussian state-space model, its predictive must be obtained by a bootstrap particle filter (svgarch_forecast), and it emerges as a scale-mixture of normals/t — heavier-tailed than any single density, for free.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import garch_scratch as G, sv_scratch as S, svgarch_forecast as F
df = pd.read_csv("sp500ret.csv"); ret = df["ret"].values.astype(float); dates = pd.to_datetime(df["date"])
ret = ret - ret.mean()
n = len(ret); ntr = 3500 # train = first 3500 days (1987-2001)
train, test_idx = ret[:ntr], np.arange(ntr, n)
mask = np.zeros(n, bool); mask[ntr:] = True
alphas = (0.01, 0.05)
print("train %d (%s..%s) test %d (%s..%s)" % (ntr, dates[0].date(), dates[ntr-1].date(),
n-ntr, dates[ntr].date(), dates[n-1].date()))
# --- estimate each model once on the training window ---
th_gn, _, _ = G.garch_mle(train, gjr=False, student=False)
th_gt, _, _ = G.garch_mle(train, gjr=False, student=True)
th_jt, _, _ = G.garch_mle(train, gjr=True, student=True)
s0 = float(np.var(train))
Dsv = S.sv_gibbs(train, niter=4000, burn=1000, seed=3) # Gaussian SV (KSC sampler)
mu, phi, sig = Dsv["mu"].mean(), Dsv["phi"].mean(), Dsv["sigma"].mean()
mu2, phi2, sig2, nu2, rho2 = F.sv_fit_pf(train, N=2000, seed=0) # SV-t-leverage (particle-filter ML, ~5 min)
print("estimated on train:")
print(" GARCH-N omega %.4f alpha %.3f beta %.3f" % (th_gn[0], th_gn[1], th_gn[2]))
print(" GARCH-t omega %.4f alpha %.3f beta %.3f nu %.1f" % (th_gt[0], th_gt[1], th_gt[2], th_gt[3]))
print(" GJR-t omega %.4f alpha %.3f beta %.3f gamma %.3f nu %.1f" % (th_jt[0], th_jt[1], th_jt[2], th_jt[3], th_jt[4]))
print(" SV mu %.2f phi %.3f sigma_eta %.3f" % (mu, phi, sig))
print(" SV-t-lev mu %.2f phi %.3f sigma %.3f nu %.1f rho %+.2f" % (mu2, phi2, sig2, nu2, rho2))
train 3500 (1987-03-10..2001-01-11) test 2023 (2001-01-12..2009-01-30) estimated on train: GARCH-N omega 0.0159 alpha 0.099 beta 0.893 GARCH-t omega 0.0061 alpha 0.051 beta 0.944 nu 5.2 GJR-t omega 0.0127 alpha 0.021 beta 0.924 gamma 0.089 nu 5.5 SV mu -0.37 phi 0.980 sigma_eta 0.174 SV-t-lev mu 0.12 phi 0.988 sigma 0.134 nu 6.9 rho -0.49
Rolling the filters forward — closed form vs particle filter¶
For GARCH the variance recursion is deterministic, so garch_predict runs it forward with the fixed parameters and reads off a scaled Student-$t$ (or Normal) predictive at each step.
For SV the latent log-variance $h_t$ is never observed, so sv_particle_filter runs a bootstrap particle filter: carry $N$ particles $\{h_t^{(i)}\}$, (i) propagate each through the state equation — for SV-t-leverage the propagation shock is correlated with the last return, $h_t=\mu+\phi(h_{t-1}-\mu)+\sigma(\rho\,u_{t-1}+\sqrt{1-\rho^2}\,\zeta)$ with $u_{t-1}=r_{t-1}e^{-h_{t-1}/2}$ — (ii) form the 1-step predictive as the average $\frac1N\sum_i f\!\big(r_t;e^{h_t^{(i)}}\big)$, a scale-mixture, (iii) weight by that likelihood and resample. The SV-$t$-leverage parameters were themselves estimated on the training window by maximizing this particle filter's likelihood (sv_fit_pf) — the recovered $\nu\approx7,\ \rho\approx-0.49$ match the fat tails and leverage found in sv_ext_R.
lps_gn, var_gn, _ = F.garch_predict(ret, th_gn, s0, gjr=False, alphas=alphas)
lps_gt, var_gt, sd_gt = F.garch_predict(ret, th_gt, s0, gjr=False, alphas=alphas)
lps_jt, var_jt, sd_jt = F.garch_predict(ret, th_jt, s0, gjr=True, alphas=alphas)
lps_sv, var_sv, svvol = F.sv_particle_filter(ret, mu, phi, sig, N=3000, seed=1, alphas=alphas, var_mask=mask)
lps_stl, var_stl, stlvol = F.sv_particle_filter(ret, mu2, phi2, sig2, nu=nu2, rho=rho2, N=3000, seed=1, alphas=alphas, var_mask=mask)
models = {"GARCH-N": (lps_gn, var_gn), "GARCH-t": (lps_gt, var_gt), "GJR-t": (lps_jt, var_jt),
"SV": (lps_sv, var_sv), "SV-t-lev": (lps_stl, var_stl)}
print("Mean 1-step log predictive score (test set, higher = better):")
for nm, (lps, _) in models.items():
print(" %-9s % .4f" % (nm, lps[test_idx].mean()))
Mean 1-step log predictive score (test set, higher = better): GARCH-N -1.4438 GARCH-t -1.4390 GJR-t -1.4238 SV -1.4376 SV-t-lev -1.4274
Reading the score. The ranking is GJR-t $>$ SV-t-lev $>$ SV $>$ GARCH-t $>$ GARCH-N. Two clean messages: (1) within each paradigm, adding fat tails and leverage helps — GJR-t beats the plain GARCHs, SV-t-leverage beats plain SV; (2) matched feature-for-feature, the two paradigms are very close — GJR-t leads SV-$t$-leverage by only 0.004 nats, and plain Gaussian SV already edges out GARCH-$t$ (its scale-mixture predictive is fat-tailed for free). The next cell asks which translates into well-calibrated risk numbers.
What is a nat? The predictive score is the natural log of the predictive density, $\ln p(r_t\mid\text{past})$, so score differences are measured in nats — the base-$e$ unit of information (the natural-log analogue of a bit; $1$ nat $\approx 1.44$ bits). A gap of $\Delta$ nats in the mean daily score means one model assigns, on average, $e^{\Delta}$ times more predictive density to what actually happened. Here $e^{0.004}\approx 1.004$ — about $0.4\%$ more density per day, a near-tie; summed over the $\sim$2000-day test set the cumulative gap is only $\approx 7$ nats (the small vertical distance between the gold and purple curves below).
rt = ret[test_idx]
print("VaR backtest (test set 2001-2009). Kupiec/Christoffersen p<0.05 = REJECT the model.")
for al in alphas:
print("\n alpha = %.0f%% viol%% n Kupiec Christ. ES model/real verdict" % (al*100))
for nm, (_, var) in models.items():
v = var[al][test_idx]
bt = F.backtest(rt, v, al)
es_model = np.nanmean(var[al][test_idx][rt < v]) if (rt < v).any() else np.nan
ok = "PASS" if (bt["p_kupiec"] > .05 and bt["p_cc"] > .05) else "REJECT"
print(" %-9s %6.2f %4d p=%.3f p=%.3f %5.2f / %5.2f %s"
% (nm, bt["rate"]*100, bt["n_viol"], bt["p_kupiec"], bt["p_cc"], es_model, bt["es"], ok))
VaR backtest (test set 2001-2009). Kupiec/Christoffersen p<0.05 = REJECT the model.
alpha = 1% viol% n Kupiec Christ. ES model/real verdict
GARCH-N 1.58 32 p=0.015 p=0.044 -2.47 / -3.07 REJECT
GARCH-t 0.94 19 p=0.781 p=0.803 -2.68 / -3.35 PASS
GJR-t 0.74 15 p=0.221 p=0.422 -2.73 / -3.48 PASS
SV 1.29 26 p=0.217 p=0.298 -2.97 / -3.67 PASS
SV-t-lev 0.69 14 p=0.141 p=0.306 -2.70 / -3.43 PASS
alpha = 5% viol% n Kupiec Christ. ES model/real verdict
GARCH-N 5.39 109 p=0.429 p=0.504 -1.90 / -2.54 PASS
GARCH-t 6.18 125 p=0.019 p=0.031 -1.82 / -2.50 REJECT
GJR-t 6.03 122 p=0.039 p=0.118 -1.79 / -2.40 REJECT
SV 6.08 123 p=0.031 p=0.095 -1.89 / -2.58 REJECT
SV-t-lev 5.24 106 p=0.623 p=0.858 -1.98 / -2.63 PASS
d = pd.to_datetime(dates.values[test_idx])
fig, ax = plt.subplots(1, 2, figsize=(14, 4.7), gridspec_kw=dict(width_ratios=[3, 2]))
base = lps_gn[test_idx]
for nm, col in [("GARCH-t","seagreen"), ("GJR-t","goldenrod"), ("SV","firebrick"), ("SV-t-lev","darkviolet")]:
ax[0].plot(d, np.cumsum(models[nm][0][test_idx] - base), color=col, lw=1.7, label=nm)
ax[0].axhline(0, color="steelblue", lw=1.2, ls="--", label="GARCH-N (baseline)")
ax[0].axvspan(pd.Timestamp("2008-09-01"), pd.Timestamp("2009-01-30"), color="0.9", zorder=0)
ax[0].set_title("Cumulative 1-step log-predictive-score advantage over GARCH-N\n(rising = better forecast; shaded = 2008 crisis)")
ax[0].set_ylabel("cumulative LPS gain vs GARCH-N"); ax[0].legend(fontsize=8, loc="upper left")
w = (d >= pd.Timestamp("2008-01-01")) & (d <= pd.Timestamp("2009-01-30")); dv = d[w]
ax[1].plot(dv, np.abs(rt[w]), color="0.8", lw=.5, label="|return|")
ax[1].fill_between(dv, stlvol["lo"][test_idx][w], stlvol["hi"][test_idx][w], color="darkviolet", alpha=.2, label="SV-t-lev 90% band")
ax[1].plot(dv, stlvol["vol"][test_idx][w], color="darkviolet", lw=1.5, label="SV-t-lev (particle filter)")
ax[1].plot(dv, sd_jt[test_idx][w], color="goldenrod", lw=1.3, label="GJR-t")
ax[1].set_title("2008 crisis: the two full models head-to-head\n(SV-t-leverage carries a credible band; GJR-t is a single line)")
ax[1].set_ylabel("daily vol (%)"); ax[1].legend(fontsize=7.5, loc="upper left")
for lab in ax[1].get_xticklabels(): lab.set_rotation(20)
plt.tight_layout(); plt.show()
Results¶
Matched feature-for-feature, the two paradigms are a near dead heat — with a twist in the risk numbers.
- Density forecasting: GJR-t by a whisker, SV-$t$-leverage a very close second. GJR-t has the best log predictive score, but SV-$t$-leverage trails by only 0.004 nats and beats every other model. Both "full" models pull away from the plain ones during 2007–08 (the vertical climb in the left panel), and in the crisis window their predictive volatilities are almost indistinguishable (right panel) — except SV-$t$-leverage carries a credible band a deterministic GARCH cannot.
- Risk calibration: SV-$t$-leverage is the only model clean at both levels. At 1% every fat-tailed / mixture model passes and Gaussian GARCH-N is rejected (thin tail under-covers). At 5% the ranking flips: GARCH-$t$, GJR-$t$ and plain SV all take too many hits (~6% violations vs the 5% target, rejected) because a low-$\nu$ Student-$t$ is peaked in the middle and puts its 5% quantile too close to zero — while SV-$t$-leverage passes both 1% and 5% (0.69% and 5.24%, Kupiec $p=0.14$ and $0.62$). Its latent-vol dynamics plus a milder $\nu$ get the whole left tail right, not just the extreme.
- Plain SV punches above its weight. Even Gaussian, symmetric SV out-scores GARCH-$t$ and passes the 1% VaR — its particle-filter predictive is a scale-mixture (fat tails without an explicit $\nu$) and it carries volatility uncertainty. In the calm mid-2000s GARCH-$t$ actually forecasts worse than Gaussian GARCH-N (its fat tails cost it when nothing happens); SV stays above the baseline throughout.
Bottom line. With the same two features, latent-state SV and observation-driven GARCH forecast the return density essentially equally well; the practical edge goes to SV-$t$-leverage on tail calibration (the only model well-behaved at both VaR levels) and on honest uncertainty (a predictive band), at the cost of a particle filter instead of a closed form. The flexibility of a stochastic volatility is not free — but here it pays for itself where it matters most, in the tail.
References¶
- Kupiec, P. (1995). Techniques for verifying the accuracy of risk measurement models. J. Derivatives.
- Christoffersen, P. (1998). Evaluating interval forecasts. International Economic Review.
- Gordon, N., Salmond, D. & Smith, A. (1993). Novel approach to nonlinear/non-Gaussian Bayesian state estimation. IEE Proc-F — the bootstrap particle filter.
- Pitt, M. & Shephard, N. (1999). Filtering via simulation: auxiliary particle filters. JASA.
- Gneiting, T. & Raftery, A. (2007). Strictly proper scoring rules, prediction, and estimation. JASA — the log predictive score.
- Omori, Y., Chib, S., Shephard, N. & Nakajima, J. (2007). Stochastic volatility with leverage. J. Econometrics.