Probabilistic models: quantile loss, DeepAR and intervals.

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:

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.

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}")

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,

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.

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.

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}")

τ = 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.

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()})

{'unconditional': (12060, 9), 'HAR scale': (12060, 9), 'oracle': (12060, 9)}

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)

        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.

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)

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).

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).

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}")

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).

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%}")

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.

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)

                  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).

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)

                    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.

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.

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%}")

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.

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.

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))

crossed pairs: 0 of 96480
days with some crossing: 0.0%
pinball: 0.174 monotone: 0.1746

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?

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())

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}

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?