Punteggio delle previsioni probabilistiche: CRPS, calibrazione PIT e DeepAR

Questo blog ha già sostenuto la necessità di prevedere distribuzioni anziché punti e ha già capito cosa fare con l'intervallo una volta ottenuto: previsione conforme per il dimensionamento della posizione consapevole del rischio deriva intervalli privi di distribuzione e la regola di dimensionamento che li consuma, e l'articolo Temporal Fusion Transformer fornisce un multi-quantile strato di uscita. Ciò che nessuno dei due copre è la parte che decide se vale la pena fidarsi di qualcosa: come si valuta una distribuzione predittiva e come si controlla che l'incertezza dichiarata sia onesta.

Questo è ciò di cui tratta questo articolo. Tre cose nello specifico:

1. CRPS: la regola di punteggio corretta per le previsioni probabilistiche e la sua esatta relazione con la perdita del flipper già riportata nell'articolo di TFT.
2. L'istogramma PIT: un sistema diagnostico di calibrazione che rileva come un modello è calibrato in modo errato, non solo se lo è.
3. DeepAR: la famiglia di modelli di campionamento autoregressivo, che appare su questo blog solo come nota a piè di pagina di riferimento e non viene mai spiegata.

Prima un pezzo di struttura, perché determina di quale macchinario hai bisogno. Una previsione puntuale fallisce in modo diverso a seconda del regime: in un trend a bassa volatilità la distribuzione condizionale è stretta e quasi simmetrica e una previsione puntuale è un buon riassunto; prima di un evento programmato è bimodale e la media condizionale si trova esattamente dove è meno probabile che il prezzo raggiunga; in una crisi predomina la coda sinistra e la media sottovaluta gravemente il ribasso. Tre percorsi recuperano l'intera distribuzione: parametrica (prevedere i parametri di una famiglia presunta: rapido, rischio di errata specificazione), basata sui quantili (prevedere una griglia fissa: senza ipotesi, discreta) e basata sul campione (percorsi Monte Carlo dal campionamento autoregressivo, dall'abbandono o dagli insiemi: flessibile, costoso). Gli ultimi due dominano la finanza proprio perché la distribuzione cambia forma tra questi regimi.

Previsioni quantiliche, in breve

La perdita del flipper e una griglia quantilica concreta sono già fornite nel livello di output quantile TFT, quindi non riformulerò la formula. L'unica cosa che vale la pena interiorizzare è la sua asimmetria: at $\tau = 0.5$ la sovra e la sotto previsione costano lo stesso e la perdita si riduce al MAE, ma a $\tau = 0.95$ la sottoprevisione è penalizzata $19\times$ più che una previsione eccessiva. Quel rapporto $\tau/(1-\tau)$ è il meccanismo: è ciò che trascina il valore adattato fino a un livello che solo il 5% delle osservazioni dovrebbe superare.

Il problema pratico che l'articolo pubblicato non menziona è il incrocio quantile: nulla lo impedisce in una perdita per quantile $\hat{q}_{0.25} > \hat{q}_{0.75}$. Con abbastanza dati e una dorsale condivisa è raro, ma accade ai quantili estremi su campioni sottili e interrompe silenziosamente qualsiasi CRPS a valle o calcolo di copertura. La soluzione economica è una sorta post-hoc del vettore quantile previsto; quella di principio è una parametrizzazione dell'output monotona (predict $q_{\min}$ più incrementi non negativi).

Un'implementazione minima da zero, utile quando si desiderano risultati quantili senza adottare un quadro di previsione:

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

L'ampiezza dell'intervallo si associa alla dimensione della posizione e noi elaboriamo tale mappatura per intero, incluso il filtro no-trade del rapporto edge, in previsione conforme per dimensionamento della posizione consapevole del rischio; IL $1/\sigma$ l'equivalenza al targeting della volatilità è derivata in targeting della volatilità con previsioni GARCH.

DeepAR: previsione probabilistica autoregressiva

DeepAR (Salinas et al., 2020) segue il percorso basato su campioni. Invece di prevedere direttamente i quantili, utilizza un RNN autoregressivo per parametrizzare una probabilità in ogni passaggio, quindi disegna distribuzioni predittive complete facendo avanzare tale probabilità.

Architettura. Ad ogni passo $t$ la rete consuma l'osservazione precedente $z_{t-1}$ (ridimensionato in base a un fattore per serie), covariate $\mathbf{x}_t$e funzionalità statiche opzionali. Un LSTM trasporta lo stato:

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

