Modelos probabilísticos: pérdida cuantílica, DeepAR e intervalos.

De un número a una distribución

Las lecciones 5.7.1 y 5.7.2 predijeron un número: la log RV de mañana. Muchas decisiones de trading necesitan más. Un VaR es un cuantil (lección 4.2.3); un stop razonable depende de cuánto se puede mover el precio en un día normal; el tamaño de una posición con volatility targeting, de la dispersión esperada. Un modelo probabilístico prevé la distribución del valor futuro, o al menos varios de sus cuantiles, y de ahí salen intervalos de predicción.

Dos intervalos que se confunden a menudo:

El problema: prever cuantiles del rendimiento de mañana de los seis activos de vol5_diario.parquet. El README de los datos lo dice claro: la media condicional del rendimiento es constante, así que no hay nada que prever en la media. Pero la dispersión cambia mucho y es previsible (es un HAR en el logaritmo de la varianza): un buen modelo probabilístico tiene que ensanchar el intervalo los días volátiles y estrecharlo los tranquilos. Datos sintéticos con la verdad conocida: vol5_verdad.parquet da la varianza integrada esperada de cada día, el mejor dato de escala posible.

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                       # posiciones de 5 % y 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()     # solo 2008-2015
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"curtosis de z_sig en la prueba: {stats.kurtosis(y):.1f}")

tramo
prueba    12060
ent       11958
val        3012
purga        12
curtosis de z_sig en la prueba: 3.5

Cada rendimiento se divide por la desviación de su activo en el entrenamiento (z), así que los seis comparten escala y se pueden modelar juntos, como en las dos lecciones anteriores. Las colas son gruesas: curtosis en exceso de 3,5 en la prueba, como en los rendimientos diarios reales. Las features son las del HAR (log RV de hoy y sus medias de 5 y 22 días, centradas por activo) más el rendimiento de hoy. Se prevén nueve cuantiles a la vez, del 1 al 99 %; el intervalo del 90 % va del cuantil 5 % al 95 %.

La pérdida cuantílica

Para entrenar un modelo que prevé el cuantil $\tau$ se usa la pérdida pinball (Koenker y Bassett, 1978): con $u = y - q$ el error,

Si el dato queda por encima del cuantil previsto ($u > 0$), cuesta $\tau u$; si queda por debajo, $(1 - \tau)|u|$. Para $\tau = 0{,}05$, quedarse corto por abajo cuesta 19 veces más que pasarse por arriba, y el mínimo se alcanza cuando solo el 5 % de los datos queda por debajo: en el cuantil. Con $\tau = 0{,}5$ es la mitad del error absoluto y da la mediana.

La forma de comprobarlo es la de la lección 5.7.1: buscar con autograd la constante que minimiza la pérdida media y compararla con el cuantil empírico.

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)   # tasa que baja a cero
    for _ in range(2000):
        opt.zero_grad()
        pinball(c, cripto, tau).backward()
        opt.step()
        paso.step()
    print(f"τ = {tau:.2f}: mínimo de la pinball {c.item():+.3f}, "
          f"cuantil empírico {np.quantile(cripto.numpy(), tau):+.3f}")

τ = 0.05: mínimo de la pinball -1.652, cuantil empírico -1.650
τ = 0.50: mínimo de la pinball +0.052, cuantil empírico +0.052
τ = 0.95: mínimo de la pinball +1.504, cuantil empírico +1.504

Coinciden; los dos milésimos del 5 % son la interpolación de np.quantile entre dos datos consecutivos, un tramo en el que la pinball es casi plana. La tasa de aprendizaje que baja a cero no es un adorno: el gradiente de la pinball va a saltos (cada dato aporta $-\tau$ o $1 - \tau$), y con una tasa fija Adam oscila alrededor del mínimo con una amplitud del orden de esa tasa. La pinball media sobre varios cuantiles es además una regla de puntuación propia (Gneiting y Raftery, 2007): en esperanza, ningún modelo puntúa mejor que el que da la distribución verdadera. Por eso sirve para comparar modelos probabilísticos, y la cobertura sola no: un intervalo puede cubrir el 90 % siendo inútil, como se verá.

Las referencias: incondicional, escala HAR, oráculo y GARCH

Antes de cualquier red, cuatro formas de prever cuantiles que marcan la escala:

