Scoring Probabilistic Forecasts: CRPS, PIT Calibration, and DeepAR

本博客已经说明了为什么应当预测分布而不是点值，也已经推导出得到区间后该如何使用它：面向风险的头寸规模共形预测推导了与分布无关的区间以及消费这些区间的规模规则，时间融合 Transformer一文则实现了多分位数输出层。但两篇文章都没有涉及决定这些结果是否值得信任的部分：如何为预测分布评分，以及如何检查它所声明的不确定性是否诚实。

这正是本文要讨论的内容，具体有三件事：

1. CRPS——概率预测的正确评分规则，以及它与 TFT 文章已经报告的 pinball 损失之间的精确关系。
2. PIT 直方图——一种校准诊断，可以读出模型是如何失准的，而不仅仅是是否失准。
3. DeepAR——自回归采样模型族；它在本博客中只作为基准测试脚注出现过，却从未被解释。

先说明一个决定所需工具的背景。点预测在不同状态下会以不同方式失败：在低波动趋势中，条件分布窄且近似对称，点预测是很好的概括；在预定事件之前，分布呈双峰，条件均值恰好落在价格最不可能到达的位置；在危机中，左尾占主导，均值会严重低估下行风险。有三条路径可以恢复完整分布——参数化（预测假定分布族的参数：速度快，但有设定错误风险）、基于分位数（预测固定网格：无分布假设，但结果离散）和基于样本（通过自回归采样、dropout 或集成生成蒙特卡洛路径：灵活但昂贵）。后两者在金融领域占主导地位，正是因为分布会随这些状态改变形状。

分位数预测简述

pinball 损失和一个具体的分位数网格已经在 TFT 分位数输出层中给出，因此这里不再重复公式。值得牢记的一点是它的不对称性：当 $\tau = 0.5$ 时，高估和低估的代价相同，损失退化为 MAE；但当 $\tau = 0.95$ 时，低估受到的惩罚是高估的 $19\times$。这个比率 $\tau/(1-\tau)$ 就是其中的机制——它把拟合值向上拉到只有 5% 的观测值会超过的水平。

已发表的文章没有提到的实际问题是分位数交叉：逐分位数的损失并不能阻止 $\hat{q}_{0.25} > \hat{q}_{0.75}$。数据足够多且共享骨干网络时，这种情况很少见，但在稀疏样本的极端分位数处确实会发生，并且会悄悄破坏后续任何 CRPS 或覆盖率计算。便宜的修复方法是对预测的分位数向量做事后排序；更原则性的做法是采用单调输出参数化（预测 $q_{\min}$ 加上非负增量）。

下面是一个从零实现的最小版本，适合想要分位数输出、但不想采用预测框架的场景：

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

区间宽度会映射为头寸规模；我们在面向风险的头寸规模共形预测中完整推导了这一映射，包括边际比率不交易过滤器；$1/\sigma$ 与波动率目标的等价关系则在使用 GARCH 预测进行波动率目标控制中推导。

DeepAR：自回归概率预测

DeepAR（Salinas 等，2020）采用基于样本的路径。它不是直接预测分位数，而是使用自回归 RNN 在每一步参数化一个似然，然后将该似然向前滚动以抽取完整的预测分布。

架构。 在每个时刻 $t$，网络接收前一个观测 $z_{t-1}$（按序列分别缩放）、协变量 $\mathbf{x}_t$ 以及可选的静态特征。LSTM 携带状态：

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

一个全连接头将 $\mathbf{h}_t$ 映射为似然参数 $\theta_t = \text{MLP}(\mathbf{h}_t)$——高斯分布对应 $(\mu_t, \sigma_t)$，Student-t 分布对应 $(\mu_t, \sigma_t, \nu_t)$。训练时在所有序列上最大化 $\sum_i \sum_t \log p(z_{i,t} \mid \theta_{i,t})$。推理时从 $p(\cdot \mid \theta_t)$ 中采样，将样本作为下一步输入反馈回去，并重复 $S$ 次得到 $S$ 条轨迹。

对交易而言有两个重要后果。似然是建模选择，因此肥尾是需要主动选择的东西，而不是祈求出现的东西——重尾使用 Student-t，双峰事件行为则使用高斯混合 $\sum_k w_k \mathcal{N}(z_t; \mu_k(\mathbf{h}_t), \sigma_k^2(\mathbf{h}_t))$。多步预测具有一致性：每条样本路径都是保留序列相关性的合理轨迹，这正是任何超过一步的预测范围所需要的。（第三个常见卖点——一个跨许多序列的全局模型胜过 $N$ 个按资产分别训练的模型——与 TFT 文章提出的数据效率论点相同，只是每序列缩放和静态协变量在这里扮演了静态编码器的角色。）

使用 GluonTS 的 DeepAR

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 控制蒙特卡洛路径的数量。在线使用时通常 100-200 个就足够；离线评估时使用 500-1000 个，因为用太少样本估计的 CRPS 会有向下偏差。

