--- title: "5.7.3 Modelos probabilísticos: cuantiles y DeepAR" description: "Pérdida cuantílica, red de cuantiles monótona, DeepAR propio y de neuralforecast, cuantiles de GJR-GARCH y calibración condicional de los intervalos." url: "https://pythonparatrading.com/curso-python-trading/nivel-5/modelos-probabilisticos-cuantiles-y-deepar/" lang: "es" image: "https://pythonparatrading.com/og/curso-python-trading/nivel-5/modelos-probabilisticos-cuantiles-y-deepar.png" --- # 5.7.3 · Modelos probabilísticos: pérdida cuantílica, DeepAR e intervalos > Para dimensionar, poner un stop o calcular un VaR no basta con la previsión media: hace falta la distribución. Esta lección predice cuantiles del rendimiento de mañana con una red entrenada con la pérdida cuantílica y con un DeepAR, los compara con los cuantiles de un GJR-GARCH y de un oráculo, y enseña a comprobar que un intervalo del 90 % cubre el 90 % también en los días que importan. - Curso: [Nivel 5: Machine learning y deep learning (Keras/PyTorch)](https://pythonparatrading.com/curso-python-trading/nivel-5/) - Módulo: 5.7 · Deep learning con PyTorch - Duración: 75 min - Requisitos: [5.7.1 · PyTorch: tensores, autograd y el bucle de entrenamiento](https://pythonparatrading.com/curso-python-trading/nivel-5/pytorch-tensores-autograd-y-bucle-de-entrenamiento/), [4.2.3 · VaR, Expected Shortfall y backtesting del VaR](https://pythonparatrading.com/curso-python-trading/nivel-4/var-expected-shortfall-y-backtesting-del-var/), [2.1.4 · Modelos de volatilidad: ARCH, GARCH, EGARCH y GJR con la librería arch](https://pythonparatrading.com/curso-python-trading/nivel-2/volatilidad-arch-garch-egarch-gjr/) ## Al terminar esta lección - Explicar la pérdida cuantílica (pinball) y comprobar que su mínimo es el cuantil empírico. - Entrenar una red que prevé varios cuantiles a la vez sin que se crucen, y un DeepAR que prevé una distribución t de Student por paso. - Construir las referencias: cuantiles incondicionales, escala HAR, GJR-GARCH con `arch` y un oráculo con la varianza verdadera. - Medir la calibración: cobertura total, cobertura por tramos de volatilidad, excesos con el test de Kupiec y pérdida pinball media. - Ajustar un intervalo con predicción conforme sobre la validación y entender por qué en series no garantiza la cobertura. ## 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: | | Intervalo de confianza | Intervalo de predicción | |---|---|---| | Sobre qué | un parámetro (la media, una pendiente) | la próxima observación | | Qué recoge | la incertidumbre de la estimación | esa incertidumbre más el ruido del propio dato | | Con más datos | se estrecha hacia cero | se estrecha hasta la dispersión condicional, no más | | En esta lección | — | el intervalo del 90 % del rendimiento de mañana | 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. ```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 # 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}") ``` ```output 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, $$ \rho_\tau(u) = \max\big(\tau u,\; (\tau - 1)\, u\big). $$ 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. *Fig. — Pérdida pinball ρτ(u) = máx(τu, (τ − 1)u) para los cuantiles 5 %, 50 % y 95 %: el cuantil 5 % castiga 19 veces más quedarse por encima del dato que por debajo.* Tres funciones en forma de V asimétrica con el vértice en u = 0. Para τ = 0,05 (destacada), la rama izquierda sube con pendiente 0,95 y la derecha con 0,05; para τ = 0,5, las dos ramas tienen pendiente 0,5; para τ = 0,95 (discontinua), es la imagen especular de τ = 0,05. 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. ```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) # 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}") ``` ```output τ = 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. ```python 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()}) ``` ```output {'incondicional': (12060, 9), 'escala HAR': (12060, 9), 'oráculo': (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") # 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) ``` ```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 ``` 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. ```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("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) ``` ```output 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). ```mermaid 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). ```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): """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}") ``` ```output 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`). ```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: cobertura del 90 % {cob_nf:.1%}") ``` ```output 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. ```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 {"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) ``` ```output 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). ```python 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) ``` ```output 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. *Fig. — Intervalo del 90 % previsto por el DeepAR propio para el rendimiento diario de la cripto en el primer semestre de 2020 (banda) y rendimientos observados (línea; destacados, los días que quedan fuera).* Una banda que se ensancha a lo largo del semestre: su anchura media pasa de 6,7 puntos porcentuales en enero a 11,8 en junio, a medida que sube la volatilidad tras el choque de marzo. 13 de los 125 rendimientos quedan fuera (10,4 %), repartidos por todo el semestre: 4 en enero, cuando la banda es más estrecha; el mayor, +13,5 % el 19 de febrero. 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. ```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["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%}") ``` ```output 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. | Modelo | Qué prevé | Colas | Coste | |---|---|---|---| | Incondicional | cuantiles fijos | empíricas | ninguno; inútil en riesgo | | Escala HAR + forma empírica | σ de mañana | empíricas tipificadas | una regresión | | GJR-GARCH-t | σ de mañana | t de Student | una estimación por activo | | Red de cuantiles | cada cuantil | aprendidas, malas en el 1 % | una red pequeña | | DeepAR | μ, σ, ν por paso | t de Student | una red recurrente | ## 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. **Solución.** ```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("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)) ``` ```output pares cruzados: 0 de 96480 días con algún cruce: 0.0% pinball: 0.174 monótona: 0.1746 ``` Con estos datos y esta red no hay ningún cruce, y la pinball es casi la misma (0,1740 frente a 0,1746): los cuantiles de un rendimiento tipificado están bien separados y la red aprende a ordenarlos. Aparecen con más cuantiles, más cercanos entre sí, menos datos o extrapolación a features extremas. La parametrización monótona no cuesta nada y elimina el problema; arreglarlo después (ordenando las salidas) también funciona, pero esconde que la red había aprendido algo incoherente. 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é? **Solución.** ```python 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()) ``` ```output 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} ``` La escala HAR usa la RV de 5 minutos, una medida mucho más precisa de la volatilidad de hoy que el rendimiento al cuadrado de un día, que es lo único que ve el GARCH. Por eso reacciona antes a los cambios de nivel (la cripto tiene uno en 2019 y saltos frecuentes) y su cobertura varía menos entre tercios (del 88,7 % al 91,8 %, frente al 88,2-93,6 % del GARCH, que además cubre de más en los días tranquilos y medios). La diferencia exacta depende del activo; lo general es que un modelo de riesgo diario mejora más con mejores datos de volatilidad (RV intradía) que con una arquitectura más compleja. 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? **Solución.** Con la tabla de la lección, los incondicionales quedan cerca del 90 % y podrían salir elegidos: es el modelo que da el mismo intervalo todos los días. La cobertura total es una condición necesaria, no suficiente. Hay que pedirle tres cosas más: cobertura condicional (por tercios de volatilidad, por activo, por régimen), independencia de los excesos (que no lleguen en rachas, el test de Christoffersen de 4.2.3) y una regla de puntuación propia (la pinball media, o el CRPS si se prevé la distribución entera) que premie la información: a igual calibración, el intervalo más estrecho es el mejor. Con esos criterios, el incondicional queda el último y la escala HAR, cerca del oráculo. ## Autoevaluación 1. ¿Qué valor constante minimiza la pérdida pinball media con τ = 0,05 sobre una muestra? - A. La media de la muestra - B. La mediana - C. El cuantil 5 % empírico de la muestra - D. La media menos 1,645 desviaciones *Respuesta C: El cuantil 5 % empírico de la muestra.* La derivada de la pérdida respecto a la constante es la fracción de observaciones por debajo menos τ: se anula cuando el 5 % de las observaciones queda por debajo. Con τ = 0,5 da la mediana; la media sale del error cuadrático. 2. Unos cuantiles incondicionales (los mismos todos los días) cubren el 90,4 % en la prueba. ¿Están bien calibrados? - A. Sí: cubren lo que prometen - B. Solo en media: cubren el 95,7 % en los días tranquilos y el 83,9 % en los volátiles, justo cuando el riesgo importa - C. No: deberían cubrir el 100 % - D. Sí, si el test de Kupiec no rechaza *Respuesta B: Solo en media: cubren el 95,7 % en los días tranquilos y el 83,9 % en los volátiles, justo cuando el riesgo importa.* La cobertura marginal no basta. Un intervalo útil cubre lo prometido condicionado a lo que se sabe hoy; la forma práctica de comprobarlo es partir la prueba por tramos de volatilidad (o por régimen) y medir la cobertura en cada uno. 3. ¿Por qué la red de cuantiles se parametriza como el primer cuantil más una suma acumulada de valores positivos? - A. Para que entrene más deprisa - B. Porque la pérdida pinball lo exige - C. Para que la mediana sea cero - D. Para que los cuantiles no se crucen: con salidas libres, el cuantil del 5 % puede quedar por encima del del 10 % en algunos días *Respuesta D: Para que los cuantiles no se crucen: con salidas libres, el cuantil del 5 % puede quedar por encima del del 10 % en algunos días.* Entrenar cada cuantil con su pérdida no garantiza el orden. `softplus` hace positivos los saltos entre cuantiles consecutivos y la suma acumulada los ordena por construcción. 4. El cuantil del 1 % de la red de cuantiles deja por debajo el 1,8 % de los días de prueba, y el del GJR-GARCH-t, el 1,0 %. ¿Por qué cuesta más a la red? - A. Porque un cuantil del 1 % se aprende con el 1 % de los datos: con 12.000 filas, unas 120 observaciones en la cola; el GARCH impone una forma paramétrica (t de Student) que extrapola la cola con todos los datos - B. Porque la red tiene pocas capas - C. Porque el GARCH ve el futuro - D. Porque la pinball no funciona en las colas *Respuesta A: Porque un cuantil del 1 % se aprende con el 1 % de los datos: con 12.000 filas, unas 120 observaciones en la cola; el GARCH impone una forma paramétrica (t de Student) que extrapola la cola con todos los datos.* Las colas son pocas observaciones. Un modelo con forma paramétrica de la cola, o una red que prevé los parámetros de una distribución (DeepAR), sufre menos que una que estima cada cuantil extremo por separado. 5. El ajuste conforme con la validación de 2016-2017 lleva la cobertura de la red de cuantiles del 90,2 % al 89,1 % en la prueba. ¿Por qué no da exactamente el 90 %? - A. Porque el código está mal - B. Porque hay pocos datos en la prueba - C. Porque la garantía de la predicción conforme exige que la validación y la prueba sean intercambiables, y una serie financiera cambia de régimen: 2016-2017 fue un periodo más tranquilo que 2018-2025 - D. Porque la predicción conforme solo funciona con modelos lineales *Respuesta C: Porque la garantía de la predicción conforme exige que la validación y la prueba sean intercambiables, y una serie financiera cambia de régimen: 2016-2017 fue un periodo más tranquilo que 2018-2025.* La predicción conforme es una garantía sin supuestos sobre el modelo, pero no sin supuestos sobre los datos. En series, sus variantes adaptativas recalculan el ajuste con una ventana móvil. ## Ideas clave - La pérdida pinball ρτ(u) = máx(τu, (τ − 1)u) se minimiza en el cuantil τ. Es una regla de puntuación propia: premia acertar la distribución, no solo cubrir. - La cobertura total engaña: los cuantiles incondicionales cubren el 90,4 % en total, pero el 83,9 % en el tercio de días más volátiles. La calibración se mide condicionada (por tramos de volatilidad). - Con datos de volatilidad casi HAR, la referencia más fuerte es simple: σ del HAR por los cuantiles empíricos de los rendimientos tipificados. Queda a la altura del oráculo y por delante del GARCH, la red de cuantiles y DeepAR. - Una red de cuantiles necesita monotonía por construcción y sufre en las colas extremas; DeepAR prevé los parámetros de una distribución por paso y sus colas salen de la forma elegida (aquí, t de Student). - La predicción conforme corrige la anchura con la validación, pero en series su garantía depende de que validación y prueba se parezcan. ## Recursos ### Documentación - [Koenker y Bassett (1978): Regression Quantiles, Econometrica 46(1)](https://doi.org/10.2307/1913643) · doi.org - [Salinas, Flunkert y Gasthaus (2017, arXiv): DeepAR, Probabilistic Forecasting with Autoregressive Recurrent Networks](https://arxiv.org/abs/1704.04110) · arxiv.org - [Romano, Patterson y Candès (2019): Conformalized Quantile Regression](https://arxiv.org/abs/1905.03222) · arxiv.org - [Gneiting y Raftery (2007): Strictly Proper Scoring Rules, Prediction, and Estimation](https://doi.org/10.1198/016214506000001437) · doi.org - [arch: modelos de volatilidad, previsión y distribuciones](https://arch.readthedocs.io/en/latest/univariate/forecasting.html) · arch.readthedocs.io ### Datos - [Datos del nivel 5: rendimientos y varianza realizada de seis activos (Parquet)](https://pythonparatrading.com/curso/datos/nivel-5/vol5_diario.parquet) · Descarga · vol5_diario.parquet - [Datos del nivel 5: varianza integrada verdadera y prevista (Parquet)](https://pythonparatrading.com/curso/datos/nivel-5/vol5_verdad.parquet) · Descarga · vol5_verdad.parquet Lección anterior: [5.7.2 · Arquitecturas de series: TCN, N-BEATS, PatchTST e iTransformer](https://pythonparatrading.com/curso-python-trading/nivel-5/tcn-n-beats-patchtst-e-itransformer/) Lección siguiente: [5.7.4 · Redes de grafos para relaciones entre activos](https://pythonparatrading.com/curso-python-trading/nivel-5/redes-de-grafos-para-activos/)