• Incondicional: los cuantiles empíricos de z_sig en el entrenamiento, los mismos todos los días.
• Escala HAR: la volatilidad de mañana prevista con el HAR de la RV, $\hat\sigma_t = e^{\hat x_{t+1}/2}$, por los cuantiles empíricos de los rendimientos divididos por ella. Separa la escala (previsible) de la forma (estable).
• Oráculo de escala: lo mismo con la varianza integrada verdadera esperada, iv_prevista. Nadie la tiene; marca el suelo.
• GJR-GARCH(1,1) con t de Student (lección 2.1.4), estimado por activo con arch en 2008-2015 y aplicado sin reestimar: sus cuantiles son $\mu + \sigma_{t+1}\, t_\nu^{-1}(\tau)$, con la t tipificada a varianza 1.

def cuantiles_por_escala(sigma, en):
    """sigma (todas las filas) por los cuantiles empíricos de
    z_sig / sigma en el entrenamiento."""
    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 = {"incondicional": np.tile(np.quantile(h.loc[ent, "z_sig"], TAU),
                              (pru.sum(), 1)),
     "escala HAR": cuantiles_por_escala(sig_har, pru),
     "oráculo": cuantiles_por_escala(sig_ora, pru)}
print({k: q.shape for k, q in Q.items()})

{'incondicional': (12060, 9), 'escala HAR': (12060, 9), 'oráculo': (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")
    # previsiones hechas en cada fecha t >= 2018 para 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

Con last_obs el modelo se estima con los datos anteriores a 2016; forecast(start=…) aplica esos parámetros filtrando todos los datos y devuelve, para cada fecha $t$, la varianza prevista para $t + 1$, que es exactamente lo que se necesita: la fila de la fecha $t$ de la prueba con su objetivo de mañana. Los rendimientos van en tanto por ciento (multiplicados por 100) porque el optimizador de arch converge mejor con esa escala, y ppf de su distribución t ya está tipificada a varianza 1. Las persistencias ($\alpha + \gamma/2 + \beta$) rondan 0,9-0,99 y los grados de libertad, entre 3,3 y 5,3: colas gruesas en todos.

Una red de cuantiles

La forma directa de prever cuantiles con deep learning: una red con una salida por cuantil, entrenada con la suma de las pérdidas pinball. Tiene un problema conocido, el cruce de cuantiles: nada obliga a que el cuantil del 5 % quede por debajo del del 10 %. La solución es construirlo: la red da el primer cuantil y los saltos entre cuantiles consecutivos, que pasan por softplus (siempre positivos) y se acumulan.

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("mejor época:", entrenar(red_q, pinball_media, datos))
red_q.eval()
with torch.no_grad():
    Q["red de cuantiles"] = red_q(datos["prueba"][0]).double().numpy()
    q_val_red = red_q(datos["val"][0]).double().numpy()
cruces = (np.diff(Q["red de cuantiles"], axis=1) <= 0).sum()
print("cuantiles cruzados en la prueba:", cruces)

mejor época: 27
cuantiles cruzados en la prueba: 0

Cero cruces, por construcción. El precio de la monotonía es nulo: la red puede representar cualquier conjunto de cuantiles ordenados.

DeepAR: una distribución por paso

DeepAR (Salinas, Flunkert y Gasthaus, 2017; Amazon) es una red recurrente que, en cada paso de la serie, recibe el valor anterior y las covariables y devuelve los parámetros de una distribución para el valor siguiente. Se entrena maximizando la verosimilitud de lo observado en todos los pasos de la ventana, y se usa de dos formas: para un paso, leyendo los cuantiles de la distribución prevista; para varios, muestreando un valor, metiéndolo como entrada del paso siguiente y repitiendo (trayectorias de Monte Carlo).

flowchart TB
  E[z hoy y log RV] --> R[LSTM]
  H[Estado previo] --> R
  R --> P[μ, σ, ν]
  P --> D[t de Student]
  D --> C[Cuantiles]
  D -. muestreo .-> E

Las diferencias con la red de cuantiles: la forma de la distribución la elige quien modela (aquí una t de Student, que deja aprender colas gruesas con un parámetro ν), y los cuantiles extremos salen de esa forma, no de estimar cada uno por separado con los pocos datos de la cola. La versión propia de abajo usa ventanas de 32 días, un LSTM de 16 unidades y la log RV como covariable; es pequeña para que entrene en segundos (en producción, más unidades, más épocas y varias semillas).

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):
    """Ventanas de LV días (z y log RV) que terminan en cada fila de
    `en`, y el objetivo de cada paso (el z del día siguiente)."""
    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):                             # todos los pasos
    mu, sigma, nu = salida
    t = torch.distributions.StudentT(nu, mu, sigma)
    return -t.log_prob(y).mean()

def nll_ultimo(salida, y):                      # solo el último
    return nll([p[:, -1] for p in salida], y[:, -1])

torch.manual_seed(0)
deepar = DeepAR()
print("mejor época:", 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"ν medio en la prueba: {nu.mean():.1f}")

mejor época: 17
ν medio en la prueba: 4.6