预测分布的其他路径

有三种替代方案值得了解，但不需要为它们各自设立章节。MC Dropout（Gal 和 Ghahramani，2016）在推理时保持 dropout 开启，并计算 $S$ 次前向传播的均值和方差；它本质上是伪装起来的近似变分推断，也是给已经训练好的模型快速添加不确定性的方式。深度集成（Lakshminarayanan 等，2017）用不同随机种子训练 $M$ 个副本，并将预测分布视为混合分布——基准测试中表现始终很强，但成本为 $M\times$，因此不适合对延迟敏感的预测范围，却适用于 4 小时或日线范围。归一化流（Rasul 等，2021）学习从简单基分布到任意目标分布的可逆映射，无需选择参数分布族即可捕捉多峰性和不对称性。三者都产生样本，因此下一节的所有内容都可以原样应用。

CRPS：预测分布的正确指标

MSE 和 MAE 用来评价点预测。它们无法告诉你一个分布是否良好，因为它们始终只观察分布的一个概括。标准替代指标是连续秩概率评分（Continuous Ranked Probability Score），也是本文最有用的单个工具。

CRPS 是预测 CDF 与退化 CDF 之间的积分平方距离，后者将全部质量放在实际发生的结果上：

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

它凭借三个性质获得了自己的位置：

• 它是正确的评分规则（Gneiting 和 Raftery，2007）：只有当预测分布等于真实分布时，期望评分才最小。故意过度自信或用人为宽大的分布进行保守对冲，都不会提高分数。以中位数计算的 MSE 没有这一性质，因此 CRPS 并非可选项。
• 它概括了 MAE。 对退化的点预测，它会退化为绝对误差，因此 CRPS 使用目标的单位——如果预测每小时对数收益，0.004 的 CRPS 可以直接与 40 个基点的平均绝对误差比较，而不是一个无法进行合理性检查的无量纲数字。
• 它在校准的约束下奖励尖锐度。 在两个校准程度相同的预测之间，更窄的那个得分更好。这也说明单独使用 CRPS 不够：糟糕的分数无法告诉你两个条件究竟是哪一个失败了，这正是下一节要解决的问题。

与分位数损失的桥梁

这是博客其他部分缺少的联系。给定分位数网格 $\mathcal{T}$ 而不是完整 CDF，CRPS 可以用 pinball 损失近似：

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

直译这一关系：TFT 文章报告的“分位数损失”和这里讨论的 CRPS 是同一个量，只差一个 2 的因子；网格越密，近似越精确。它们不是相互竞争的指标，没有理由同时报告两者。

这里有一个会造成问题的单位注意事项。GluonTS 的 mean_wQuantileLoss 是同一个平均 pinball 损失，但又除以目标绝对值之和进行归一化，所以它是无量纲的，可在不同尺度的资产之间比较。标准 CRPS 不做归一化，保持收益单位。不要在标签写着“CRPS”的地方打印 mean_wQuantileLoss——本文的草稿版本正是这样做的，这是比较两个不在同一尺度上的数字的简单方法。

从样本计算 CRPS

对于蒙特卡洛样本 $\{y^{(s)}\}_{s=1}^{S}$（DeepAR、集成模型和流模型），使用能量形式：

$\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(S^2)$，当评估数千个预测时会成为瓶颈。先排序可以将其化简：对于升序次序统计量 $y_{(1)} \le \dots \le y_{(S)}$，双重求和等于 $2\sum_k (2k - S - 1)\, y_{(k)}$，因此整体复杂度为由排序主导的 $O(S \log S)$。

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)

对于同一预测分布，两种形式应当高度一致；如果不一致，应怀疑分位数交叉或样本过少。在生产环境中，properscoring.crps_ensemble(observation, samples) 是经过充分测试的即插即用实现。

校准：声明的不确定性诚实吗？

良好的 CRPS 不能保证区间的含义与其声明相符。一个名义 90% 区间只覆盖 70% 结果的模型过于自信；在读取区间宽度的规模规则中，它会恰好在不该加杠杆时让你加大杠杆。TFT 文章提出了这一警告，并规定用共形预测修复；下面说明如何实际发现和诊断这种失败。

PIT 直方图

对每个观测值 $y_t$，计算它在该时刻预测 CDF 下的分位数：

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

如果模型经过校准，$\{u_t\}$ 在 $[0,1]$ 上服从均匀分布。这个诊断的力量在于，偏差的形状会直接指出失败模式：

• U 形： 过度自信——太多质量落在尾部，分布过窄。
• 驼峰形： 信心不足——观测值聚集在中心附近，分布过宽。
• 左偏： 模型系统性高估。
• 右偏： 模型系统性低估。

需要避免与联合风险的 copula 模型混淆：该模型也使用概率积分变换，但那里是边际变换，即把 GARCH-EVT 边际转换为伪均匀观测值的预处理步骤，以便对其拟合 copula。这里的变换应用于样本外预测，均匀性是正在检验的结果，而不是被制造出来的输入。数学相同，但推断方向相反。

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

