Pontuação de previsões probabilísticas: CRPS, calibração PIT e DeepAR

Este blog já defendeu a previsão de distribuições em vez de pontos e já descobriu o que fazer com o intervalo quando você o tiver: previsão conforme para dimensionamento de posição com reconhecimento de risco deriva intervalos livres de distribuição e a regra de dimensionamento que os consome, e o artigo Transformador de fusão temporal envia um artigo camada de saída multi-quantil. O que nenhum deles cobre é a parte que decide se vale a pena confiar em algo: como você avalia uma distribuição preditiva e como verifica se sua incerteza declarada é honesta.

É disso que trata este artigo. Três coisas especificamente:

1. CRPS — a regra de pontuação adequada para previsões probabilísticas e sua relação exata com a perda de pinball que o artigo do TFT já relata.
2. O histograma PIT — um diagnóstico de calibração que lê como um modelo está mal calibrado, e não apenas se está.
3. DeepAR — a família de modelos de amostragem autorregressiva, que aparece neste blog apenas como uma nota de rodapé de referência e nunca é explicada.

Uma peça de estrutura primeiro, porque ela determina de qual maquinário você precisa. Uma previsão pontual falha de forma diferente por regime: numa tendência de baixa volatilidade a distribuição condicional é estreita e quase simétrica e uma previsão pontual é um bom resumo; antes de um evento programado é bimodal, e a média condicional fica exatamente onde o preço tem menos probabilidade de cair; em uma crise, a cauda esquerda domina e a média subestima muito o lado negativo. Três rotas recuperam a distribuição completa — paramétrica (prever os parâmetros de uma família assumida: rápido, risco de especificação incorreta), baseada em quantil (prever uma grade fixa: livre de suposições, discreta) e baseada em amostra (caminhos de Monte Carlo de amostragem autorregressiva, abandono ou conjuntos: flexível, caro). Os dois últimos dominam nas finanças precisamente porque a distribuição muda de forma entre esses regimes.

Previsões de quantis, resumidamente

A perda de pinball e uma grade de quantil concreta já são fornecidas na camada de saída de quantil do TFT, então não vou reformular a fórmula. A única coisa que vale a pena internalizar é a sua assimetria: pelo menos $\tau = 0.5$ a previsão excessiva e insuficiente custam o mesmo e a perda reduz-se ao MAE, mas pelo menos $\tau = 0.95$ subpredição é penalizada $19\times$ mais do que uma previsão exagerada. Essa proporção $\tau/(1-\tau)$ é o mecanismo - é o que arrasta o valor ajustado até um nível que apenas 5% das observações deveriam exceder.

O problema prático que o artigo publicado não menciona é o cruzamento de quantis: nada em uma perda por quantil impede $\hat{q}_{0.25} > \hat{q}_{0.75}$. Com dados suficientes e um backbone compartilhado, isso é raro, mas acontece nos quantis extremos em amostras finas e quebra silenciosamente qualquer CRPS downstream ou cálculo de cobertura. A solução barata é uma espécie post-hoc do vetor quantil previsto; o princípio é uma parametrização de saída monotônica (prever $q_{\min}$ mais incrementos não negativos).

Uma implementação mínima do zero, útil quando você deseja resultados quantílicos sem adotar uma estrutura de previsão:

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

A largura do intervalo é mapeada para o tamanho da posição, e trabalhamos esse mapeamento por completo — incluindo o filtro sem negociação de proporção de borda — em previsão conforme para dimensionamento de posição ciente de risco; o $1/\sigma$ a equivalência à meta de volatilidade é derivada em meta de volatilidade com previsões GARCH.

DeepAR: previsão probabilística autoregressiva

DeepAR (Salinas et al., 2020) segue o caminho baseado em amostra. Em vez de prever quantis diretamente, ele usa um RNN autoregressivo para parametrizar uma probabilidade em cada etapa e, em seguida, desenha distribuições preditivas completas, avançando essa probabilidade.

Arquitetura. Em cada etapa $t$ a rede consome a observação anterior $z_{t-1}$ (escalonado por um fator por série), covariáveis $\mathbf{x}_t$e recursos estáticos opcionais. Um LSTM carrega o estado:

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

