Puntuación de pronósticos probabilísticos: CRPS, calibración PIT y DeepAR

Este blog ya ha defendido la predicción de distribuciones en lugar de puntos, y ya ha resuelto qué hacer con el intervalo una vez que lo tenga: predicción conforme para el dimensionamiento de posiciones conscientes del riesgo deriva intervalos libres de distribución y la regla de dimensionamiento que los consume, y el artículo Temporal Fusion Transformer envía un capa de salida multicuantil. Lo que ninguno de los dos cubre es la parte que decide si vale la pena confiar en algo: cómo se califica una distribución predictiva y cómo se verifica que la incertidumbre declarada sea honesta.

De eso se trata este artículo. Tres cosas concretamente:

1. CRPS: la regla de puntuación adecuada para pronósticos probabilísticos y su relación exacta con la pérdida de pinball que ya informa el artículo de TFT.
2. El histograma PIT: un diagnóstico de calibración que indica cómo un modelo está mal calibrado, no solo si lo está.
3. DeepAR: la familia de modelos de muestreo autorregresivo, que aparece en este blog solo como una nota al pie de referencia y nunca se explica.

Primero una pieza del marco, porque determina qué maquinaria necesita. Un pronóstico puntual falla de manera diferente según el régimen: en una tendencia de baja volatilidad la distribución condicional es estrecha y casi simétrica y un pronóstico puntual es un buen resumen; antes de un evento programado es bimodal, y la media condicional se ubica exactamente donde el precio tiene menos probabilidades de aterrizar; en una crisis la cola izquierda domina y la media subestima gravemente la desventaja. Tres rutas recuperan la distribución completa: paramétrica (predice los parámetros de una familia supuesta: rápido, riesgo de especificación errónea), basada en cuantiles (predice una cuadrícula fija: libre de suposiciones, discreta) y basada en muestras (rutas de Monte Carlo a partir de muestreo autorregresivo, abandono o conjuntos: flexibles, costosas). Los dos últimos dominan en las finanzas precisamente porque la distribución cambia de forma entre esos regímenes.

Pronósticos cuantiles, brevemente

La pérdida de pinball y una cuadrícula de cuantiles concreta ya se proporcionan en la [capa de salida de cuantiles TFT] (/en/blog/post/temporal-fusion-transformer-trading), por lo que no reformularé la fórmula. Lo único que vale la pena interiorizar es su asimetría: al menos $\tau = 0.5$ La sobrepredicción y la subpredicción cuestan lo mismo y la pérdida se reduce a MAE, pero al menos $\tau = 0.95$ la subpredicción es penalizada $19\times$ más que una sobrepredicción. esa proporción $\tau/(1-\tau)$ es el mecanismo: es lo que arrastra el valor ajustado hasta un nivel que sólo el 5% de las observaciones debería superar.

El problema práctico que el artículo publicado no menciona es el cruce de cuantiles: nada en una pérdida por cuantil evita $\hat{q}_{0.25} > \hat{q}_{0.75}$. Con suficientes datos y una red troncal compartida, esto es raro, pero ocurre en los cuantiles extremos de muestras delgadas, y silenciosamente rompe cualquier CRPS posterior o cálculo de cobertura. La solución barata es una especie post-hoc del vector cuantil predicho; el de principios es una parametrización de salida monótona (predecir $q_{\min}$ más incrementos no negativos).

Una implementación mínima desde cero, útil cuando desea resultados cuantiles sin adoptar un marco de pronóstico:

import torch
import torch.nn as nn

class QuantileRegressionNet(nn.Module):
    """Multi-quantile forecasting network for financial returns."""

    def __init__(self, input_dim: int, hidden_dim: int = 128,
                 quantiles: list[float] = [0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99]):
        super().__init__()
        self.quantiles = quantiles
        self.backbone = nn.Sequential(
            nn.Linear(input_dim, hidden_dim), nn.ReLU(), nn.Dropout(0.2),
            nn.Linear(hidden_dim, hidden_dim), nn.ReLU(), nn.Dropout(0.2),
        )
        self.heads = nn.ModuleList([nn.Linear(hidden_dim, 1) for _ in quantiles])

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        h = self.backbone(x)
        out = torch.cat([head(h) for head in self.heads], dim=-1)
        return torch.sort(out, dim=-1).values