Una testa densa di mappe $\mathbf{h}_t$ ai parametri di verosimiglianza $\theta_t = \text{MLP}(\mathbf{h}_t)$ — $(\mu_t, \sigma_t)$ per gaussiano, $(\mu_t, \sigma_t, \nu_t)$ per Studente-t. La formazione massimizza $\sum_i \sum_t \log p(z_{i,t} \mid \theta_{i,t})$ in tutte le serie. Per deduzione campione da $p(\cdot \mid \theta_t)$, reinserire il campione come input successivo e ripetere $S$ volte per ottenere $S$ traiettorie.

Due conseguenze sono importanti per il trading. La probabilità è una scelta di modello, quindi le code spesse sono qualcosa che selezioni piuttosto che qualcosa che speri: Student-t per code pesanti o una miscela gaussiana $\sum_k w_k \mathcal{N}(z_t; \mu_k(\mathbf{h}_t), \sigma_k^2(\mathbf{h}_t))$ per il comportamento degli eventi bimodali. Le previsioni a più passaggi sono coerenti: ciascun percorso campione è una traiettoria plausibile che preserva la correlazione seriale, che è ciò di cui hai bisogno per qualsiasi orizzonte più lungo di un passaggio. (Il terzo consueto punto di forza: un modello globale per molte serie $N$ modelli per asset: è lo stesso argomento sull'efficienza dei dati avanzato nell'articolo TFT, con il ridimensionamento per serie e le covariate statiche che svolgono il ruolo dei codificatori statici.)

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 controlla il conteggio dei percorsi Monte Carlo. Per l'uso dal vivo solitamente sono sufficienti 100-200; per la valutazione offline utilizzare 500-1000, poiché il CRPS stimato da un numero troppo basso di campioni è distorto da un valore basso.

Altri percorsi verso una distribuzione predittiva

Vale la pena conoscere tre alternative, che però non necessitano di sezioni a parte. MC Dropout (Gal e Ghahramani, 2016) mantiene attivo il dropout durante l'inferenza e prende la media e la varianza attraverso $S$ passaggi in avanti; è un'inferenza variazionale approssimativa sotto mentite spoglie ed è il modo più veloce per fissare l'incertezza su un modello che hai già addestrato. Insiemi profondi (Lakshminarayanan et al., 2017) si allenano $M$ copie da semi diversi e trattano la distribuzione predittiva come una miscela - costantemente forte nei benchmark, a $M\times$ il costo, che lo esclude per orizzonti sensibili alla latenza ma non per 4 ore o giornalieri. I flussi normalizzanti (Rasul et al., 2021) apprendono una mappa invertibile da una semplice densità di base a un obiettivo arbitrario, catturando multimodalità e asimmetria senza scegliere una famiglia parametrica. Tutti e tre producono campioni, quindi tutto nella sezione successiva si applica a loro senza modifiche.

CRPS: la giusta metrica per una distribuzione predittiva

Previsioni punteggio MSE e MAE. Non possono dirti se una distribuzione era buona, perché ne guardano solo un riassunto. La sostituzione standard è il Punteggio di probabilità classificato continuo ed è l'elemento più utile in questo articolo.

Il CRPS è la distanza quadrata integrata tra la CDF prevista e la CDF degenere che attribuisce tutta la massa a ciò che è realmente accaduto:

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

Tre proprietà gli valgono il suo posto:

• È una regola di punteggio corretta (Gneiting e Raftery, 2007): minimizzata nell'aspettativa solo quando la distribuzione prevista è uguale a quella vera. Né l’eccessiva sicurezza intenzionale né la copertura con una distribuzione artificialmente ampia migliorano il tuo punteggio. MSE-on-the-median non ha tale proprietà, motivo per cui la CRPS non è facoltativa.
• Generalizza MAE. Per una previsione puntuale degenerata crolla in errore assoluto, quindi CRPS vive nelle unità del target: prevedendo rendimenti logaritmici orari, un CRPS di 0,004 è direttamente paragonabile a un errore assoluto medio di 40 bps, piuttosto che a un numero senza unità che non è possibile controllare.
• Premia la nitidezza soggetta a calibrazione. Tra due previsioni ugualmente calibrate, quella più ristretta ottiene risultati migliori. Questo è anche il motivo per cui il CRPS da solo non è sufficiente: un punteggio negativo non dice quale delle due condizioni è fallita, a cui serve la sezione successiva.

Il ponte verso la perdita quantilica

Questa è la connessione che manca al resto del blog. Data una griglia quantilica $\mathcal{T}$ piuttosto che un CDF completo, il CRPS è approssimato dalla perdita del flipper:

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

Leggilo letteralmente: la "perdita quantile" riportata nell'articolo del TFT e il CRPS qui discusso sono la stessa quantità fino a un fattore 2, con l'approssimazione che si restringe man mano che la griglia si densifica. Non sono parametri concorrenti e non vi è motivo di segnalarli entrambi.

