確率予測のスコアリング：CRPS、PIT キャリブレーション、DeepAR

このブログでは、点ではなく分布を予測する意義をすでに論じ、分布が得られた後にその区間をどう使うかも扱ってきた。リスクを考慮したポジションサイジングのためのコンフォーマル予測 は分布に依存しない区間と、それを使うサイジングルールを導出し、Temporal Fusion Transformer の記事は複数分位点の出力層を実装している。どちらも扱っていないのは、それらを信頼する価値があるかを決める部分、つまり 予測分布をどうスコアリングし、示された不確実性が正直かをどう確認するか である。

この記事のテーマはそこにある。具体的には次の三つだ。

1. CRPS — 確率予測の適正スコアリングルールであり、TFT の記事ですでに報告しているピンボール損失との正確な関係。
2. PIT ヒストグラム — モデルが誤キャリブレーションかどうかだけでなく、どのように誤っているかを読み取る診断。
3. DeepAR — 自己回帰サンプリングによるモデル群。このブログではベンチマークの脚注に登場するだけで、説明したことがない。

まず、必要な仕組みを決める前提を一つ。点予測の失敗はレジームによって異なる。低ボラティリティのトレンドでは条件付き分布が狭くほぼ対称なので、点予測は十分な要約になる。予定イベントの前は二峰性になり、条件付き平均は価格が 最も 到達しそうにない場所にちょうど位置する。危機局面では左裾が支配的になり、平均は下落リスクを大きく過小評価する。完全な分布を取り戻す道は三つある——パラメトリック（仮定した分布族のパラメータを予測：高速だがミススペシフィケーションのリスク）、分位点ベース（固定グリッドを予測：仮定不要だが離散的）、サンプルベース（自己回帰サンプリング、ドロップアウト、アンサンブルからモンテカルロ経路を生成：柔軟だが高コスト）。後二者が金融で優勢なのは、まさにこれらのレジーム間で分布の形が変わるからだ。

分位点予測の概要

ピンボール損失と具体的な分位点グリッドは 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 et al., 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 で十分だが、サンプルが少ないと推定 CRPS が低く偏るため、オフライン評価では 500〜1000 を使う。

予測分布を得るその他の方法

知っておく価値はあるが、独立した節を必要としない代替手法が三つある。MC Dropout（Gal and Ghahramani, 2016）は推論時もドロップアウトを有効にし、$S$ 回の順伝播にわたる平均と分散を取る。これは姿を変えた近似変分推論で、すでに学習済みのモデルに不確実性を追加する最速の方法だ。ディープアンサンブル（Lakshminarayanan et al., 2017）は異なるシードで $M$ 個のコピーを学習し、予測分布を混合分布として扱う。ベンチマークでは一貫して強いがコストは $M\times$ で、レイテンシー重視の期間には不向きでも 4 時間足や日足には使える。正規化フロー（Rasul et al., 2021）は単純な基底密度から任意の対象分布への可逆写像を学習し、分布族を選ばずに多峰性と非対称性を捉える。三つともサンプルを生成するため、次節の内容はそのまま適用できる。

CRPS：予測分布に適した指標

MSE と MAE は点予測を評価する。しかし 分布 がよかったかは分からない。見ているのが分布の一つの要約だけだからだ。標準的な代替が 連続ランク確率スコア（CRPS） であり、この記事で最も有用な指標である。

CRPS は、予測 CDF と、実際に起きた値にすべての質量を置く縮退 CDF との積分二乗距離である。

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

CRPS がその地位を得る理由は三つある。

• 適正スコアリングルールである。（Gneiting and Raftery, 2007）期待値が最小になるのは、予測分布が真の分布と一致するときだけだ。意図的な過信も、人工的に広い分布でのヘッジもスコアを改善しない。中央値に対する MSE にはこの性質がないため、CRPS は任意ではない。
• MAE を一般化する。 縮退した点予測では絶対誤差に縮退するので、CRPS は ターゲットの単位 を持つ。時間足の対数リターンを予測するなら、CRPS 0.004 は、健全性を確認できない無次元の数ではなく、平均絶対誤差 40 bps と直接比較できる。
• キャリブレーションを条件にシャープネスを評価する。 同じようにキャリブレーションされた二つの予測なら、狭い方が高いスコアになる。CRPS だけでは不十分なのもこのためだ。悪いスコアだけでは、二つの条件のどちらが失敗したか分からない。その診断が次節の役割になる。