def pinball_loss(predictions: torch.Tensor, targets: torch.Tensor,
                 quantiles: list[float]) -> torch.Tensor:
    """predictions: (batch, n_quantiles); targets: (batch, 1)."""
    errors = targets - predictions
    tau = torch.tensor(quantiles, device=predictions.device).unsqueeze(0)
    return torch.max(tau * errors, (tau - 1) * errors).mean()

El ancho del intervalo se asigna al tamaño de la posición, y trabajamos esa asignación en su totalidad, incluido el filtro de no comercio de relación de borde, en predicción conforme para el tamaño de la posición consciente del riesgo; el $1/\sigma$ la equivalencia con el objetivo de volatilidad se deriva en objetivo de volatilidad con pronósticos GARCH.

DeepAR: pronóstico probabilístico autorregresivo

DeepAR (Salinas et al., 2020) toma la ruta basada en muestras. En lugar de predecir cuantiles directamente, utiliza un RNN autorregresivo para parametrizar una probabilidad en cada paso y luego dibuja distribuciones predictivas completas haciendo avanzar esa probabilidad.

Arquitectura. En cada paso $t$ la red consume la observación anterior $z_{t-1}$ (escalado por un factor por serie), covariables $\mathbf{x}_t$y funciones estáticas opcionales. Un LSTM lleva el estado:

$\mathbf{h}_t = \text{LSTM}(\mathbf{h}_{t-1}, [z_{t-1}, \mathbf{x}_t])$

Una cabeza densa mapas $\mathbf{h}_t$ a los parámetros de probabilidad $\theta_t = \text{MLP}(\mathbf{h}_t)$ — $(\mu_t, \sigma_t)$ para gaussiano, $(\mu_t, \sigma_t, \nu_t)$ para Estudiante-t. El entrenamiento maximiza $\sum_i \sum_t \log p(z_{i,t} \mid \theta_{i,t})$ en todas las series. En la inferencia usted muestra de $p(\cdot \mid \theta_t)$, retroalimente la muestra como la siguiente entrada y repita $S$ tiempos para conseguir $S$ trayectorias.

Dos consecuencias son importantes para el comercio. La probabilidad es una elección de modelo, por lo que las colas gruesas son algo que usted selecciona en lugar de algo que espera: Student-t para colas pesadas o una mezcla gaussiana. $\sum_k w_k \mathcal{N}(z_t; \mu_k(\mathbf{h}_t), \sigma_k^2(\mathbf{h}_t))$ para el comportamiento de eventos bimodales. Los pronósticos de varios pasos son coherentes: cada ruta de muestra es una trayectoria plausible que preserva la correlación serial, que es lo que se necesita para cualquier horizonte de más de un paso. (El tercer punto de venta habitual: un modelo global en muchas series supera $N$ modelos por activo: es el mismo argumento de eficiencia de datos presentado en el artículo de TFT, con el escalado por serie y las covariables estáticas desempeñando el papel de codificadores estáticos allí).

DeepAR con GluonTS

import pandas as pd
from gluonts.dataset.pandas import PandasDataset
from gluonts.torch.model.deepar import DeepAREstimator
from gluonts.torch.distributions import StudentTOutput
from gluonts.evaluation import make_evaluation_predictions, Evaluator


def prepare_crypto_dataset(returns_df: pd.DataFrame, freq: str = "h"):
    """Wide return frame (index=datetime, columns=assets) -> GluonTS dataset.

    from_long_dataframe wants LONG format: one row per (timestamp, item_id)
    with the target in a column. Passing a DataFrame of dicts does not work.
    """
    long_df = (
        returns_df.stack()
        .rename("target")
        .rename_axis(index=["timestamp", "item_id"])
        .reset_index()
        .dropna(subset=["target"])
    )
    return PandasDataset.from_long_dataframe(
        long_df, target="target", item_id="item_id",
        timestamp="timestamp", freq=freq,
    )


