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 bps 的平均絕對誤差比較，而不是一個無法進行合理性檢查的無單位數字。
• 它在校準的前提下獎勵尖銳度。 在兩個校準程度相同的預測之間，較窄者分數較好。這也是為什麼 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 沒有正規化，仍以報酬單位表示。不要把 mean_wQuantileLoss 印在標示為「CRPS」的欄位下——本文草稿版正是如此，這很容易造成兩個不在同一尺度上的數字被比較。

從樣本計算 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：使用自回歸循環網路的機率預測。" International Journal of Forecasting，36(3)，1181-1191。
• Koenker, R. & Bassett, G. (1978)。"迴歸分位數。" Econometrica，46(1)，33-50。
• Gneiting, T. & Raftery, A. E. (2007)。"嚴格適當的評分規則、預測與估計。" Journal of the American Statistical Association，102(477)，359-378。
• Gneiting, T., Balabdaoui, F., & Raftery, A. E. (2007)。"機率預測、校準與尖銳度。" Journal of the Royal Statistical Society: Series 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