--- title: "5.7.3 Probabilistic models: quantiles and DeepAR" description: "Quantile loss, a monotone quantile network, DeepAR from scratch and in neuralforecast, GJR-GARCH quantiles and conditional calibration of intervals." url: "https://pythonparatrading.com/en/python-trading-course/level-5/probabilistic-forecasting-quantile-loss-deepar/" lang: "en" image: "https://pythonparatrading.com/og/en/python-trading-course/level-5/probabilistic-forecasting-quantile-loss-deepar.png" --- # 5.7.3 · Probabilistic models: quantile loss, DeepAR and intervals > To size a position, set a stop or compute a VaR, the mean forecast is not enough: you need the distribution. This lesson predicts quantiles of tomorrow's return with a network trained on the quantile loss and with a DeepAR, compares them with the quantiles of a GJR-GARCH and of an oracle, and shows how to check that a 90% interval covers 90% on the days that matter too. - Course: [Level 5: Machine learning and deep learning (Keras/PyTorch)](https://pythonparatrading.com/en/python-trading-course/level-5/) - Module: 5.7 · Deep learning with PyTorch - Duration: 75 min - Prerequisites: [5.7.1 · PyTorch: tensors, autograd and the training loop](https://pythonparatrading.com/en/python-trading-course/level-5/pytorch-tensors-autograd-training-loop/), [4.2.3 · VaR, Expected Shortfall and VaR backtesting](https://pythonparatrading.com/en/python-trading-course/level-4/var-expected-shortfall-and-var-backtesting/), [2.1.4 · Volatility models: ARCH, GARCH, EGARCH and GJR with the arch library](https://pythonparatrading.com/en/python-trading-course/level-2/volatility-models-arch-garch-egarch-gjr/) ## By the end of this lesson - Explain the quantile (pinball) loss and check that its minimum is the empirical quantile. - Train a network that forecasts several quantiles at once without them crossing, and a DeepAR that forecasts a Student's t distribution per step. - Build the benchmarks: unconditional quantiles, HAR scale, GJR-GARCH with `arch` and an oracle with the true variance. - Measure calibration: total coverage, coverage by volatility bucket, exceedances with the Kupiec test and mean pinball loss. - Adjust an interval with conformal prediction on the validation set and understand why, in time series, it does not guarantee coverage. ## From a number to a distribution Lessons 5.7.1 and 5.7.2 predicted a number: tomorrow's log RV. Many trading decisions need more. A VaR is a quantile (lesson 4.2.3); a sensible stop depends on how far the price can move on a normal day; the size of a position under volatility targeting, on the expected dispersion. A **probabilistic model** forecasts the distribution of the future value, or at least several of its quantiles, and prediction intervals come from there. Two intervals that are often confused: | | Confidence interval | Prediction interval | |---|---|---| | About what | a parameter (the mean, a slope) | the next observation | | What it captures | the uncertainty of the estimate | that uncertainty plus the noise of the data point itself | | With more data | it narrows toward zero | it narrows down to the conditional dispersion, no further | | In this lesson | — | the 90% interval of tomorrow's return | The problem: forecast quantiles of tomorrow's **return** for the six assets of `vol5_diario.parquet`. The data README says it clearly: the conditional mean of the return is constant, so there is nothing to forecast in the mean. But the **dispersion** changes a lot and is predictable (it is a HAR on the logarithm of variance): a good probabilistic model has to widen the interval on volatile days and narrow it on calm ones. Synthetic data with a known truth: `vol5_verdad.parquet` gives the expected integrated variance of each day, the best possible scale data. ```python import copy import logging import warnings import numpy as np import pandas as pd import torch from torch import nn from scipy import stats torch.manual_seed(0) torch.use_deterministic_algorithms(True) TAU = np.array([0.01, 0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95, 0.99]) I05, I95 = 1, 7 # positions of 5% and 95% v = pd.read_parquet("vol5_diario.parquet").merge( pd.read_parquet("vol5_verdad.parquet")[ ["fecha", "activo", "iv_prevista"]], on=["fecha", "activo"]) def filas(g): x = np.log(g["rv_5min"]) return pd.DataFrame({ "fecha": g["fecha"], "activo": g["activo"], "r": g["rendimiento"], "d": x, "w": x.rolling(5).mean(), "m": x.rolling(22).mean(), "d_sig": x.shift(-1), "r_sig": g["rendimiento"].shift(-1), "iv_sig": g["iv_prevista"].shift(-1), "fecha_y": g["fecha"].shift(-1)}) h = (pd.concat([filas(g) for _, g in v.groupby("activo")]) .dropna().reset_index(drop=True)) h["tramo"] = np.select( [h["fecha_y"] < "2016-01-01", (h["fecha"] >= "2016-01-01") & (h["fecha_y"] < "2018-01-01"), h["fecha"] >= "2018-01-01"], ["ent", "val", "prueba"], "purga") ent, val, pru = (h["tramo"] == k for k in ("ent", "val", "prueba")) escala = h[ent].groupby("activo")["r"].std() # 2008-2015 only centro = h[ent].groupby("activo")["d"].mean() h["z"] = h["r"] / h["activo"].map(escala) h["z_sig"] = h["r_sig"] / h["activo"].map(escala) for c in ["d", "w", "m", "d_sig"]: h[c] = h[c] - h["activo"].map(centro) y = h.loc[pru, "z_sig"].to_numpy() print(h["tramo"].value_counts().to_string()) print(f"kurtosis of z_sig on the test set: {stats.kurtosis(y):.1f}") ``` ```output tramo prueba 12060 ent 11958 val 3012 purga 12 kurtosis of z_sig on the test set: 3.5 ``` Each return is divided by the standard deviation of its asset in training (`z`), so the six share a scale and can be modeled together, as in the two previous lessons (`ent`, `val` and `prueba` are, as there, training, validation and test). The tails are fat: excess kurtosis of 3.5 on the test set, as in real daily returns. The features are those of the HAR (today's log RV and its 5- and 22-day means, centered per asset) plus today's return. Nine quantiles are forecast at once, from 1% to 99%; the 90% interval goes from the 5% quantile to the 95% one. ## The quantile loss To train a model that forecasts the $\tau$ quantile, the **pinball loss** is used (Koenker and Bassett, 1978): with $u = y - q$ the error, $$ \rho_\tau(u) = \max\big(\tau u,\; (\tau - 1)\, u\big). $$ If the data point lies above the forecast quantile ($u > 0$), it costs $\tau u$; if it lies below, $(1 - \tau)|u|$. For $\tau = 0.05$, falling short on the downside costs 19 times more than overshooting on the upside, and the minimum is reached when only 5% of the data lie below: at the quantile. With $\tau = 0.5$ it is half the absolute error and gives the median. *Fig. — Pinball loss ρτ(u) = max(τu, (τ − 1)u) for the 5%, 50% and 95% quantiles: the 5% quantile penalises sitting above the data 19 times more than below it.* Three asymmetric V-shaped functions with the vertex at u = 0. For τ = 0.05 (highlighted), the left branch rises with slope 0.95 and the right with 0.05; for τ = 0.5 both branches have slope 0.5; for τ = 0.95 (dashed), it is the mirror image of τ = 0.05. The way to check it is the one from lesson 5.7.1: find with autograd the constant that minimizes the mean loss and compare it with the empirical quantile. ```python cripto = torch.tensor(h.loc[ent & (h["activo"] == "cripto"), "z_sig"].to_numpy()) def pinball(q, y, tau): u = y - q return torch.maximum(tau * u, (tau - 1) * u).mean() for tau in (0.05, 0.50, 0.95): c = torch.zeros((), dtype=torch.float64, requires_grad=True) opt = torch.optim.Adam([c], lr=0.05) paso = torch.optim.lr_scheduler.LambdaLR( opt, lambda k: 1 - k / 2000) # rate that decays to zero for _ in range(2000): opt.zero_grad() pinball(c, cripto, tau).backward() opt.step() paso.step() print(f"τ = {tau:.2f}: pinball minimum {c.item():+.3f}, " f"empirical quantile {np.quantile(cripto.numpy(), tau):+.3f}") ``` ```output τ = 0.05: pinball minimum -1.652, empirical quantile -1.650 τ = 0.50: pinball minimum +0.052, empirical quantile +0.052 τ = 0.95: pinball minimum +1.504, empirical quantile +1.504 ``` They agree; the two thousandths at 5% are `np.quantile`'s interpolation between two consecutive data points, a stretch where the pinball loss is almost flat. The learning rate that decays to zero is not decoration: the pinball gradient moves in jumps (each data point contributes $-\tau$ or $1 - \tau$), and with a fixed rate Adam oscillates around the minimum with an amplitude of the order of that rate. The mean pinball loss over several quantiles is also a **proper scoring rule** (Gneiting and Raftery, 2007): in expectation, no model scores better than the one that gives the true distribution. That is why it serves to compare probabilistic models, and coverage alone does not: an interval can cover 90% while being useless, as will be seen. ## The benchmarks: unconditional, HAR scale, oracle and GARCH Before any network, four ways of forecasting quantiles that set the scale: - **Unconditional**: the empirical quantiles of `z_sig` in training, the same every day. - **HAR scale**: tomorrow's volatility forecast with the RV HAR, $\hat\sigma_t = e^{\hat x_{t+1}/2}$, times the empirical quantiles of the returns divided by it. It separates the scale (predictable) from the shape (stable). - **Scale oracle**: the same with the true expected integrated variance, `iv_prevista`. Nobody has it; it marks the floor. - **GJR-GARCH(1,1) with Student's t** (lesson 2.1.4), estimated per asset with `arch` on 2008-2015 and applied without re-estimation: its quantiles are $\mu + \sigma_{t+1}\, t_\nu^{-1}(\tau)$, with the t standardized to unit variance. ```python def cuantiles_por_escala(sigma, en): """sigma (all rows) times the empirical quantiles of z_sig / sigma in training.""" forma = np.quantile(h.loc[ent, "z_sig"] / sigma[ent], TAU) return sigma[en].to_numpy()[:, None] * forma def har_x(en): return np.c_[np.ones(en.sum()), h.loc[en, ["d", "w", "m"]]] b = np.linalg.lstsq(har_x(ent), h.loc[ent, "d_sig"], rcond=None)[0] sig_har = pd.Series(np.exp(0.5 * har_x(h["tramo"] != "") @ b)) sig_ora = np.sqrt(h["iv_sig"]) / h["activo"].map(escala) Q = {"unconditional": np.tile(np.quantile(h.loc[ent, "z_sig"], TAU), (pru.sum(), 1)), "HAR scale": cuantiles_por_escala(sig_har, pru), "oracle": cuantiles_por_escala(sig_ora, pru)} print({k: q.shape for k, q in Q.items()}) ``` ```output {'unconditional': (12060, 9), 'HAR scale': (12060, 9), 'oracle': (12060, 9)} ``` ```python from arch import arch_model q_garch = np.zeros((pru.sum(), len(TAU))) param = {} for a in sorted(h["activo"].unique()): r = v[v["activo"] == a].set_index("fecha")["rendimiento"] * 100 modelo = arch_model(r, mean="Constant", vol="GARCH", p=1, o=1, q=1, dist="t") res = modelo.fit(last_obs="2016-01-01", disp="off") # forecasts made on each date t >= 2018 for t + 1 var = res.forecast(start="2018-01-01", horizon=1, reindex=False).variance["h.1"] t_tip = modelo.distribution.ppf(TAU, res.params[["nu"]]) filas_a = (h.loc[pru, "activo"] == a).to_numpy() sigma = np.sqrt(var.reindex(h.loc[pru, "fecha"][filas_a])) q_garch[filas_a] = ((res.params["mu"] + sigma.to_numpy()[:, None] * t_tip) / 100 / escala[a]) param[a] = res.params[["alpha[1]", "gamma[1]", "beta[1]", "nu"]] Q["GJR-GARCH-t"] = q_garch pd.DataFrame(param).T.round(3) ``` ```output alpha[1] gamma[1] beta[1] nu accion 0.000 0.045 0.964 3.431 bono 0.023 0.051 0.856 5.276 cripto 0.016 0.034 0.946 3.297 divisa 0.025 0.016 0.942 4.355 indice 0.015 0.021 0.962 4.806 oro 0.046 -0.046 0.886 4.907 ``` With `last_obs` the model is estimated on the data before 2016; `forecast(start=…)` applies those parameters while filtering all the data and returns, for each date $t$, the variance forecast for $t + 1$, which is exactly what is needed: the row of test date $t$ with its target for tomorrow. Returns are in percent (multiplied by 100) because `arch`'s optimizer converges better at that scale, and the `ppf` of its t distribution is already standardized to unit variance. The persistences ($\alpha + \gamma/2 + \beta$) are around 0.9-0.99 and the degrees of freedom between 3.3 and 5.3: fat tails in all of them. ## A quantile network The direct way to forecast quantiles with deep learning: a network with one output per quantile, trained with the sum of the pinball losses. It has a known problem, **quantile crossing**: nothing forces the 5% quantile to stay below the 10% one. The solution is to build it in: the network gives the first quantile and the **jumps** between consecutive quantiles, which go through `softplus` (always positive) and are accumulated. ```python FE = ["d", "w", "m", "z"] med, des = h.loc[ent, FE].mean(), h.loc[ent, FE].std() def tensores(en): X = torch.tensor(((h.loc[en, FE] - med) / des).to_numpy(), dtype=torch.float32) return X, torch.tensor(h.loc[en, "z_sig"].to_numpy(), dtype=torch.float32) datos = {k: tensores(en) for k, en in [("ent", ent), ("val", val), ("prueba", pru)]} tau_t = torch.tensor(TAU, dtype=torch.float32) def pinball_media(q, y): u = y[:, None] - q return torch.maximum(tau_t * u, (tau_t - 1) * u).mean() class RedCuantiles(nn.Module): def __init__(self, n, k=len(TAU), oculta=32): super().__init__() self.red = nn.Sequential( nn.Linear(n, oculta), nn.ReLU(), nn.Linear(oculta, oculta), nn.ReLU(), nn.Linear(oculta, k)) def forward(self, x): o = self.red(x) saltos = nn.functional.softplus(o[:, 1:]) return torch.cat([o[:, :1], o[:, :1] + torch.cumsum(saltos, dim=1)], 1) def entrenar(modelo, perdida, datos, epocas=30, lr=1e-3, paciencia=5, lote=256, perdida_val=None): perdida_val = perdida_val or perdida (Xe, Ye), (Xv, Yv) = datos["ent"], datos["val"] opt = torch.optim.Adam(modelo.parameters(), lr=lr) gen = torch.Generator().manual_seed(0) mejor, sin_mejora = (np.inf, None, -1), 0 for epoca in range(epocas): modelo.train() orden = torch.randperm(len(Xe), generator=gen) for i in range(0, len(Xe), lote): j = orden[i:i + lote] opt.zero_grad() perdida(modelo(Xe[j]), Ye[j]).backward() opt.step() modelo.eval() with torch.no_grad(): pv = perdida_val(modelo(Xv), Yv).item() if pv < mejor[0]: mejor = (pv, copy.deepcopy(modelo.state_dict()), epoca) sin_mejora = 0 else: sin_mejora += 1 if sin_mejora >= paciencia: break modelo.load_state_dict(mejor[1]) return mejor[2] torch.manual_seed(0) red_q = RedCuantiles(len(FE)) print("best epoch:", entrenar(red_q, pinball_media, datos)) red_q.eval() with torch.no_grad(): Q["quantile network"] = red_q(datos["prueba"][0]).double().numpy() q_val_red = red_q(datos["val"][0]).double().numpy() cruces = (np.diff(Q["quantile network"], axis=1) <= 0).sum() print("crossed quantiles on the test set:", cruces) ``` ```output best epoch: 27 crossed quantiles on the test set: 0 ``` Zero crossings, by construction. The price of monotonicity is nil: the network can represent any set of ordered quantiles. ## DeepAR: one distribution per step **DeepAR** (Salinas, Flunkert and Gasthaus, 2017; Amazon) is a recurrent network that, at each step of the series, receives the previous value and the covariates and returns the **parameters of a distribution** for the next value. It is trained by maximizing the likelihood of what was observed at every step of the window, and it is used in two ways: for one step, reading the quantiles of the forecast distribution; for several, sampling a value, feeding it in as the input of the next step and repeating (Monte Carlo paths). ```mermaid flowchart TB E[z today and log RV] --> R[LSTM] H[Previous state] --> R R --> P[μ, σ, ν] P --> D[Student's t] D --> C[Quantiles] D -. sampling .-> E ``` The differences from the quantile network: the shape of the distribution is chosen by the modeler (here a **Student's t**, which lets fat tails be learned through a parameter ν), and the extreme quantiles come from that shape, not from estimating each one separately with the few data points in the tail. The version of our own below uses 32-day windows, a 16-unit LSTM and log RV as a covariate; it is small so that it trains in seconds (in production, more units, more epochs and several seeds). ```python LV = 32 h["xs"] = h["d"] / h.loc[ent, "d"].std() act = h["activo"].to_numpy() zz, xx = h["z"].to_numpy(np.float32), h["xs"].to_numpy(np.float32) zs = h["z_sig"].to_numpy(np.float32) def secuencias(en): """Windows of LV days (z and log RV) ending at each row of `en`, and the target of each step (the next day's z).""" fin = [i for i in np.flatnonzero(en.to_numpy()) if i - LV + 1 >= 0 and act[i - LV + 1] == act[i]] S = np.stack([np.c_[zz[i - LV + 1:i + 1], xx[i - LV + 1:i + 1]] for i in fin]) Y = np.stack([zs[i - LV + 1:i + 1] for i in fin]) return torch.tensor(S), torch.tensor(Y) seq = {k: secuencias(en) for k, en in [("ent", ent), ("val", val), ("prueba", pru)]} assert len(seq["prueba"][0]) == pru.sum() class DeepAR(nn.Module): def __init__(self, entradas=2, oculta=16): super().__init__() self.lstm = nn.LSTM(entradas, oculta, batch_first=True) self.cabeza = nn.Linear(oculta, 3) def forward(self, s): # (n, LV, 2) o = self.cabeza(self.lstm(s)[0]) # (n, LV, 3) sigma = nn.functional.softplus(o[..., 1]) + 1e-3 nu = 2 + nn.functional.softplus(o[..., 2]) return o[..., 0], sigma, nu def nll(salida, y): # all steps mu, sigma, nu = salida t = torch.distributions.StudentT(nu, mu, sigma) return -t.log_prob(y).mean() def nll_ultimo(salida, y): # only the last one return nll([p[:, -1] for p in salida], y[:, -1]) torch.manual_seed(0) deepar = DeepAR() print("best epoch:", entrenar(deepar, nll, seq, epocas=20, paciencia=4, perdida_val=nll_ultimo)) deepar.eval() with torch.no_grad(): mu, sigma, nu = (p[:, -1].double().numpy() for p in deepar(seq["prueba"][0])) Q["DeepAR"] = (mu[:, None] + sigma[:, None] * stats.t.ppf(TAU[None, :], nu[:, None])) print(f"mean ν on the test set: {nu.mean():.1f}") ``` ```output best epoch: 17 mean ν on the test set: 4.6 ``` - Training uses the likelihood of the 32 steps of each window (that is how DeepAR trains); validation uses only that of the last one, which is the one used for forecasting and the only one whose target falls entirely in 2016-2017. - `torch.distributions.StudentT` gives the log-likelihood with a gradient; it has no quantile function, so the quantiles come from `scipy.stats.t.ppf` with the forecast parameters. With a normal, `Normal.icdf` would do it in torch. - The mean ν it learns is around 4-5, in line with the degrees of freedom the GARCH estimates. The library version is `neuralforecast.models.DeepAR`, with the same data format as lesson 5.7.2, the distribution in the loss (`DistributionLoss`) and log RV as a historical covariate (`hist_exog_list`): it is only known up to today. Its intervals come from sampling the forecast distribution (`DeepAR-lo-90`, `DeepAR-hi-90`). ```python from neuralforecast import NeuralForecast from neuralforecast.models import DeepAR as DeepARNF from neuralforecast.losses.pytorch import DistributionLoss, MQLoss for n in ("pytorch_lightning", "lightning.pytorch", "lightning_fabric", "lightning", "neuralforecast"): logging.getLogger(n).setLevel(logging.ERROR) warnings.filterwarnings("ignore", module="pytorch_lightning") warnings.filterwarnings("ignore", module="lightning") c_lrv = v[v["fecha"] < "2016-01-01"].groupby("activo")["rv_5min"] \ .apply(lambda s: np.log(s).mean()) larga = pd.DataFrame({ "unique_id": v["activo"], "ds": v.groupby("activo").cumcount(), "y": v["rendimiento"] / v["activo"].map(escala), "lrv": np.log(v["rv_5min"]) - v["activo"].map(c_lrv)}) nf = NeuralForecast(models=[DeepARNF( h=1, input_size=LV, lstm_n_layers=1, lstm_hidden_size=16, lstm_dropout=0.0, hist_exog_list=["lrv"], loss=DistributionLoss(distribution="StudentT", level=[90]), valid_loss=MQLoss(level=[90]), max_steps=300, val_check_steps=50, early_stop_patience_steps=3, batch_size=32, windows_batch_size=256, scaler_type=None, random_seed=0, accelerator="cpu", deterministic=True, logger=False, enable_progress_bar=False, enable_model_summary=False, enable_checkpointing=False)], freq=1) n_pru = pru.sum() // 6 cv = nf.cross_validation(larga, val_size=val.sum() // 6 + 1, test_size=n_pru, n_windows=None, step_size=1, refit=False) cv = cv.sort_values(["unique_id", "ds"]) assert np.allclose(cv["y"].to_numpy(), y, atol=1e-6) cob_nf = ((cv["y"] >= cv["DeepAR-lo-90"]) & (cv["y"] <= cv["DeepAR-hi-90"])).mean() print(f"neuralforecast DeepAR: 90% coverage {cob_nf:.1%}") ``` ```output neuralforecast DeepAR: 90% coverage 88.9% ``` The library version covers 88.9%: with 300 training steps it falls somewhat short, and its quantiles come from sampling the distribution, not from the formula. Tuning `max_steps`, the network size and the number of samples is the normal work with a library; the code of our own above serves to know what it does inside. ## Calibration: total and conditional Three measures for each model, on the 2018-2025 test set: - **Coverage of the 90% interval** and **exceedances of the 5% quantile** (the fraction of days below it; it should be 5%), with the p-value of the Kupiec test from lesson 4.2.3 (proportion of exceedances versus the promised one). - **Mean pinball loss** over the nine quantiles: the proper rule that summarizes the whole distribution. - **Mean width** of the 90% interval: at equal coverage, narrower is more informative. ```python def kupiec(excesos, n, p): x = excesos pi = x / n lr = -2 * ((n - x) * np.log(1 - p) + x * np.log(p) - (n - x) * np.log(1 - pi) - x * np.log(pi)) return stats.chi2.sf(lr, 1) def resumen(q): u = y[:, None] - q exc = (y < q[:, I05]).sum() return {"coverage 90": np.mean((y >= q[:, I05]) & (y <= q[:, I95])), "exceed. 5": exc / len(y), "p Kupiec": kupiec(exc, len(y), 0.05), "below 1%": np.mean(y < q[:, 0]), "pinball": np.maximum(TAU * u, (TAU - 1) * u).mean(), "width 90": np.mean(q[:, I95] - q[:, I05])} tabla = pd.DataFrame({k: resumen(q) for k, q in Q.items()}).T tabla.round(4) ``` ```output coverage 90 exceed. 5 p Kupiec below 1% pinball width 90 unconditional 0.9040 0.0480 0.3129 0.0113 0.1771 3.1451 HAR scale 0.8999 0.0522 0.2626 0.0105 0.1737 3.0368 oracle 0.8993 0.0516 0.4296 0.0116 0.1734 3.0311 GJR-GARCH-t 0.9021 0.0539 0.0523 0.0104 0.1755 3.0782 quantile network 0.9022 0.0478 0.2740 0.0177 0.1746 3.0642 DeepAR 0.9011 0.0554 0.0076 0.0107 0.1741 3.0504 ``` On **total** coverage they all look good: the six cover between 89.9% and 90.4%. Even the unconditional quantiles, which give the same interval on a calm day in 2019 as in March 2020. The pinball loss does rank them: the oracle is the best (0.1734), the HAR scale is almost level with it (0.1737), and DeepAR (0.1741), the quantile network (0.1746) and the GJR-GARCH (0.1755) sit between them and the unconditional quantiles (0.1771). The exceedances of the 5% quantile are close to 5% for all of them, although Kupiec rejects DeepAR's (5.5%, p = 0.008) and leaves the GARCH at the limit (p = 0.05); with 12,060 days, half a point is already detectable. And in the 1% tail, the quantile network leaves 1.8% of the days below: estimating an extreme quantile from its own data is hard. The way to see what fails is to split the test set by **volatility bucket**: the thirds of days according to the true expected variance (which no model uses; on real data, that day's realized RV). ```python tercio = pd.qcut(sig_ora[pru].to_numpy(), 3, labels=["low", "medium", "high"]) cob = pd.DataFrame({ k: pd.Series((y >= q[:, I05]) & (y <= q[:, I95])) .groupby(np.asarray(tercio)).mean() for k, q in Q.items()}).T[["low", "medium", "high"]] cob.round(3) ``` ```output low medium high unconditional 0.957 0.915 0.839 HAR scale 0.905 0.904 0.891 oracle 0.887 0.904 0.906 GJR-GARCH-t 0.938 0.905 0.864 quantile network 0.915 0.905 0.886 DeepAR 0.905 0.905 0.893 ``` The figure at the top draws it. The unconditional quantiles cover 95.7% of calm days and 83.9% of volatile ones: their 90% total is the average of an interval that is too wide and another that is too narrow, and it fails just when risk matters. The HAR scale, the quantile network and DeepAR stay within 1.5 points of 90% in the three thirds. The GJR-GARCH, which sees volatility only through squared returns, reacts more slowly: 93.8% on calm days and 86.4% on volatile ones. Not even the oracle nails the three thirds (88.7% in the low one): once the scale is known, the shape of standardized returns, estimated on 2008-2015, is not exactly the same in every bucket, and shocks (March 2020 in the index, November 2022 in crypto) cannot be forecast with the previous day's information. *Fig. — 90% interval forecast by the hand-written DeepAR for the crypto's daily return in the first half of 2020 (band) and observed returns (line; highlighted, the days outside).* A band that widens over the half-year: its mean width goes from 6.7 percentage points in January to 11.8 in June as volatility rises after the March shock. 13 of the 125 returns fall outside (10.4%), spread over the whole period: 4 in January, when the band is narrowest; the largest, +13.5% on 19 February. The figure shows DeepAR's 90% interval for crypto in the first half of 2020, converted to returns in percent. The volatility shock of March (a jump in variance, not an extreme return that day) does not show on the day it happens: the interval widens as the RV registers it, from a mean width of 6.7 points in January to 11.8 in June. Thirteen of the 125 days fall outside (10.4%), spread over the whole half-year; four in January, when the interval is narrowest, and the largest, a +13.5% on February 19. That is what an honest conditional model does: it reacts to information, it does not guess it. ## Conformal prediction: adjusting with the validation set **Conformal prediction** corrects an interval so that it covers what it promises, without assumptions about the model. The version for quantiles (CQR; Romano, Patterson and Candès, 2019): on a calibration set the model has not seen, you measure how far each data point falls outside the interval, $E_i = \max(\hat q_{lo} - y_i,\; y_i - \hat q_{hi})$ (negative if it falls inside), and widen (or narrow) the interval with the $\lceil (n+1)(1-\alpha) \rceil / n$ quantile of those values. If calibration and test are **exchangeable**, coverage on the test set is guaranteed. ```python yv = datos["val"][1].double().numpy() lo, hi = q_val_red[:, I05], q_val_red[:, I95] E = np.maximum(lo - yv, yv - hi) n = len(E) ajuste = np.quantile(E, np.ceil((n + 1) * 0.90) / n) q_red = Q["quantile network"] cob_val = np.mean((yv >= lo) & (yv <= hi)) cob_antes = np.mean((y >= q_red[:, I05]) & (y <= q_red[:, I95])) cob_cqr = np.mean((y >= q_red[:, I05] - ajuste) & (y <= q_red[:, I95] + ajuste)) print(f"validation coverage: {cob_val:.1%}; adjustment {ajuste:+.3f}") print(f"test: before {cob_antes:.1%}, after {cob_cqr:.1%}") ``` ```output validation coverage: 91.1%; adjustment -0.058 test: before 90.2%, after 89.1% ``` On the validation set the network covered somewhat more than 90%, so the adjustment is negative: it narrows the interval. On the test set, where the network already covered almost exactly 90%, the narrowed interval falls short. It is not a failure of the method: its guarantee requires validation and test to be alike, and 2016-2017 was a calmer period than 2018-2025. In time series, the **adaptive** variants (which recompute the adjustment with a rolling window of recent errors) come closer to the promise. The general lesson: no statistical correction replaces measuring coverage in the period and the buckets in which the model is going to be used. | Model | What it forecasts | Tails | Cost | |---|---|---|---| | Unconditional | fixed quantiles | empirical | none; useless for risk | | HAR scale + empirical shape | tomorrow's σ | standardized empirical | one regression | | GJR-GARCH-t | tomorrow's σ | Student's t | one estimation per asset | | Quantile network | each quantile | learned, poor at 1% | a small network | | DeepAR | μ, σ, ν per step | Student's t | a recurrent network | ## Exercises 1. Remove the monotone parameterization from the quantile network (have `forward` return `self.red(x)` as is), train it the same way and count how many pairs of consecutive quantiles cross on the test set. **Solution.** ```python class RedLibre(RedCuantiles): def forward(self, x): return self.red(x) torch.manual_seed(0) libre = RedLibre(len(FE)) entrenar(libre, pinball_media, datos) libre.eval() with torch.no_grad(): q_libre = libre(datos["prueba"][0]).double().numpy() cruzados = np.diff(q_libre, axis=1) <= 0 print("crossed pairs:", cruzados.sum(), "of", cruzados.size) print("days with some crossing:", f"{cruzados.any(axis=1).mean():.1%}") print("pinball:", round(resumen(q_libre)["pinball"], 4), "monotone:", round(tabla.loc["quantile network", "pinball"], 4)) ``` ```output crossed pairs: 0 of 96480 days with some crossing: 0.0% pinball: 0.174 monotone: 0.1746 ``` With these data and this network there is no crossing at all, and the pinball loss is almost the same (0.1740 versus 0.1746): the quantiles of a standardized return are well separated and the network learns to order them. Crossings appear with more quantiles closer together, less data or extrapolation to extreme features. The monotone parameterization costs nothing and removes the problem; fixing it afterward (sorting the outputs) also works, but it hides the fact that the network had learned something incoherent. 2. Repeat the table of 90% coverage by volatility third for crypto only, and compare the HAR scale with the GJR-GARCH-t. Which one adapts better and why? **Solution.** ```python es_cripto = (h.loc[pru, "activo"] == "cripto").to_numpy() t_c = pd.qcut(sig_ora[pru].to_numpy()[es_cripto], 3, labels=["low", "medium", "high"]) for k in ("HAR scale", "GJR-GARCH-t", "DeepAR"): q, yc = Q[k][es_cripto], y[es_cripto] dentro = pd.Series((yc >= q[:, I05]) & (yc <= q[:, I95])) print(f"{k:<12}", dentro.groupby(np.asarray(t_c)).mean() .reindex(["low", "medium", "high"]).round(3).to_dict()) ``` ```output HAR scale {'low': 0.887, 'medium': 0.918, 'high': 0.901} GJR-GARCH-t {'low': 0.936, 'medium': 0.936, 'high': 0.882} DeepAR {'low': 0.897, 'medium': 0.912, 'high': 0.897} ``` The HAR scale uses the 5-minute RV, a much more precise measure of today's volatility than a single day's squared return, which is all the GARCH sees. That is why it reacts sooner to level shifts (crypto has one in 2019 and frequent jumps) and its coverage varies less across thirds (from 88.7% to 91.8%, versus 88.2-93.6% for the GARCH, which also over-covers on calm and medium days). The exact difference depends on the asset; the general point is that a daily risk model improves more with better volatility data (intraday RV) than with a more complex architecture. 3. A colleague proposes measuring the quality of the 90% interval with total coverage alone and choosing the model whose coverage is closest to 90%. Which model would they choose with the lesson's table, and what would you tell them? **Solution.** With the lesson's table, the unconditional quantiles come close to 90% and could end up chosen: it is the model that gives the same interval every day. Total coverage is a necessary condition, not a sufficient one. Three more things have to be asked of it: conditional coverage (by volatility third, by asset, by regime), independence of exceedances (that they do not come in clusters, the Christoffersen test of 4.2.3) and a proper scoring rule (the mean pinball loss, or the CRPS if the whole distribution is forecast) that rewards information: at equal calibration, the narrowest interval is the best. With those criteria, the unconditional model comes last and the HAR scale, close to the oracle. ## Self-assessment 1. Which constant value minimizes the mean pinball loss with τ = 0.05 over a sample? - A. The sample mean - B. The median - C. The empirical 5% quantile of the sample - D. The mean minus 1.645 standard deviations *Answer C: The empirical 5% quantile of the sample.* The derivative of the loss with respect to the constant is the fraction of observations below it minus τ: it vanishes when 5% of the observations lie below. With τ = 0.5 it gives the median; the mean comes from squared error. 2. Some unconditional quantiles (the same every day) cover 90.4% on the test set. Are they well calibrated? - A. Yes: they cover what they promise - B. Only on average: they cover 95.7% on calm days and 83.9% on volatile ones, just when risk matters - C. No: they should cover 100% - D. Yes, if the Kupiec test does not reject *Answer B: Only on average: they cover 95.7% on calm days and 83.9% on volatile ones, just when risk matters.* Marginal coverage is not enough. A useful interval covers what it promises conditional on what is known today; the practical way to check it is to split the test set by volatility bucket (or by regime) and measure coverage in each one. 3. Why is the quantile network parameterized as the first quantile plus a cumulative sum of positive values? - A. So that it trains faster - B. Because the pinball loss requires it - C. So that the median is zero - D. So that the quantiles do not cross: with free outputs, the 5% quantile can end up above the 10% one on some days *Answer D: So that the quantiles do not cross: with free outputs, the 5% quantile can end up above the 10% one on some days.* Training each quantile with its own loss does not guarantee the order. `softplus` makes the jumps between consecutive quantiles positive and the cumulative sum orders them by construction. 4. The 1% quantile of the quantile network leaves 1.8% of the test days below it, and that of the GJR-GARCH-t, 1.0%. Why is it harder for the network? - A. Because a 1% quantile is learned from 1% of the data: with 12,000 rows, about 120 observations in the tail; the GARCH imposes a parametric form (Student's t) that extrapolates the tail using all the data - B. Because the network has few layers - C. Because the GARCH sees the future - D. Because the pinball loss does not work in the tails *Answer A: Because a 1% quantile is learned from 1% of the data: with 12,000 rows, about 120 observations in the tail; the GARCH imposes a parametric form (Student's t) that extrapolates the tail using all the data.* Tails are few observations. A model with a parametric tail shape, or a network that forecasts the parameters of a distribution (DeepAR), suffers less than one that estimates each extreme quantile separately. 5. The conformal adjustment with the 2016-2017 validation set takes the quantile network's coverage from 90.2% to 89.1% on the test set. Why does it not give exactly 90%? - A. Because the code is wrong - B. Because there is little test data - C. Because the guarantee of conformal prediction requires validation and test to be exchangeable, and a financial series changes regime: 2016-2017 was a calmer period than 2018-2025 - D. Because conformal prediction only works with linear models *Answer C: Because the guarantee of conformal prediction requires validation and test to be exchangeable, and a financial series changes regime: 2016-2017 was a calmer period than 2018-2025.* Conformal prediction is a guarantee without assumptions about the model, but not without assumptions about the data. In time series, its adaptive variants recompute the adjustment over a rolling window. ## Key points - The pinball loss ρτ(u) = max(τu, (τ − 1)u) is minimized at the τ quantile. It is a proper scoring rule: it rewards getting the distribution right, not just covering. - Total coverage is misleading: unconditional quantiles cover 90.4% overall, but 83.9% on the most volatile third of days. Calibration is measured conditionally (by volatility bucket). - With almost-HAR volatility data, the strongest benchmark is simple: the HAR σ times the empirical quantiles of standardized returns. It matches the oracle and beats the GARCH, the quantile network and DeepAR. - A quantile network needs monotonicity by construction and struggles in the extreme tails; DeepAR forecasts the parameters of a distribution per step and its tails come from the chosen shape (here, Student's t). - Conformal prediction corrects the width with the validation set, but in time series its guarantee depends on validation and test being alike. ## Resources ### Documentation - [Koenker and Bassett (1978): Regression Quantiles, Econometrica 46(1)](https://doi.org/10.2307/1913643) · doi.org - [Salinas, Flunkert and Gasthaus (2017, arXiv): DeepAR, Probabilistic Forecasting with Autoregressive Recurrent Networks](https://arxiv.org/abs/1704.04110) · arxiv.org - [Romano, Patterson and Candès (2019): Conformalized Quantile Regression](https://arxiv.org/abs/1905.03222) · arxiv.org - [Gneiting and Raftery (2007): Strictly Proper Scoring Rules, Prediction, and Estimation](https://doi.org/10.1198/016214506000001437) · doi.org - [arch: volatility models, forecasting and distributions](https://arch.readthedocs.io/en/latest/univariate/forecasting.html) · arch.readthedocs.io ### Data - [Level 5 data: returns and realized variance of six assets (Parquet)](https://pythonparatrading.com/curso/datos/nivel-5/vol5_diario.parquet) · Download · vol5_diario.parquet - [Level 5 data: true and forecast integrated variance (Parquet)](https://pythonparatrading.com/curso/datos/nivel-5/vol5_verdad.parquet) · Download · vol5_verdad.parquet Previous lesson: [5.7.2 · Time series architectures: TCN, N-BEATS, PatchTST and iTransformer](https://pythonparatrading.com/en/python-trading-course/level-5/tcn-n-beats-patchtst-itransformer/) Next lesson: [5.7.4 · Graph neural networks for asset relationships](https://pythonparatrading.com/en/python-trading-course/level-5/graph-neural-networks-for-assets/)