estimator = DeepAREstimator(
    prediction_length=24,    # 24 hours ahead
    context_length=168,      # one week of hourly data
    freq="h",
    num_layers=2,
    hidden_size=64,
    dropout_rate=0.1,
    lr=1e-3,
    batch_size=64,
    trainer_kwargs={"max_epochs": 50},
    distr_output=StudentTOutput(),
)

predictor = estimator.train(training_data=train_dataset)

forecast_it, ts_it = make_evaluation_predictions(
    dataset=test_dataset, predictor=predictor, num_samples=500,
)
forecasts, actuals = list(forecast_it), list(ts_it)

evaluator = Evaluator(quantiles=[0.05, 0.25, 0.5, 0.75, 0.95])
agg_metrics, item_metrics = evaluator(actuals, forecasts)
print(f"mean_wQuantileLoss: {agg_metrics['mean_wQuantileLoss']:.4f}")

num_samples controla el recuento de rutas de Monte Carlo. Para uso en vivo, normalmente son suficientes 100-200; para la evaluación fuera de línea utilice 500-1000, ya que el CRPS estimado a partir de muy pocas muestras tiene un sesgo bajo.

Otras rutas hacia una distribución predictiva

Vale la pena conocer tres alternativas pero que no necesitan apartados propios. MC Dropout (Gal y Ghahramani, 2016) mantiene el abandono activo en la inferencia y toma la media y la varianza $S$ pases hacia adelante; es una inferencia variacional aproximada disfrazada y es la forma más rápida de incorporar incertidumbre a un modelo que ya entrenó. Conjuntos profundos (Lakshminarayanan et al., 2017) entrenan $M$ copias de diferentes semillas y tratar la distribución predictiva como una mezcla: consistentemente fuerte en los puntos de referencia, al menos $M\times$ el coste, lo que lo descarta para horizontes sensibles a la latencia pero no para 4h ni diario. Normalizing flows (Rasul et al., 2021) aprende un mapa invertible desde una densidad base simple hasta un objetivo arbitrario, capturando la multimodalidad y la asimetría sin elegir una familia paramétrica. Los tres producen muestras, por lo que todo lo que se describe en la siguiente sección se aplica a ellos sin cambios.

CRPS: la métrica adecuada para una distribución predictiva

Pronósticos de puntos MSE y MAE. No pueden decirte si una distribución fue buena, porque solo miran un resumen de la misma. El reemplazo estándar es el Puntuación de probabilidad clasificada continua, y es lo más útil de este artículo.

CRPS es la distancia al cuadrado integrada entre la CDF predicha y la CDF degenerada que pone toda la masa en lo que realmente sucedió:

$\text{CRPS}(F, y) = \int_{-\infty}^{\infty} \left(F(z) - \mathbb{1}[y \leq z]\right)^2 dz$

Tres propiedades le ganan su lugar:

• Es una regla de puntuación adecuada (Gneiting y Raftery, 2007): minimizada en la expectativa sólo cuando la distribución predicha es igual a la verdadera. Ni el exceso de confianza deliberado ni la cobertura con una distribución artificialmente amplia mejoran su puntuación. MSE en la mediana no tiene esa propiedad, por lo que CRPS no es opcional.
• Generaliza MAE. Para un pronóstico puntual degenerado, colapsa hasta un error absoluto, por lo que CRPS vive en las unidades del objetivo: al pronosticar retornos de registros horarios, un CRPS de 0,004 es directamente comparable a un error absoluto medio de 40 pb, en lugar de un número sin unidades que no se puede verificar.
• Recompensa la nitidez sujeta a calibración. Entre dos pronósticos igualmente calibrados, el más estrecho obtiene mejor puntuación. Ésa es también la razón por la que el CRPS por sí solo es insuficiente: una mala puntuación no indica cuál de las dos condiciones falló, y para eso sirve la siguiente sección.

El puente hacia la pérdida cuantil

Esta es la conexión que falta en el resto del blog. Dada una grilla cuantil $\mathcal{T}$ en lugar de una CDF completa, CRPS se aproxima mediante la pérdida de pinball:

$\text{CRPS} \approx \frac{2}{|\mathcal{T}|} \sum_{\tau \in \mathcal{T}} \mathcal{L}_\tau(y, \hat{q}_\tau)$