Uma cabeça densa mapeia $\mathbf{h}_t$ para parâmetros de probabilidade $\theta_t = \text{MLP}(\mathbf{h}_t)$ — $(\mu_t, \sigma_t)$ para Gaussiano, $(\mu_t, \sigma_t, \nu_t)$ para estudante-t. O treinamento maximiza $\sum_i \sum_t \log p(z_{i,t} \mid \theta_{i,t})$ em todas as séries. Na inferência você amostra de $p(\cdot \mid \theta_t)$, alimente a amostra como a próxima entrada e repita $S$ vezes para conseguir $S$ trajetórias.

Duas consequências são importantes para a negociação. A probabilidade é uma escolha de modelagem, então caudas grossas são algo que você seleciona e não algo que você espera — Student-t para caudas pesadas ou uma mistura gaussiana $\sum_k w_k \mathcal{N}(z_t; \mu_k(\mathbf{h}_t), \sigma_k^2(\mathbf{h}_t))$ para comportamento de evento bimodal. As previsões em várias etapas são coerentes: cada caminho de amostra é uma trajetória plausível que preserva a correlação serial, que é o que você precisa para qualquer horizonte maior que uma etapa. (O terceiro argumento de venda habitual – um modelo global em muitas séries $N$ modelos por ativo — é o mesmo argumento de eficiência de dados apresentado no artigo do TFT, com escalonamento por série e covariáveis estáticas desempenhando o papel dos codificadores estáticos.)

DeepAR com 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 a contagem de caminhos de Monte Carlo. Para uso ao vivo, 100-200 geralmente é suficiente; para avaliação off-line, use 500-1000, uma vez que a CRPS estimada a partir de poucas amostras tem um viés baixo.

Outras rotas para uma distribuição preditiva

Vale a pena conhecer três alternativas, mas não precisam de seções próprias. MC Dropout (Gal e Ghahramani, 2016) mantém o abandono ativo na inferência e calcula a média e a variância entre $S$ passes para frente; é uma inferência variacional aproximada disfarçada e é a maneira mais rápida de inserir a incerteza em um modelo que você já treinou. Conjuntos profundos (Lakshminarayanan et al., 2017) treinam $M$ cópias de sementes diferentes e tratam a distribuição preditiva como uma mistura - consistentemente forte em benchmarks, em $M\times$ o custo, o que exclui horizontes sensíveis à latência, mas não para 4h ou diariamente. Normalização de fluxos (Rasul et al., 2021) aprende um mapa invertível de uma densidade de base simples para um alvo arbitrário, capturando multimodalidade e assimetria sem escolher uma família paramétrica. Todos os três produzem amostras, então tudo na próxima seção se aplica a eles inalterado.

CRPS: a métrica certa para uma distribuição preditiva

Previsões de pontos de pontuação MSE e MAE. Eles não podem dizer se uma distribuição foi boa, porque eles só olham para um resumo dela. A substituição padrão é a Pontuação de Probabilidade Classificada Contínua, e é a coisa mais útil neste artigo.

CRPS é a distância quadrada integrada entre o CDF previsto e o CDF degenerado que coloca toda a massa no que realmente aconteceu:

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

Três propriedades conquistam seu lugar:

• É uma regra de pontuação adequada (Gneiting e Raftery, 2007): minimizada na expectativa apenas quando a distribuição prevista é igual à verdadeira. Nem o excesso de confiança deliberado nem a cobertura com uma distribuição artificialmente ampla melhoram a sua pontuação. O MSE-on-the-median não possui tal propriedade, razão pela qual o CRPS não é opcional.
• Ele generaliza o MAE. Para uma previsão de ponto degenerado, ele cai para o erro absoluto, então o CRPS vive nas unidades do alvo — prevendo retornos de log por hora, um CRPS de 0,004 é diretamente comparável a um erro absoluto médio de 40 bps, em vez de um número sem unidade que você não pode verificar a sanidade.
• Recompensa a nitidez sujeita a calibração. Entre duas previsões igualmente calibradas, as pontuações mais estreitas melhoram. É também por isso que o CRPS por si só é insuficiente: uma pontuação ruim não indica qual das duas condições falhou, e é para isso que serve a próxima seção.

A ponte para a perda de quantil

Esta é a conexão que falta ao resto do blog. Dada uma grade quantílica $\mathcal{T}$ em vez de um CDF completo, o CRPS é aproximado pela perda de pinball:

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

Leia isso literalmente: a "perda quantil" relatada no artigo do TFT e no CRPS discutido aqui são a mesma quantidade até um fator de 2, com a aproximação diminuindo à medida que a grade se densifica. Não são métricas concorrentes e não há razão para relatar ambas.

