Repository Public Research preview
1"""One-day-ahead VaR / ES forecasts and backtests (percent returns).23VaR_p(t+1) = mu_{t+1} + sqrt(h_{t+1}) q_p ; ES_p(t+1) = mu_{t+1} + sqrt(h_{t+1}) e_p, where q_p / e_p are the p-quantile4and lower-tail mean of the model's standardized innovation law (normal, t, ged, skewt via wp46.models.distributions). For5distribution-free models ('empirical'), q_p and e_p are the empirical quantile / tail mean of the standardized in-window6residuals z = (r - mu)/sqrt(sigma2_filtered) supplied by the caller (the driver stores them per refit as `emp_q`).7Backtests: Kupiec (1995) unconditional coverage LR, Christoffersen (1998) conditional coverage LR, Engle & Manganelli8(2004) dynamic quantile test with 4 lagged hits and the VaR as regressors, Fissler & Ziegel (2016) FZ0 joint loss9(Patton, Ziegel & Chen 2019 parameterization): FZ0 = -1/(alpha ES) 1{r<=VaR} (VaR - r) + VaR/ES + log(-ES) - 1, valid for ES < 0."""10from __future__ import annotations1112import numpy as np13import pandas as pd14from scipy import stats1516from wp46.models.distributions import expected_shortfall, ppf1718LEVELS = (0.01, 0.05)192021def var_es_forecasts(mu: np.ndarray, h1: np.ndarray, dist: str, dp1: np.ndarray, dp2: np.ndarray, levels=LEVELS,22 emp_q: dict | None = None) -> dict:23 """Vectorized over origins; distribution parameters may vary by origin (refits)."""24 out = {}25 s = np.sqrt(np.maximum(h1, 1e-12))26 for p in levels:27 q = np.full(len(mu), np.nan); e = np.full(len(mu), np.nan)28 if dist == "empirical":29 if emp_q is not None:30 q[:] = emp_q[p]["q"]; e[:] = emp_q[p]["es"]31 elif dist == "normal":32 q[:] = stats.norm.ppf(p); e[:] = -stats.norm.pdf(stats.norm.ppf(p)) / p33 else:34 dp1c = np.where(np.isfinite(dp1), np.maximum(dp1, 2.051 if dist in ("t", "skewt") else 0.31), 5.0) # bound-safe shape35 params = np.column_stack([dp1c, dp2]) if dist == "skewt" else dp1c[:, None]36 uniq, inv = np.unique(np.round(params, 10), axis=0, return_inverse=True)37 for k, prm in enumerate(uniq):38 m = inv.ravel() == k39 q[m] = ppf(dist, p, prm)[0]40 e[m] = expected_shortfall(dist, p, prm)41 out[p] = {"var": mu + s * q, "es": mu + s * e}42 return out434445def kupiec(hits: np.ndarray, p: float) -> dict:46 n = len(hits); x = int(hits.sum())47 if n == 0:48 return dict(lr=np.nan, p=np.nan, n=0, x=0, rate=np.nan)49 pi = x / n50 with np.errstate(divide="ignore", invalid="ignore"):51 l0 = x * np.log(p) + (n - x) * np.log(1 - p)52 l1 = (x * np.log(pi) if x > 0 else 0.0) + ((n - x) * np.log(1 - pi) if x < n else 0.0)53 lr = -2 * (l0 - l1)54 return dict(lr=float(lr), p=float(stats.chi2.sf(lr, 1)), n=n, x=x, rate=pi)555657def christoffersen_cc(hits: np.ndarray, p: float) -> dict:58 h = hits.astype(int)59 if len(h) < 2:60 return dict(lr_ind=np.nan, p_ind=np.nan, lr_cc=np.nan, p_cc=np.nan)61 n00 = np.sum((h[:-1] == 0) & (h[1:] == 0)); n01 = np.sum((h[:-1] == 0) & (h[1:] == 1))62 n10 = np.sum((h[:-1] == 1) & (h[1:] == 0)); n11 = np.sum((h[:-1] == 1) & (h[1:] == 1))63 p01 = n01 / max(n00 + n01, 1); p11 = n11 / max(n10 + n11, 1); pi = (n01 + n11) / max(len(h) - 1, 1)6465 def ll(a, b, c, d, q0, q1):66 with np.errstate(divide="ignore", invalid="ignore"):67 t = 0.068 for cnt, pr in ((a, 1 - q0), (b, q0), (c, 1 - q1), (d, q1)):69 if cnt > 0:70 t += cnt * np.log(pr)71 return t72 l_ind = ll(n00, n01, n10, n11, p01, p11); l_0 = ll(n00, n01, n10, n11, pi, pi)73 lr_ind = -2 * (l_0 - l_ind)74 lr_uc = kupiec(hits, p)["lr"]75 lr_cc = lr_uc + lr_ind76 return dict(lr_ind=float(lr_ind), p_ind=float(stats.chi2.sf(lr_ind, 1)), lr_cc=float(lr_cc), p_cc=float(stats.chi2.sf(lr_cc, 2)))777879def engle_manganelli_dq(hits: np.ndarray, var: np.ndarray, p: float, lags: int = 4) -> dict:80 h = hits.astype(float) - p81 n = len(h)82 if n <= lags + 10:83 return dict(stat=np.nan, p=np.nan)84 X = np.column_stack([np.ones(n - lags)] + [h[lags - j - 1:n - j - 1] for j in range(lags)] + [var[lags:]])85 y = h[lags:]86 XtX = X.T @ X87 try:88 b = np.linalg.solve(XtX, X.T @ y)89 except np.linalg.LinAlgError:90 return dict(stat=np.nan, p=np.nan)91 stat = float(b @ XtX @ b / (p * (1 - p)))92 return dict(stat=stat, p=float(stats.chi2.sf(stat, X.shape[1])))939495def fz0_loss(r: np.ndarray, var: np.ndarray, es: np.ndarray, p: float) -> np.ndarray:96 """FZ0 is undefined for ES >= 0 and numerically meaningless for |ES| below 1e-3 percent (degenerate variance forecasts):97 such origins are returned as NaN and excluded from the mean (their count is reported as `fz0_invalid`)."""98 es = np.where(es < -1e-3, es, np.nan)99 hit = (r <= var).astype(float)100 with np.errstate(invalid="ignore", divide="ignore"):101 return -hit * (var - r) / (p * es) + var / es + np.log(-es) - 1.0102103104def backtest_table(r: np.ndarray, fc: dict, p: float) -> dict:105 var, es = fc[p]["var"], fc[p]["es"]106 m = np.isfinite(r) & np.isfinite(var) & np.isfinite(es)107 r, var, es = r[m], var[m], es[m]108 hits = (r <= var)109 k = kupiec(hits, p); c = christoffersen_cc(hits, p); d = engle_manganelli_dq(hits, var, p)110 fz = fz0_loss(r, var, es, p)111 return dict(n=int(len(r)), hit_rate=k["rate"], kupiec_p=k["p"], cc_p=c["p_cc"], dq_p=d["p"], fz0_mean=float(np.nanmean(fz)),112 fz0_invalid=int(np.sum(~np.isfinite(fz))),113 es_shortfall_ratio=float(np.mean(r[hits]) / np.mean(es[hits])) if hits.sum() > 0 else np.nan)114