Un avvertimento sull'unità, perché morde. GluonTS mean_wQuantileLoss è la stessa perdita media del flipper normalizzata dalla somma dei valori target assoluti, che la rende adimensionale e comparabile tra asset di diversa scala. Il CRPS vero e proprio non è normalizzato e rimane nelle unità di ritorno. Non stampare mean_wQuantileLoss sotto un'etichetta che dice "CRPS" - la bozza di questo articolo ha fatto esattamente questo, ed è un modo semplice per confrontare due numeri che non sono sulla stessa scala.

CRPS dai campioni

Con i campioni di Monte Carlo $\{y^{(s)}\}_{s=1}^{S}$ (DeepAR, insiemi, flussi), utilizza la forma energetica:

$\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')}|$

Il primo termine premia la precisione, il secondo penalizza la sovradispersione. Scritto ingenuamente il secondo termine lo è $O(S^2)$ e diventa il collo di bottiglia quando ottieni migliaia di previsioni. L'ordinamento prima lo comprime: per le statistiche in ordine crescente $y_{(1)} \le \dots \le y_{(S)}$, la doppia somma è uguale $2\sum_k (2k - S - 1)\, y_{(k)}$, quindi il tutto è $O(S \log S)$ dominato dal genere.

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)

Entrambe le forme dovrebbero concordare strettamente sulla stessa distribuzione predittiva; in caso contrario, sospettare l'incrocio dei quantili o un numero insufficiente di campioni. Nella produzione, properscoring.crps_ensemble(observation, samples) è un drop-in ben collaudato.

Calibrazione: l'incertezza dichiarata è onesta?

Un buon CRPS non garantisce che gli intervalli significhino ciò che affermano. Un modello il cui intervallo nominale del 90% cattura il 70% dei risultati è troppo sicuro di sé e, in una regola di dimensionamento che legge l’ampiezza dell’intervallo, ti aumenterà esattamente quando non dovrebbe. L'articolo del TFT riporta l'avvertimento e prescrive la previsione conforme come soluzione; quello che segue è il modo in cui effettivamente rilevi e diagnostichi l'errore.

L'istogramma PIT

Per ogni osservazione $y_t$, calcola il suo quantile sotto la CDF prevista di quel passaggio:

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

Se il modello è calibrato, $\{u_t\}$ è uniforme $[0,1]$. Il potere diagnostico risiede nel fatto che la forma della deviazione indica la modalità di guasto:

• A forma di U: troppo sicuro di sé: troppa massa finisce nelle code, la distribuzione è troppo ristretta.
• A forma di gobba: poco sicuro di sé: le osservazioni si concentrano vicino al centro, la distribuzione è troppo ampia.
• Asimmetria a sinistra: il modello prevede sistematicamente una previsione eccessiva.
• Asimmetrica a destra: il modello sottostima sistematicamente.

Una nota per evitare confusione con copula models for joint rischio, che utilizza anche la trasformazione integrale di probabilità: esiste una trasformazione marginale, un passaggio di preelaborazione che converte i marginali GARCH-EVT in osservazioni pseudo-uniformi in modo che una copula possa essere adattata ad essi. In questo caso la trasformazione viene applicata alle previsioni fuori campione e l'uniformità è il risultato testato, non un input prodotto. Stessa matematica, direzione di inferenza opposta.

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

Due avvertenze sul test KS. Presuppone valori PIT indipendenti e le previsioni multi-orizzonte sovrapposte sono fortemente autocorrelate, quindi considera il valore p come un indicatore approssimativo e la forma dell'istogramma come la prova reale. E con un numero sufficiente di osservazioni il test rifiuta l’uniformità per un’errata calibrazione troppo piccola per avere importanza; la dimensione dell’effetto è ciò che dovrebbe guidare una decisione.

Verifica della copertura

Una diagnosi più grossolana ma più direttamente attuabile: fai il $\alpha$Gli intervalli % contengono $\alpha$% dei risultati?

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 regola empirica: i gap inferiori a 2 punti percentuali rappresentano rumore in dimensioni campione tipiche, mentre i gap superiori a 3-5 punti sono un vero problema di calibrazione che verrà visualizzato nel dimensionamento della posizione. Eseguilo per orizzonte, non in pool: la copertura quasi sempre diminuisce man mano che l'orizzonte si estende e il pooling la nasconde.

Quando la calibrazione non riesce