Uma advertência de unidade, porque morde. GluonTS mean_wQuantileLoss é a mesma perda média de pinball normalizada pela soma dos valores-alvo absolutos, o que a torna adimensional e comparável entre ativos de diferentes escalas. O CRPS propriamente dito não é normalizado e permanece em unidades de retorno. Não imprima mean_wQuantileLoss sob um rótulo que diz “CRPS” – a versão preliminar deste artigo fez exatamente isso e é uma maneira fácil de comparar dois números que não estão na mesma escala.

CRPS de amostras

Com amostras de Monte Carlo $\{y^{(s)}\}_{s=1}^{S}$ (DeepAR, conjuntos, fluxos), use a forma de energia:

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

O primeiro termo recompensa a precisão, o segundo penaliza a superdispersão. Escrito ingenuamente, o segundo termo é $O(S^2)$ e se torna o gargalo quando você obtém milhares de previsões. Classificar primeiro o recolhe: para estatísticas de ordem crescente $y_{(1)} \le \dots \le y_{(S)}$, a soma dupla é igual $2\sum_k (2k - S - 1)\, y_{(k)}$, então a coisa toda é $O(S \log S)$ dominado pelo tipo.

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 as formas devem concordar estreitamente com a mesma distribuição preditiva; caso contrário, suspeite de cruzamento de quantis ou de poucas amostras. Na produção, properscoring.crps_ensemble(observation, samples) é um drop-in bem testado.

Calibração: a incerteza declarada é honesta?

Um bom CRPS não garante que os intervalos signifiquem o que afirmam. Um modelo cujo intervalo nominal de 90% captura 70% dos resultados é excessivamente confiante e, em uma regra de dimensionamento que indica a largura do intervalo, ele o alavancará precisamente quando não deveria. O artigo do TFT declara o aviso e prescreve a previsão conforme como solução; o que se segue é como você realmente detecta e diagnostica a falha.

O histograma PIT

Para cada observação $y_t$, calcule seu quantil sob o CDF previsto para essa etapa:

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

Se o modelo estiver calibrado, $\{u_t\}$ é uniforme em $[0,1]$. O poder de diagnóstico é que a forma do desvio nomeia o modo de falha:

• Em forma de U: excesso de confiança — muita massa cai nas caudas, a distribuição é muito estreita.
• Em forma de corcunda: falta de confiança — as observações agrupam-se perto do centro, a distribuição é muito ampla.
• Inclinação à esquerda: o modelo superestima sistematicamente.
• Inclinação para a direita: o modelo subestima sistematicamente.

Uma nota para evitar confusão com modelos de cópula para risco conjunto, que também usa a transformação integral de probabilidade: aí está uma transformação marginal, uma etapa de pré-processamento que converte marginais GARCH-EVT em observações pseudouniformes para que uma cópula possa ser ajustada a elas. Aqui a transformação é aplicada a previsões fora da amostra e a uniformidade é o resultado que está sendo testado, e não um insumo sendo fabricado. Mesma matemática, direção oposta de inferência.

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

Duas advertências sobre o teste KS. Ele assume valores PIT independentes, e as previsões multihorizontes sobrepostas são fortemente autocorrelacionadas - portanto, trate o valor p como um sinalizador aproximado e o formato do histograma como a evidência real. E com observações suficientes, o teste rejeita a uniformidade devido a um erro de calibração demasiado pequeno para ter importância; o tamanho do efeito é o que deve orientar uma decisão.

Verificação de cobertura

Um diagnóstico mais grosseiro, mas mais diretamente acionável: faça o $\alpha$% intervalos contêm $\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

A regra prática: lacunas abaixo de 2 pontos percentuais são ruído em tamanhos de amostra típicos, e lacunas além de 3 a 5 pontos são um problema real de calibração que aparecerá no dimensionamento de posição. Execute-o por horizonte, não em pool – a cobertura quase sempre se degrada à medida que o horizonte se estende, e o pool a oculta.

Quando a calibração falha

Três reparos padrão, em ordem crescente de garantia. A escala de temperatura divide o parâmetro de escala previsto por um valor aprendido $T$ ajustado em dados retidos - um parâmetro, trivialmente barato, corrige excesso/falta de confiança uniforme, mas nada dependente da forma. Recalibração isotônica mapeia os níveis de quantis previstos para frequências observadas monotonicamente, o que lida com a distorção de forma, mas precisa de um conjunto de calibração razoavelmente grande. A previsão conforme envolve qualquer modelo e oferece garantias de cobertura de amostras finitas; o algoritmo split-conformal completo, a classificação exata da estatística da ordem e as armadilhas de interpolação e fixação que um resumo de uma linha convida estão todos em previsão conforme para negociação.