• El entrenamiento usa la verosimilitud de los 32 pasos de cada ventana (así entrena DeepAR); la validación, solo la del último, que es el que se usa para prever y el único cuyo objetivo cae entero en 2016-2017.
• torch.distributions.StudentT da la log-verosimilitud con gradiente; no tiene función cuantil, así que los cuantiles salen de scipy.stats.t.ppf con los parámetros previstos. Con una normal, Normal.icdf lo haría en torch.
• El ν medio que aprende ronda 4-5, en línea con los grados de libertad que estima el GARCH.

La versión de librería es neuralforecast.models.DeepAR, con el mismo formato de datos que la lección 5.7.2, la distribución en la pérdida (DistributionLoss) y la log RV como covariable histórica (hist_exog_list): solo se conoce hasta hoy. Sus intervalos salen de muestrear la distribución prevista (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: cobertura del 90 % {cob_nf:.1%}")

neuralforecast DeepAR: cobertura del 90 % 88.9%

La versión de librería cubre el 88,9 %: con 300 pasos de entrenamiento queda algo corta, y sus cuantiles salen de muestrear la distribución, no de la fórmula. Ajustar max_steps, el tamaño de la red y las muestras es el trabajo normal con una librería; el código propio de arriba sirve para saber qué hace por dentro.

Calibración: total y condicional

Tres medidas para cada modelo, en la prueba 2018-2025:

• Cobertura del intervalo del 90 % y excesos del cuantil 5 % (la fracción de días por debajo; debería ser 5 %), con el p-valor del test de Kupiec de la lección 4.2.3 (proporción de excesos frente a la prometida).
• Pinball media sobre los nueve cuantiles: la regla propia que resume la distribución entera.
• Anchura media del intervalo del 90 %: a igual cobertura, más estrecho es más informativo.

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 {"cobertura 90": np.mean((y >= q[:, I05])
                                    & (y <= q[:, I95])),
            "excesos 5": exc / len(y),
            "p Kupiec": kupiec(exc, len(y), 0.05),
            "debajo del 1 %": np.mean(y < q[:, 0]),
            "pinball": np.maximum(TAU * u, (TAU - 1) * u).mean(),
            "anchura 90": np.mean(q[:, I95] - q[:, I05])}

tabla = pd.DataFrame({k: resumen(q) for k, q in Q.items()}).T
tabla.round(4)

                  cobertura 90  excesos 5  p Kupiec  debajo del 1 %  pinball  anchura 90
incondicional           0.9040     0.0480    0.3129          0.0113   0.1771      3.1451
escala HAR              0.8999     0.0522    0.2626          0.0105   0.1737      3.0368
oráculo                 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
red de cuantiles        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

En la cobertura total todos parecen buenos: los seis cubren entre el 89,9 % y el 90,4 %. Hasta los cuantiles incondicionales, que dan el mismo intervalo un día tranquilo de 2019 que en marzo de 2020. La pinball sí ordena: el oráculo es el mejor (0,1734), la escala HAR queda casi a su altura (0,1737), y DeepAR (0,1741), la red de cuantiles (0,1746) y el GJR-GARCH (0,1755) quedan entre ellos y los incondicionales (0,1771). Los excesos del cuantil 5 % están cerca del 5 % en todos, aunque Kupiec rechaza los de DeepAR (5,5 %, p = 0,008) y deja al GARCH en el límite (p = 0,05); con 12.060 días, medio punto ya es detectable. Y en la cola del 1 %, la red de cuantiles deja por debajo el 1,8 % de los días: estimar un cuantil extremo con sus propios datos es difícil. La forma de ver qué falla es partir la prueba por tramos de volatilidad: los tercios de días según la varianza verdadera esperada (que ningún modelo usa; en datos reales, la RV realizada de ese día).

tercio = pd.qcut(sig_ora[pru].to_numpy(), 3,
                 labels=["baja", "media", "alta"])
cob = pd.DataFrame({
    k: pd.Series((y >= q[:, I05]) & (y <= q[:, I95]))
       .groupby(np.asarray(tercio)).mean()
    for k, q in Q.items()}).T[["baja", "media", "alta"]]
cob.round(3)

                   baja  media   alta
incondicional     0.957  0.915  0.839
escala HAR        0.905  0.904  0.891
oráculo           0.887  0.904  0.906
GJR-GARCH-t       0.938  0.905  0.864
red de cuantiles  0.915  0.905  0.886
DeepAR            0.905  0.905  0.893

La figura de la cabecera lo dibuja. Los cuantiles incondicionales cubren el 95,7 % de los días tranquilos y el 83,9 % de los volátiles: su 90 % total es la media de un intervalo demasiado ancho y otro demasiado estrecho, y falla justo cuando el riesgo importa. La escala HAR, la red de cuantiles y DeepAR quedan a menos de 1,5 puntos del 90 % en los tres tercios. El GJR-GARCH, que ve la volatilidad solo a través de los rendimientos al cuadrado, reacciona más despacio: 93,8 % en los días tranquilos y 86,4 % en los volátiles. Ni el oráculo clava los tres tercios (88,7 % en el bajo): conocida la escala, la forma de los rendimientos tipificados, estimada en 2008-2015, no es exactamente la misma en todos los tramos, y los choques (marzo de 2020 en el índice, noviembre de 2022 en la cripto) no se pueden prever con la información del día anterior.

La figura muestra el intervalo del 90 % de DeepAR para la cripto en el primer semestre de 2020, pasado a rendimientos en tanto por ciento. El choque de volatilidad de marzo (un salto en la varianza, no un rendimiento extremo ese día) no se ve el día en que ocurre: el intervalo se va ensanchando a medida que la RV lo registra, de 6,7 puntos de anchura media en enero a 11,8 en junio. Trece de los 125 días quedan fuera (10,4 %), repartidos por todo el semestre; cuatro en enero, cuando el intervalo es más estrecho, y el mayor, un +13,5 % el 19 de febrero. Es lo que hace un modelo condicional honesto: reacciona a la información, no la adivina.

Predicción conforme: ajustar con la validación

La predicción conforme corrige un intervalo para que cubra lo prometido sin supuestos sobre el modelo. La versión para cuantiles (CQR; Romano, Patterson y Candès, 2019): en un conjunto de calibración que el modelo no ha visto, se mide cuánto se sale cada dato del intervalo, $E_i = \max(\hat q_{lo} - y_i,\; y_i - \hat q_{hi})$ (negativo si queda dentro), y se ensancha (o estrecha) el intervalo con el cuantil $\lceil (n+1)(1-\alpha) \rceil / n$ de esos valores. Si calibración y prueba son intercambiables, la cobertura de la prueba queda garantizada.

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["red de cuantiles"]
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"cobertura en validación: {cob_val:.1%}; ajuste {ajuste:+.3f}")
print(f"prueba: antes {cob_antes:.1%}, después {cob_cqr:.1%}")