Léalo literalmente: la "pérdida cuantil" informada en el artículo de TFT y el CRPS analizado aquí son la misma cantidad hasta un factor de 2, y la aproximación se estrecha a medida que la red se densifica. No son métricas competitivas y no hay motivo para informar ambas.

Una unidad de advertencia, porque muerde. GluonTS mean_wQuantileLoss es esa misma pérdida promedio de pinball normalizada por la suma de valores objetivo absolutos, lo que la hace adimensional y comparable entre activos de diferente escala. CRPS propiamente dicho no está normalizado y permanece en unidades de retorno. no imprimir mean_wQuantileLoss debajo de una etiqueta que dice "CRPS": la versión preliminar de este artículo hizo exactamente eso, y es una manera fácil de comparar dos números que no están en la misma escala.

CRPS de muestras

Con muestras de Montecarlo $\{y^{(s)}\}_{s=1}^{S}$ (DeepAR, conjuntos, flujos), utilice la forma de energía:

$\text{CRPS}(F, y) = \frac{1}{S} \sum_{s=1}^{S} |y^{(s)} - y| - \frac{1}{2S^2} \sum_{s=1}^{S} \sum_{s'=1}^{S} |y^{(s)} - y^{(s')}|$

El primer término premia la precisión, el segundo penaliza la dispersión excesiva. Escrito ingenuamente, el segundo término es $O(S^2)$ y se convierte en el cuello de botella cuando se obtienen miles de pronósticos. La clasificación primero lo colapsa: para estadísticas de orden ascendente $y_{(1)} \le \dots \le y_{(S)}$, la doble suma es igual $2\sum_k (2k - S - 1)\, y_{(k)}$, entonces todo es $O(S \log S)$ dominado por el género.

import numpy as np

def crps_empirical(samples: np.ndarray, observation: float) -> float:
    """CRPS from Monte Carlo samples, O(n log n) via order statistics."""
    n = len(samples)
    mae = np.mean(np.abs(samples - observation))
    x = np.sort(samples)
    k = np.arange(1, n + 1)
    dispersion = np.sum((2 * k - n - 1) * x) / n**2
    return mae - dispersion


def crps_quantile(quantile_predictions: np.ndarray,
                  quantile_levels: np.ndarray,
                  observation: float) -> float:
    """CRPS approximation from a quantile grid (2x mean pinball loss)."""
    errors = observation - quantile_predictions
    pinball = np.where(errors >= 0,
                       quantile_levels * errors,
                       (quantile_levels - 1) * errors)
    return 2.0 * np.mean(pinball)

Ambas formas deberían coincidir estrechamente en la misma distribución predictiva; si no lo hacen, sospeche de un cruce de cuantiles o de muy pocas muestras. En producción, properscoring.crps_ensemble(observation, samples) es un complemento bien probado.

Calibración: ¿Es honesta la incertidumbre declarada?

Un buen CRPS no garantiza que los intervalos signifiquen lo que dicen. Un modelo cuyo intervalo nominal del 90 % captura el 70 % de los resultados es demasiado confiado, y en una regla de tamaño que lee el ancho del intervalo, lo impulsará precisamente cuando no debería hacerlo. El artículo de TFT establece la advertencia y prescribe la predicción conforme como solución; lo que sigue es cómo se detecta y diagnostica realmente la falla.

El histograma PIT

Para cada observación $y_t$, calcule su cuantil según la CDF prevista de ese paso:

$u_t = \hat{F}_t(y_t)$

Si el modelo está calibrado, $\{u_t\}$ es uniforme en $[0,1]$. El poder de diagnóstico es que la forma de la desviación nombra el modo de falla:

• En forma de U: demasiado confiado: demasiada masa cae en las colas, la distribución es demasiado estrecha.
• En forma de joroba: poco seguro: las observaciones se agrupan cerca del centro, la distribución es demasiado amplia.
• Sesgado a la izquierda: el modelo sobrepredice sistemáticamente.
• Sesgado a la derecha: el modelo subpredice sistemáticamente.

