Value at Risk — models and backtests¶
An application of the volatility models: RiskMetrics · GARCH-t · GJR-t · Historical Simulation¶
The 1-day Value at Risk at level $\alpha$ is the quantile $\mathrm{VaR}_t^\alpha$ with $\Pr(r_t < \mathrm{VaR}_t^\alpha\mid\mathcal F_{t-1})=\alpha$. For a conditionally-scaled model $r_t=\sigma_t\varepsilon_t$ this is $\mathrm{VaR}_t^\alpha=\sigma_t\,q_\alpha$, where $q_\alpha$ is the $\alpha$-quantile of the (standardized) innovation. We compute it four ways and backtest each with the seminal tests:
- Kupiec (1995) — unconditional coverage $LR_{uc}$: is the violation rate equal to $\alpha$? $\ \sim\chi^2_1$.
- Christoffersen (1998) — independence $LR_{ind}$ and conditional coverage $LR_{cc}=LR_{uc}+LR_{ind}$: are violations at the right rate and not clustered? $\ \sim\chi^2_1,\chi^2_2$.
- Expected Shortfall — the average loss beyond VaR (the coherent tail measure).
Out-of-sample design: estimate on the first 3500 days (1987–2001), then filter the conditional volatility forward and backtest on the held-out 2001–2009 window (which contains the 2008 crisis).
The four VaR models¶
Each model produces a one-day-ahead conditional quantile $\mathrm{VaR}_t^\alpha$ (a negative number, in %). A realized return worse than it is a violation, recorded by the hit indicator $I_t=\mathbf 1[\,r_t<\mathrm{VaR}_t^\alpha\,]$. The four differ in how they model volatility and the tail:
1. Historical Simulation (HS). The empirical $\alpha$-quantile of the training returns — one constant number reused every day. It assumes nothing about the distribution, but it is unconditional: it ignores volatility clustering, so it reacts slowly and lets violations pile up in turbulent stretches.
2. RiskMetrics (J.P. Morgan, 1996). The industry benchmark. Volatility is an exponentially-weighted moving average (EWMA) of squared returns, $$\sigma_t^2=\lambda\,\sigma_{t-1}^2+(1-\lambda)\,r_{t-1}^2,\qquad \lambda=0.94\ \text{(the standard daily decay)},$$ and returns are taken conditionally Normal, so $\mathrm{VaR}_t^\alpha=\sigma_t\,\Phi^{-1}(\alpha)$ (e.g. $\Phi^{-1}(0.01)=-2.33$). Equivalently it is an IGARCH(1,1) with $\omega=0,\ \alpha=1-\lambda,\ \beta=\lambda$ — a GARCH whose persistence is pinned at exactly 1 (shocks to variance never decay) and with no parameters to estimate ($\lambda$ is fixed by convention). It tracks volatility well, but its Gaussian tail understates extreme losses at the 1% level.
3. GARCH(1,1)-t. The workhorse of the arc: $\sigma_t^2=\omega+\alpha r_{t-1}^2+\beta\sigma_{t-1}^2$ with Student-t innovations. Volatility mean-reverts (estimated persistence $<1$), and the fat-tailed quantile $q_\alpha=\sqrt{(\nu-2)/\nu}\;t_\nu^{-1}(\alpha)$ sizes the tail correctly.
4. GJR-GARCH-t. Adds the leverage term, so bad news raises the volatility forecast faster — the VaR widens more sharply going into a sell-off.
For the three conditional models, $\mathrm{VaR}_t^\alpha=\sigma_t\,q_\alpha$. Parameters are estimated on the training window only; the volatility is then filtered forward through the test window using past returns (so there is no look-ahead), and the backtests run on the test window.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy.stats import norm, t as tdist, chi2
import garch_scratch as G
d = pd.read_csv("sp500ret.csv", parse_dates=["date"]); dates = d["date"].values
ntr = 3500
r = d["ret"].values - d["ret"].values[:ntr].mean() # demean by the training mean (no look-ahead)
r2, _ = G.prepare(r); s0 = float(np.var(r[:ntr])); neg = (r < 0).astype(float)
te = slice(ntr, len(r)); rte = r[te]
print("train %d (%s..%s) test %d (%s..%s)" %
(ntr, str(dates[0])[:7], str(dates[ntr-1])[:7], len(rte), str(dates[ntr])[:7], str(dates[-1])[:7]))
# --- fit the conditional-vol models on the training window, then FILTER over the full series ---
mGt, *_ = G.garch_mle(r[:ntr], student=True) # GARCH-t theta=[o,a,b,nu]
mJt, *_ = G.garch_mle(r[:ntr], gjr=True, student=True) # GJR-t theta=[o,a,b,g,nu]
volGt = np.sqrt(G.garch_var(mGt, r2, s0))
volJt = np.sqrt(G.garch_var(mJt, r2, s0, neg=neg, gjr=True))
lam = 0.94; s2 = np.empty(len(r)); s2[0] = s0 # RiskMetrics EWMA
for t in range(1, len(r)): s2[t] = lam * s2[t - 1] + (1 - lam) * r[t - 1] ** 2
volEW = np.sqrt(s2)
def stdt_q(a, nu): return np.sqrt((nu - 2) / nu) * tdist.ppf(a, nu)
def make_var(alpha):
return {
"Historical Sim": np.full(len(r), np.quantile(r[:ntr], alpha)), # unconditional, constant
"RiskMetrics-N": volEW * norm.ppf(alpha),
"GARCH-t": volGt * stdt_q(alpha, mGt[3]),
"GJR-t": volJt * stdt_q(alpha, mJt[4]),
}
print("GARCH-t nu %.1f GJR-t: gamma %.3f nu %.1f" % (mGt[3], mJt[3], mJt[4]))
train 3500 (1987-03..2001-01) test 2023 (2001-01..2009-01) GARCH-t nu 5.2 GJR-t: gamma 0.084 nu 5.5
Backtesting: are the VaR forecasts calibrated?¶
Take the hit sequence $\{I_t\}$ on the test window. A correct $\alpha$-VaR needs two properties (Christoffersen, 1998):
Unconditional coverage — Kupiec (1995). The violation rate should equal $\alpha$. With $x$ violations in $n$ days ($\hat p=x/n$), the likelihood-ratio (proportion-of-failures) test $$LR_{uc}=-2\ln\frac{(1-\alpha)^{\,n-x}\,\alpha^{\,x}}{(1-\hat p)^{\,n-x}\,\hat p^{\,x}}\ \sim\ \chi^2_1$$ rejects when there are too many (model under-forecasts risk) or too few (too conservative) violations.
Independence — Christoffersen (1998). Violations must not cluster: a hit today should not raise the chance of a hit tomorrow. From the $2\times2$ transition counts $n_{ij}$ of $\{I_t\}$, an LR test pits a first-order Markov chain (transition probabilities $\pi_{01},\pi_{11}$) against the i.i.d. null, giving $LR_{ind}\sim\chi^2_1$. A model that fails to react to volatility bunches its violations in crises and is rejected here even when the average rate looks fine — this is the failure Kupiec alone cannot see.
Conditional coverage combines them: $LR_{cc}=LR_{uc}+LR_{ind}\sim\chi^2_2$. Below we report the $p$-values; $p>0.05$ = "not rejected" (the model passes the test).
Expected Shortfall (ES). VaR says how often you breach; ES says how bad the breach is — the average loss given a violation, $\mathrm{ES}_t^\alpha=\mathbb E[\,r_t\mid r_t<\mathrm{VaR}_t^\alpha\,]$. It is the coherent risk measure (Artzner et al., 1999) and the current Basel market-risk standard; the ES/VaR column below is the realized tail-loss relative to the VaR level, a calibration check on tail size.
# --- backtests: Kupiec (1995) uncond. coverage + Christoffersen (1998) independence / cond. coverage ---
def _ll(k, N):
if N == 0: return 0.0
ph = k / N; s = 0.0
if k > 0: s += k * np.log(ph)
if N - k > 0: s += (N - k) * np.log(1 - ph)
return s
def kupiec(h, p):
n = len(h); x = int(h.sum())
LR = -2 * ((x * np.log(p) + (n - x) * np.log(1 - p)) - _ll(x, n))
return LR, 1 - chi2.cdf(LR, 1)
def christoffersen(h, p):
hi = h.astype(int); n00 = n01 = n10 = n11 = 0
for i in range(1, len(hi)):
a, b = hi[i - 1], hi[i]
n00 += (a == 0) & (b == 0); n01 += (a == 0) & (b == 1)
n10 += (a == 1) & (b == 0); n11 += (a == 1) & (b == 1)
ll_alt = _ll(n01, n00 + n01) + _ll(n11, n10 + n11)
ll_null = _ll(n01 + n11, n00 + n01 + n10 + n11)
LR_ind = -2 * (ll_null - ll_alt)
LR_uc, _ = kupiec(h, p); LR_cc = LR_uc + LR_ind
return (1 - chi2.cdf(LR_ind, 1)), LR_cc, (1 - chi2.cdf(LR_cc, 2))
for alpha in (0.01, 0.05):
V = make_var(alpha)
print("\n=== %d%% VaR backtest (test window, n=%d) ===" % (int(alpha*100), len(rte)))
print("%-16s %7s %7s %9s %9s %8s" % ("model", "viol", "rate%", "Kupiec_p", "Christ_p", "ES/VaR"))
for name, var in V.items():
vt = var[te]; hit = rte < vt
_, pk = kupiec(hit, alpha); _, _, pcc = christoffersen(hit, alpha)
es = rte[hit].mean() if hit.any() else np.nan
es_ratio = es / vt[hit].mean() if hit.any() else np.nan
print("%-16s %7d %7.2f %9.3f %9.3f %8.2f" %
(name, hit.sum(), 100*hit.mean(), pk, pcc, es_ratio))
print(" (target rate %.1f%%; p>0.05 => not rejected)" % (100*alpha))
=== 1% VaR backtest (test window, n=2023) === model viol rate% Kupiec_p Christ_p ES/VaR Historical Sim 56 2.77 0.000 0.000 1.55 RiskMetrics-N 44 2.17 0.000 0.000 1.23 GARCH-t 19 0.94 0.781 0.803 1.26 GJR-t 16 0.79 0.327 0.544 1.27 (target rate 1.0%; p>0.05 => not rejected) === 5% VaR backtest (test window, n=2023) === model viol rate% Kupiec_p Christ_p ES/VaR Historical Sim 179 8.85 0.000 0.000 1.69 RiskMetrics-N 126 6.23 0.014 0.046 1.36 GARCH-t 131 6.48 0.004 0.010 1.37 GJR-t 126 6.23 0.014 0.048 1.34 (target rate 5.0%; p>0.05 => not rejected)
# 1% VaR over the crisis-containing test window: adaptive (GARCH-t) vs static (Historical Sim)
V = make_var(0.01); dt = pd.to_datetime(dates[te])
fig, ax = plt.subplots(2, 1, figsize=(12, 7), sharex=True)
ax[0].plot(dt, rte, color="0.75", lw=.4, label="return")
for name, col in [("GARCH-t", "firebrick"), ("RiskMetrics-N", "steelblue"), ("Historical Sim", "seagreen")]:
ax[0].plot(dt, V[name][te], color=col, lw=1.0, label="%s 1%% VaR" % name)
ax[0].set_ylabel("return / VaR (%)"); ax[0].legend(fontsize=8, ncol=2); ax[0].set_ylim(-12, 6)
ax[0].set_title("1% VaR through 2001-2009: model VaR adapts to volatility; Historical Sim is static")
# violations: GARCH-t (adaptive, spread out) vs Historical Sim (clustered in 2008)
for k, (name, col) in enumerate([("GARCH-t", "firebrick"), ("Historical Sim", "seagreen")]):
hit = rte < V[name][te]
ax[1].plot(dt[hit], np.full(hit.sum(), 1 - k), "|", color=col, ms=12,
label="%s violations (n=%d)" % (name, hit.sum()))
ax[1].set_yticks([0, 1]); ax[1].set_yticklabels(["Historical Sim", "GARCH-t"]); ax[1].set_ylim(-0.5, 1.5)
ax[1].set_title("Where the 1% violations fall — static VaR clusters them in the crisis (fails Christoffersen)")
ax[1].legend(fontsize=8, loc="upper left")
plt.tight_layout(); plt.savefig("var_backtest.png", dpi=120, bbox_inches="tight"); plt.show()
Value at Risk in dollars — and the size of the breaches¶
VaR is defined as a return quantile, but risk is managed in money. For a $1,000,000 position in the index, the daily P&L is $\text{notional}\times r_t/100$ and the dollar VaR is $\text{notional}\times\mathrm{VaR}_t^\alpha/100$. The top panel shows the realized P&L against the 1% dollar VaR (violations = a daily loss deeper than the VaR line); the bottom panel shows, for each violation day, how many dollars the loss ran past the VaR — the realized dollar Expected Shortfall, i.e. the losses a VaR limit does not cap.
import matplotlib.ticker as mtick
NOTIONAL = 1_000_000
V1 = make_var(0.01); dt = pd.to_datetime(dates[te])
pnl = NOTIONAL * rte / 100.0
varG = NOTIONAL * V1["GARCH-t"][te] / 100.0
hit = rte < V1["GARCH-t"][te]
usd = mtick.StrMethodFormatter("\\${x:,.0f}") # literal $ (escaped for matplotlib mathtext)
fig, ax = plt.subplots(2, 1, figsize=(12, 7.4), sharex=True)
ax[0].fill_between(dt, 0, pnl, where=(pnl < 0), color="lightcoral", alpha=.55, step="mid", label="daily loss")
ax[0].fill_between(dt, 0, pnl, where=(pnl >= 0), color="0.9", step="mid")
ax[0].plot(dt, varG, color="firebrick", lw=1.1, label="GARCH-t 1% VaR")
ax[0].scatter(dt[hit], pnl[hit], color="darkred", s=24, zorder=5, label="violation (%d days)" % hit.sum())
ax[0].axhline(0, color="k", lw=.5); ax[0].yaxis.set_major_formatter(usd)
ax[0].set_ylabel("daily profit / loss (P&L)"); ax[0].legend(fontsize=8, loc="lower left")
wi = int(np.argmin(pnl))
ax[0].set_title("\\$%s S&P position: daily P&L vs the 1%% VaR (worst day \\$%s on %s)"
% (f"{NOTIONAL:,}", f"{pnl[wi]:,.0f}", str(dt[wi].date())))
exc = pnl[hit] - varG[hit] # dollars lost beyond the VaR (negative)
ax[1].bar(dt[hit], exc, width=7, color="firebrick")
ax[1].axhline(0, color="k", lw=.5); ax[1].yaxis.set_major_formatter(usd)
ax[1].set_ylabel("loss beyond VaR")
ax[1].set_title("Breach severity - dollars lost PAST the VaR on violation days (the realized \\$ Expected Shortfall)")
ax[1].text(0.005, 0.06, "on breach days: total \\$%s beyond VaR, mean \\$%s, worst \\$%s"
% (f"{exc.sum():,.0f}", f"{exc.mean():,.0f}", f"{exc.min():,.0f}"),
transform=ax[1].transAxes, fontsize=8)
plt.tight_layout(); plt.savefig("var_dollar.png", dpi=120, bbox_inches="tight"); plt.show()
Cumulative exceptions — calibration over time¶
The single most diagnostic VaR picture: each model's running count of 1% violations against the expected line $0.01\times(\text{days elapsed})$. A well-calibrated model tracks the diagonal; a model that under-forecasts risk drifts above it, and one that clusters its violations does so in steps — both pathologies visible at a glance, and both concentrated in 2008.
fig, ax = plt.subplots(figsize=(11.5, 4.7))
n = len(rte); dt = pd.to_datetime(dates[te]); V1 = make_var(0.01)
ax.plot(dt, 0.01 * np.arange(1, n + 1), "k--", lw=1.3, label="expected (1% × days)")
for name, col in [("Historical Sim", "seagreen"), ("RiskMetrics-N", "steelblue"),
("GARCH-t", "firebrick"), ("GJR-t", "purple")]:
cum = np.cumsum(rte < V1[name][te])
ax.plot(dt, cum, color=col, lw=1.6, label="%s (%d total)" % (name, cum[-1]))
ax.set_ylabel("cumulative 1% violations"); ax.legend(fontsize=8, loc="upper left")
ax.set_title("Cumulative VaR exceptions vs expected — static models (HS, RiskMetrics) run above the line and jump in 2008")
plt.tight_layout(); plt.savefig("var_cumexc.png", dpi=120, bbox_inches="tight"); plt.show()
Results¶
- Historical Simulation fails outright. Its constant, unconditional VaR ignores volatility clustering: it takes far too many hits overall (2.8% vs the 1% target — rejected by Kupiec) and they bunch in 2008 (rejected by Christoffersen conditional coverage). The clustering is exactly the failure mode conditional-coverage tests were built to catch, and the seminal argument for conditional risk models.
- RiskMetrics (EWMA) adapts to volatility and does far better on clustering, but its Gaussian tail understates extreme losses at 1% (too many violations in turbulent stretches).
- GARCH-t / GJR-t pass both Kupiec and Christoffersen cleanly at 1% (rates 0.94% / 0.79%, all p > 0.3): the conditional volatility spreads violations out over time and the Student-t tail sizes the 1% quantile correctly — you need both. At 5%, though, even these models reject (GARCH-t Kupiec p = 0.004, GJR-t 0.014; rates ~6.2–6.5% vs the 5% target): the 5% quantile is a shoulder rather than a deep tail, so the fat-tail advantage fades and small miscalibration in the fitted volatility level shows through. The 1% level is where the models cleanly separate.
- Expected Shortfall (the tail-average loss) is reported alongside; the model VaRs whose ES/VaR ratio is stable are the coherent, well-calibrated ones.
This is the risk-management application of the volatility arc — and the VaR half of the planned sv_vs_garch head-to-head, to which SV-t VaR will be added. Next in the applications phase: option pricing (Heston–Nandi closed-form GARCH, Duan Monte-Carlo, Heston stochastic-vol).