关于 KS 检验有两个注意事项。它假设 PIT 值相互独立，而重叠的多预测范围预测具有很强的自相关，因此应把 p 值视为粗略提示，把直方图形状视为真正证据。并且，观测值足够多时，检验会拒绝那些小到无关紧要的失准；真正应当驱动决策的是效应量。

覆盖率检查

一个更粗略、但更直接可操作的诊断是：$\alpha$% 的区间是否包含 $\alpha$% 的结果？

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

经验法则是：在典型样本量下，小于 2 个百分点的差距属于噪声，超过 3-5 个百分点的差距则是真实的校准问题，并会在头寸规模中体现出来。应按预测范围分别运行，而不要合并——覆盖率几乎总会随预测范围延长而下降，合并会把这一点隐藏起来。

校准失败时

这里有三种标准修复，保证程度依次递增。温度缩放将预测的尺度参数除以在留出数据上拟合的学习参数 $T$——只需一个参数，成本极低，可以修复均匀的过度自信或信心不足，但无法处理依赖形状的问题。保序重新校准将预测分位数水平单调映射到观测频率，可以处理形状失真，但需要相当大的校准集。共形预测包装任意模型并提供有限样本覆盖保证；完整的 split-conformal 算法、精确的次序统计量排名，以及单行概述容易掩盖的插值和截断陷阱，都在面向交易的共形预测中。

实际注意事项

非平稳性。 校准会随波动率状态变化而漂移，因此固定的校准集会逐渐失效。专门为此设计的机制是自适应共形推断，它更新在线未覆盖率，即使面对对抗序列也能提供长期覆盖保证——参见面向交易的共形预测中的 ACI 部分。

计算成本。 DeepAR 生成 500 条样本路径的成本约为单次前向传播的 $500\times$。日内工作应优先使用分位数回归（一次传播得到所有分位数）或参数化头（一次预测 $\mu, \sigma$）；将自回归采样留给 4 小时和日线预测范围。

状态变化。 将预测器与状态检测器配对，并保留每个状态的校准参数——使用 HMM 进行状态检测提供了完整的实现和回测。

多变量预测。 每种资产的边际分布不足以应对投资组合风险；真正重要的是联合尾部，而联合风险的 copula 模型量化了独立性假设会在多大程度上低估它。

下游使用者。 VaR 和预期短缺可以直接从基于样本的预测中得到，分别对应分位数和条件尾部均值——定义和蒙特卡洛方法见联合风险的 copula 模型。Kelly 规模是唯一不会自动得到的东西：从预测区间推导 Kelly 比例，需要额外假设该区间内部的分布，而共形预测文章明确反对把区间比率直接接到 $f^*$ 上。关于该比例实际需要什么，请参见策略的 Kelly 准则。

结论

概率预测的评估工具链很短，而且不可妥协：用 CRPS 评分，因为它是正确评分规则且使用目标单位；用 PIT 直方图 诊断，因为它指出失败模式，而不只是发出失败信号；用逐预测范围覆盖率做决策，因为这个数字会映射到头寸规模。任何被报告为“分位数损失”的东西都只是 CRPS 乘以一个 2 的因子，因此这里其实只有一个指标，而不是两个。

令人不适的一点是，所有这些工具都是用来发现模型比你希望的更差。这正是它们的意义。一个在不知道时诚实扩大区间的模型，严格来说比始终自信地保持窄区间的模型更有用；区分二者的唯一方法就是正确地给它们评分。

参考文献

• Salinas, D.、Flunkert, V.、Gasthaus, J. 和 Januschowski, T.（2020）。"DeepAR：使用自回归循环网络进行概率预测。" 国际预测杂志，36(3)，1181-1191。
• Koenker, R. 和 Bassett, G.（1978）。"回归分位数。" 计量经济学，46(1)，33-50。
• Gneiting, T. 和 Raftery, A. E.（2007）。"严格正确的评分规则、预测与估计。" 美国统计学会杂志，102(477)，359-378。
• Gneiting, T.、Balabdaoui, F. 和 Raftery, A. E.（2007）。"概率预测、校准与尖锐度。" 英国皇家统计学会杂志：B 辑，69(2)，243-268。
• Gal, Y. 和 Ghahramani, Z.（2016）。"作为贝叶斯近似的 Dropout：表示深度学习中的模型不确定性。" ICML。
• Lakshminarayanan, B.、Pritzel, A. 和 Blundell, C.（2017）。"使用深度集成进行简单且可扩展的预测不确定性估计。" NeurIPS。
• Rasul, K.、Sheikh, A.-S.、Schuster, I.、Bergmann, U. 和 Vollgraf, R.（2021）。"通过条件归一化流进行多元概率时间序列预测。" ICLR。
• GluonTS：Python 中的概率时间序列建模。https://ts.gluon.ai