Una nota para evitar confusiones con modelos de cópula para riesgo conjunto, que también utiliza la transformación integral de probabilidad: ahí hay una transformación marginal, un paso de preprocesamiento que convierte los marginales de GARCH-EVT en observaciones pseudouniformes para que se les pueda ajustar una cópula. Aquí la transformación se aplica a pronósticos fuera de la muestra y la uniformidad es el resultado que se prueba, no un insumo que se fabrica. Mismas matemáticas, dirección opuesta de inferencia.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import kstest

def pit_calibration_check(forecasts, actuals, n_bins=20):
    """PIT histogram + KS test against uniform.

    forecasts: list of SampleForecast objects (GluonTS)
    """
    pit_values = []
    for forecast, actual in zip(forecasts, actuals):
        samples = forecast.samples          # (n_samples, prediction_length)
        h = forecast.prediction_length
        for t in range(h):
            obs = actual.values[-h + t]
            pit_values.append(np.mean(samples[:, t] <= obs))

    pit_values = np.array(pit_values)
    ks_stat, p_value = kstest(pit_values, "uniform")

    fig, ax = plt.subplots(figsize=(8, 4))
    ax.hist(pit_values, bins=n_bins, density=True, alpha=0.7, edgecolor="black")
    ax.axhline(y=1.0, color="red", linestyle="--", label="Perfect calibration")
    ax.set_xlabel("PIT value")
    ax.set_ylabel("Density")
    ax.set_title(f"PIT Histogram (KS={ks_stat:.3f}, p={p_value:.3f})")
    ax.legend()
    plt.tight_layout()
    return pit_values, ks_stat, p_value

Dos advertencias sobre el examen KS. Asume valores PIT independientes y los pronósticos multihorizontes superpuestos están fuertemente autocorrelacionados; por lo tanto, trate el valor p como una señal aproximada y la forma del histograma como la evidencia real. Y con suficientes observaciones, la prueba rechaza la uniformidad por una mala calibración demasiado pequeña para importar; el tamaño del efecto es lo que debería impulsar una decisión.

Verificación de cobertura

Un diagnóstico más aproximado pero más directamente procesable: hacer lo $\alpha$% intervalos contienen $\alpha$% de resultados?

import numpy as np
import pandas as pd

def coverage_table(forecasts, actuals, levels=(0.50, 0.80, 0.90, 0.95)):
    """Empirical coverage at multiple nominal levels."""
    results = {}
    for level in levels:
        lower_q = (1 - level) / 2
        upper_q = 1 - lower_q
        covered = total = 0
        for forecast, actual in zip(forecasts, actuals):
            samples = forecast.samples
            h = forecast.prediction_length
            for t in range(h):
                obs = actual.values[-h + t]
                lo = np.quantile(samples[:, t], lower_q)
                hi = np.quantile(samples[:, t], upper_q)
                covered += int(lo <= obs <= hi)
                total += 1
        empirical = covered / total
        results[f"{int(level*100)}% interval"] = {
            "nominal": level, "empirical": empirical,
            "gap": empirical - level,
        }
    return pd.DataFrame(results).T

La regla general: las brechas inferiores a 2 puntos porcentuales son ruido en tamaños de muestra típicos, y las brechas superiores a 3-5 puntos son un problema de calibración real que se manifestará en el dimensionamiento de la posición. Ejecútelo por horizonte, no agrupado: la cobertura casi siempre se degrada a medida que se extiende el horizonte y la agrupación la oculta.

Cuando falla la calibración

Tres reparaciones estándar, en orden creciente de garantía. Escalado de temperatura divide el parámetro de escala previsto por un valor aprendido $T$ ajustado a datos retenidos: un parámetro, trivialmente barato, corrige el exceso o la falta de confianza uniforme, pero nada que dependa de la forma. La recalibración isotónica asigna niveles de cuantiles predichos a frecuencias observadas de forma monótona, lo que maneja la distorsión de la forma pero necesita un conjunto de calibración razonablemente grande. Predicción conforme envuelve cualquier modelo y ofrece garantías de cobertura de muestras finitas; el algoritmo completo de división conforme, el rango exacto de la estadística de orden y las trampas de interpolación y sujeción que invita un resumen de una línea se encuentran en predicción conforme para el comercio.

Consideraciones prácticas