Tre riparazioni standard, in ordine crescente di garanzia. Ridimensionamento della temperatura divide il parametro di scala previsto per un valore appreso $T$ adattato su dati trattenuti: un parametro, banalmente economico, fissa un livello di sovra/sottoconfidenza uniforme ma non dipende dalla forma. La ricalibrazione isotonica mappa i livelli quantili previsti alle frequenze osservate in modo monotono, il che gestisce la distorsione della forma ma richiede un set di calibrazione ragionevolmente ampio. Previsione conforme avvolge qualsiasi modello e fornisce garanzie di copertura a campione finito; l'algoritmo completamente diviso-conforme, l'esatto rango statistico dell'ordine e le trappole di interpolazione e bloccaggio invitate da un riepilogo di una riga sono tutti in previsione conforme per il trading.

Considerazioni pratiche

Non stazionarietà. La calibrazione si sposta al variare dei regimi di volatilità, quindi un set di calibrazione fisso decade. Il meccanismo progettato esattamente per questo è Adaptive Conformal Inference, che aggiorna un livello di errore di copertura online e porta una garanzia di copertura a lungo termine anche in sequenze contraddittorie - vedere la sezione ACI di previsione conforme per il trading.

Costo computazionale. DeepAR con 500 percorsi campione costa all'incirca $500\times$ un unico passaggio in avanti. Per il lavoro intraday preferire la regressione quantile (tutti i quantili in un unico passaggio) o una testa parametrica (predict $\mu, \sigma$ una volta); riservare il campionamento autoregressivo per orizzonti di 4 ore e giornalieri.

Modifiche di regime. Associa il previsore a un rilevatore di regime e mantieni i parametri di calibrazione per regime: rilevamento del regime con HMM ha l'implementazione e il backtest funzionanti.

Previsioni multivariate. Le distribuzioni marginali per asset non sono sufficienti per il rischio di portafoglio; la coda congiunta è ciò che conta e modelli copula per rischio congiunto quantifica quanto i presupposti di indipendenza lo sottostimino.

Consumatori a valle. Il VaR e il deficit atteso emergono direttamente da una previsione basata su campioni come quantile e media della coda condizionale: le definizioni e la ricetta Monte Carlo si trovano in modelli copula per rischio congiunto. Il dimensionamento di Kelly è l'unica cosa che non è gratuita: derivare una frazione di Kelly da un intervallo di previsione richiede un'ipotesi aggiuntiva sulla distribuzione all'interno di quell'intervallo, e l'articolo conforme si oppone esplicitamente all'inserimento di un rapporto di intervallo su $f^*$. Vedi il criterio di Kelly per le strategie per ciò che la frazione effettivamente richiede.

Conclusione

Lo stack di valutazione per una previsione probabilistica è breve e non negoziabile: CRPS per il punteggio, perché è corretto e risiede nelle unità del tuo obiettivo; l'istogramma PIT per la diagnosi, perché nomina la modalità di guasto anziché semplicemente segnalarne una; copertura per orizzonte per la decisione, perché questo è il numero che corrisponde alla dimensione della posizione. Tutto ciò che viene riportato come "perdita quantile" è CRPS fino a un fattore 2, quindi qui c'è una metrica, non due.

La parte scomoda è che tutto questo è un meccanismo per scoprire che un modello è peggiore di quanto speravi. Questo è il punto. Un modello che allarga onestamente i propri intervalli quando non sa è strettamente più utile di uno che rimane fiduciosamente ristretto, e l’unico modo per distinguerli è assegnargli un punteggio adeguato.

Riferimenti

• Salinas, D., Flunkert, V., Gasthaus, J., & Januschowski, T. (2020). "DeepAR: previsione probabilistica con reti ricorrenti autoregressive." Giornale internazionale di previsione, 36(3), 1181-1191.
• Koenker, R. & Bassett, G. (1978). "Quantili di regressione". Econometrica, 46(1), 33-50.
• Gneiting, T. & Raftery, A. E. (2007). "Regole di punteggio, previsione e stima strettamente corrette". Giornale dell'American Statistical Association, 102(477), 359-378.
• Gneiting, T., Balabdaoui, F., & Raftery, A. E. (2007). "Previsioni probabilistiche, calibrazione e nitidezza." Giornale della Royal Statistical Society: Serie B, 69(2), 243-268.
• Gal, Y. e Ghahramani, Z. (2016). "L'abbandono come approssimazione bayesiana: rappresentazione dell'incertezza del modello nell'apprendimento profondo". ICML.
• Lakshminarayanan, B., Pritzel, A., & Blundell, C. (2017). "Stima predittiva dell'incertezza semplice e scalabile utilizzando deep ensemble." NeurIPS.
• Rasul, K., Sheikh, A.-S., Schuster, I., Bergmann, U., & Vollgraf, R. (2021). "Previsione probabilistica multivariata di serie temporali tramite flussi di normalizzazione condizionati." ICLR.
• GluonTS: modellazione probabilistica di serie temporali in Python. https://ts.gluon.ai