Considerações Práticas

Não estacionariedade. A calibração varia à medida que os regimes de volatilidade mudam, portanto, um conjunto de calibração fixo decai. O mecanismo projetado exatamente para isso é a Inferência Conforme Adaptativa, que atualiza um nível de cobertura incorreta on-line e carrega uma garantia de cobertura de longo prazo, mesmo sob sequências adversárias — consulte a seção ACI de previsão conforme para negociação.

Custo computacional. DeepAR com 500 caminhos de amostra custa aproximadamente $500\times$ um único passe para frente. Para trabalho intradiário, prefira regressão quantílica (todos os quantis em uma passagem) ou cabeça paramétrica (prever $\mu, \sigma$ uma vez); reserve amostragem autorregressiva para horizontes de 4h e diários.

Mudanças de regime. Combine o previsor com um detector de regime e mantenha os parâmetros de calibração por regime — detecção de regime com HMM tem a implementação e o backtest trabalhados.

Previsões multivariadas. As distribuições marginais por ativo não são suficientes para o risco do portfólio; a cauda conjunta é o que importa, e modelos de cópula para risco conjunto quantifica o quanto as suposições de independência a subestimam.

Consumidores downstream. O VaR e o déficit esperado saem diretamente de uma previsão baseada em amostra como um quantil e uma média de cauda condicional — as definições e a receita de Monte Carlo estão em modelos de cópula para risco conjunto. O dimensionamento de Kelly é a única coisa que não vem de graça: derivar uma fração de Kelly a partir de um intervalo de previsão requer uma suposição extra sobre a distribuição dentro desse intervalo, e o artigo conforme argumenta explicitamente contra fixar uma proporção de intervalo em $f^*$. Consulte o critério de Kelly para estratégias para saber o que a fração realmente exige.

Conclusão

A pilha de avaliação para uma previsão probabilística é curta e inegociável: CRPS para a pontuação, porque é adequada e reside nas unidades do seu alvo; o histograma PIT para o diagnóstico, porque ele nomeia o modo de falha em vez de apenas sinalizar um; cobertura por horizonte para a decisão, porque esse é o número que mapeia o tamanho da posição. Qualquer coisa relatada como “perda quantílica” é CRPS até um fator de 2, portanto há uma métrica aqui, não duas.

A parte desconfortável é que tudo isso é um mecanismo para descobrir que um modelo é pior do que você esperava. Esse é o ponto. Um modelo que alarga honestamente os seus intervalos quando não sabe é estritamente mais útil do que aquele que se mantém estreito e confiante, e a única forma de os distinguir é pontuando-os adequadamente.

Referências

• Salinas, D., Flunkert, V., Gasthaus, J., & Januschowski, T. (2020). "DeepAR: Previsão probabilística com redes recorrentes autorregressivas." Jornal Internacional de Previsão, 36(3), 1181-1191.
• Koenker, R. & Bassett, G. (1978). "Quantis de regressão." Econométrica, 46(1), 33-50.
• Gneiting, T. & Raftery, AE (2007). "Regras de pontuação, previsão e estimativa estritamente adequadas." Jornal da Associação Estatística Americana, 102(477), 359-378.
• Gneiting, T., Balabdaoui, F., & Raftery, AE (2007). "Previsões probabilísticas, calibração e nitidez." Journal of the Royal Statistical Society: Série B, 69(2), 243-268.
• Gal, Y. & Ghahramani, Z. (2016). "Abandono como aproximação bayesiana: representando a incerteza do modelo no aprendizado profundo." ICML.
• Lakshminarayanan, B., Pritzel, A. e Blundell, C. (2017). "Estimativa de incerteza preditiva simples e escalonável usando conjuntos profundos." NeurIPS.
• Rasul, K., Sheikh, A.-S., Schuster, I., Bergmann, U., & Vollgraf, R. (2021). "Previsão probabilística multivariada de séries temporais por meio de fluxos de normalização condicionada." ICLR.
• GluonTS: Modelagem Probabilística de Séries Temporais em Python. https://ts.gluon.ai