No estacionariedad. La calibración varía a medida que cambian los regímenes de volatilidad, por lo que un conjunto de calibración fijo decae. El mecanismo diseñado exactamente para esto es la Inferencia Conforme Adaptativa, que actualiza un nivel de cobertura errónea en línea y ofrece una garantía de cobertura a largo plazo incluso en secuencias adversas; consulte la sección ACI de predicción conforme para el comercio.

Costo computacional. DeepAR con 500 rutas de muestra cuesta aproximadamente $500\times$ un solo pase hacia adelante. Para el trabajo intradiario, prefiera la regresión por cuantiles (todos los cuantiles en una sola pasada) o una cabeza paramétrica (predecir $\mu, \sigma$ una vez); reservar muestreo autorregresivo para horizontes de 4 h y diarios.

Cambios de régimen. Empareje el pronosticador con un detector de régimen y mantenga los parámetros de calibración por régimen: detección de régimen con HMM tiene la implementación y el backtest trabajados.

Previsiones multivariadas. Las distribuciones marginales por activo no son suficientes para el riesgo de cartera; la cola conjunta es lo que importa, y modelos de cópula para riesgo conjunto cuantifica hasta qué punto los supuestos de independencia lo subestiman.

Consumidores intermedios. El VaR y el déficit esperado surgen directamente de un pronóstico basado en muestras como un cuantil y una media de cola condicional: las definiciones y la receta de Monte Carlo se encuentran en modelos de cópula para riesgo conjunto. El tamaño de Kelly es lo único que no es gratuito: derivar una fracción de Kelly a partir de un intervalo de pronóstico requiere una suposición adicional sobre la distribución dentro de ese intervalo, y el artículo conforme argumenta explícitamente en contra de incluir una relación de intervalo en $f^*$. Consulte el criterio de Kelly para estrategias para conocer lo que realmente requiere la fracción.

Conclusión

La pila de evaluación para un pronóstico probabilístico es corta y no negociable: CRPS para la puntuación, porque es adecuada y reside en las unidades de su objetivo; el histograma PIT para el diagnóstico, porque nombra el modo de falla en lugar de simplemente marcar uno; cobertura por horizonte para la decisión, porque ese es el número que se relaciona con el tamaño de la posición. Todo lo que se informa como "pérdida cuantil" es CRPS hasta un factor de 2, por lo que aquí hay una métrica, no dos.

Lo incómodo es que todo esto es maquinaria para descubrir que un modelo es peor de lo que esperabas. Ese es el punto. Un modelo que sinceramente amplía sus intervalos cuando no sabe es estrictamente más útil que uno que se mantiene confiadamente estrecho, y la única forma de diferenciarlos es puntuarlos adecuadamente.

Referencias

• Salinas, D., Flunkert, V., Gasthaus, J. y Januschowski, T. (2020). "DeepAR: pronóstico probabilístico con redes recurrentes autorregresivas". Revista Internacional de Previsión, 36(3), 1181-1191.
• Koenker, R. y Bassett, G. (1978). "Cuantiles de regresión". Econométrica, 46(1), 33-50.
• Gneiting, T. y Raftery, AE (2007). "Reglas de puntuación, predicción y estimación estrictamente adecuadas". Revista de la Asociación Estadounidense de Estadística, 102(477), 359-378.
• Gneiting, T., Balabdaoui, F. y Raftery, AE (2007). "Predicciones probabilísticas, calibración y nitidez". Revista de la Royal Statistical Society: Serie B, 69(2), 243-268.
• Gal, Y. y Ghahramani, Z. (2016). "El abandono como aproximación bayesiana: representación de la incertidumbre del modelo en el aprendizaje profundo". ICML.
• Lakshminarayanan, B., Pritzel, A. y Blundell, C. (2017). "Estimación de la incertidumbre predictiva simple y escalable utilizando conjuntos profundos". NeurIPS.
• Rasul, K., Sheikh, A.-S., Schuster, I., Bergmann, U. y Vollgraf, R. (2021). "Pronóstico de series temporales probabilísticas multivariadas mediante flujos de normalización condicionados". ICLR.
• GluonTS: Modelado probabilístico de series temporales en Python. https://ts.gluon.ai