cobertura en validación: 91.1%; ajuste -0.058
prueba: antes 90.2%, después 89.1%

En la validación la red cubría algo más del 90 %, así que el ajuste es negativo: estrecha el intervalo. En la prueba, donde la red ya cubría casi exactamente el 90 %, el intervalo estrechado se queda por debajo. No es un fallo del método: su garantía exige que validación y prueba se parezcan, y 2016-2017 fue un periodo más tranquilo que 2018-2025. En series, las variantes adaptativas (que recalculan el ajuste con una ventana móvil de los errores recientes) se acercan más a la promesa. La lección general: ninguna corrección estadística sustituye a medir la cobertura en el periodo y en los tramos en que se va a usar el modelo.

Ejercicios

1. Quita la parametrización monótona de la red de cuantiles (que forward devuelva self.red(x) tal cual), entrénala igual y cuenta cuántos pares de cuantiles consecutivos se cruzan en la prueba.

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("pares cruzados:", cruzados.sum(), "de", cruzados.size)
print("días con algún cruce:", f"{cruzados.any(axis=1).mean():.1%}")
print("pinball:", round(resumen(q_libre)["pinball"], 4),
      "monótona:", round(tabla.loc["red de cuantiles", "pinball"], 4))

pares cruzados: 0 de 96480
días con algún cruce: 0.0%
pinball: 0.174 monótona: 0.1746

2. Repite la tabla de cobertura del 90 % por tercios de volatilidad solo para la cripto, y compara la escala HAR con el GJR-GARCH-t. ¿Cuál se adapta mejor y por qué?

es_cripto = (h.loc[pru, "activo"] == "cripto").to_numpy()
t_c = pd.qcut(sig_ora[pru].to_numpy()[es_cripto], 3,
              labels=["baja", "media", "alta"])
for k in ("escala HAR", "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(["baja", "media", "alta"]).round(3).to_dict())

escala HAR   {'baja': 0.887, 'media': 0.918, 'alta': 0.901}
GJR-GARCH-t  {'baja': 0.936, 'media': 0.936, 'alta': 0.882}
DeepAR       {'baja': 0.897, 'media': 0.912, 'alta': 0.897}

3. Un compañero propone medir la calidad del intervalo del 90 % solo con la cobertura total y elegir el modelo cuya cobertura quede más cerca del 90 %. ¿Qué modelo elegiría con la tabla de la lección y qué le contestarías?