分位点損失への橋渡し

これはブログの他の記事に欠けている接続だ。完全な CDF ではなく分位点グリッド $\mathcal{T}$ が与えられた場合、CRPS はピンボール損失で近似できる。

$\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 は同じ平均ピンボール損失をターゲット絶対値の合計で 正規化 したもので、無次元になり、規模の異なる資産間で比較できる。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 字型：過信——裾に落ちる質量が多すぎ、分布が狭すぎる。
• 山型：過小信頼——観測値が中心付近に集まり、分布が広すぎる。
• 左歪み：モデルは体系的に過大予測している。
• 右歪み：モデルは体系的に過小予測している。

共同リスクのコピュラモデル との混同を避けるために注記しておく。そこでも確率積分変換を使うが、あちらは 周辺変換 であり、GARCH-EVT の周辺分布を擬似一様な観測値に変えてコピュラを適合する前処理だ。ここでは変換をサンプル外予測に適用し、一様性は作り出す入力ではなく 検証する結果 である。数学は同じだが、推論の方向が逆になる。

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 アルゴリズム、正確な順序統計量の順位、そして一行の要約が招く補間・クリッピングの落とし穴は、トレーディングのためのコンフォーマル予測 にすべて記載している。

実務上の考慮事項

非定常性。 ボラティリティのレジームが変わるとキャリブレーションもずれるため、固定したキャリブレーション集合は劣化する。このために設計されたのが Adaptive Conformal Inference で、オンラインのミスカバレッジ水準を更新し、敵対的な系列に対しても長期カバレッジ保証を維持する——トレーディングのためのコンフォーマル予測 の ACI 節を参照。

計算コスト。 DeepAR は 500 本のサンプル経路で、順伝播一回のおよそ $500\times$ のコストになる。日中取引では、全分位点を一度に出せる分位点回帰か、$\mu, \sigma$ を一度予測するパラメトリックヘッドを優先し、自己回帰サンプリングは 4 時間足と日足のホライズンに限定する。

レジーム変化。 予測器をレジーム検出器と組み合わせ、レジームごとのキャリブレーションパラメータを維持する——HMM によるレジーム検出 に実装とバックテストがある。

多変量予測。 資産ごとの周辺分布だけではポートフォリオリスクに不十分だ。重要なのは共同の裾であり、共同リスクのコピュラモデル が、独立性の仮定による過小評価の大きさを定量化する。

下流の利用者。 VaR と期待ショートフォールは、サンプルベース予測からそれぞれ分位点と条件付き裾平均として直接得られる——定義とモンテカルロのレシピは 共同リスクのコピュラモデル にある。ケリー・サイジングだけは無料では得られない。予測区間からケリー比率を導くには、その区間内の分布について追加の仮定が必要であり、コンフォーマル予測の記事は区間比率を $f^*$ に取り付けることに明確に反対している。比率に本当に必要なものは 戦略のポジションサイジングにおけるケリー基準 を参照。

結論

確率予測の評価スタックは短く、譲れない。スコアには CRPS——適正で、ターゲットの単位を持つから。診断には PIT ヒストグラム——単に警告するだけでなく失敗モードを示すから。意思決定には ホライズン別カバレッジ——ポジションサイズに対応する数値だから。「分位点損失」として報告されるものは係数 2 を除けば CRPS なので、ここにある指標は二つではなく一つだ。

不都合な点は、これらすべてがモデルが期待より悪いと知るための仕組みだということだ。それが目的である。分からないときに正直に区間を広げるモデルは、自信満々に狭い区間を保つモデルより明らかに有用であり、両者を見分ける唯一の方法は適切にスコアリングすることだ。

参考文献

• 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). "ドロップアウトをベイズ近似として用いる：深層学習におけるモデル不確実性の表現。" 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