Scoring Probabilistic Forecasts: CRPS, PIT Calibration, and DeepAR

Blog ini telah menjelaskan alasan memprediksi distribusi dan bukan titik, serta apa yang dilakukan dengan interval setelah memilikinya: conformal prediction untuk position sizing yang sadar risiko menurunkan interval distribution-free dan aturan sizing yang menggunakannya, sementara artikel Temporal Fusion Transformer mengirimkan output layer multi-quantile. Yang belum dibahas adalah bagian yang menentukan apakah semuanya layak dipercaya: bagaimana menilai predictive distribution dan memeriksa apakah ketidakpastian yang dinyatakan benar.

Itulah topik artikel ini. Tiga hal secara khusus:

1. CRPS — proper scoring rule untuk forecast probabilistik, dan hubungannya yang tepat dengan pinball loss yang telah dilaporkan artikel TFT.
2. Histogram PIT — diagnostic kalibrasi yang menunjukkan bagaimana model salah kalibrasi, bukan hanya apakah ia salah.
3. DeepAR — keluarga model autoregressive-sampling, yang di blog ini hanya muncul sebagai catatan kaki benchmark dan tidak pernah dijelaskan.

Satu framing lebih dulu karena menentukan machinery yang dibutuhkan. Point forecast gagal secara berbeda menurut regime: dalam tren volatilitas rendah, conditional distribution sempit dan hampir simetris sehingga point forecast merupakan ringkasan yang baik; sebelum event terjadwal, distribusinya bimodal dan conditional mean berada tepat di tempat harga paling kecil kemungkinannya mendarat; dalam krisis, left tail mendominasi dan mean sangat meremehkan downside. Tiga jalur memulihkan distribusi penuh — parametrik (memprediksi parameter keluarga yang diasumsikan: cepat, risiko misspecification), berbasis quantile (memprediksi grid tetap: bebas asumsi, diskret), dan berbasis sample (jalur Monte Carlo dari autoregressive sampling, dropout atau ensemble: fleksibel, mahal). Dua yang terakhir dominan di finance karena bentuk distribusi memang berubah antar-regime.

Forecast Quantile secara Singkat

Pinball loss dan grid quantile konkret sudah diberikan dalam quantile output layer TFT, jadi formulanya tidak akan diulang. Hal yang perlu diingat adalah asimetrinya: pada $\tau = 0.5$, over- dan under-prediction berbiaya sama dan loss menjadi MAE, tetapi pada $\tau = 0.95$, under-prediction dikenai penalti $19\times$ lebih besar daripada over-prediction. Rasio $\tau/(1-\tau)$ adalah mekanismenya — rasio itu menarik nilai fit ke tingkat yang hanya boleh dilampaui 5% observasi.

Kerumitan praktis yang tidak disebut artikel terbitan adalah quantile crossing: tidak ada yang mencegah $\hat{q}_{0.25} > \hat{q}_{0.75}$ dalam loss per-quantile. Dengan data cukup dan backbone bersama hal ini jarang terjadi, tetapi dapat muncul pada quantile ekstrem dalam sample tipis dan diam-diam merusak CRPS atau coverage downstream. Perbaikan murah adalah post-hoc sort vektor quantile yang diprediksi; yang principled adalah monotonic output parameterization (prediksi $q_{\min}$ ditambah increment non-negatif).

Implementasi minimal dari nol, berguna ketika ingin output quantile tanpa memakai forecasting framework:

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

Interval width dipetakan ke position size, dan pemetaan itu — termasuk edge-ratio no-trade filter — diturunkan lengkap dalam conformal prediction untuk position sizing sadar risiko; kesetaraan $1/\sigma$ dengan volatility targeting diturunkan dalam volatility targeting dengan forecast GARCH.

DeepAR: Forecast Probabilistik Autoregresif

DeepAR (Salinas et al., 2020) mengambil jalur berbasis sample. Alih-alih memprediksi quantile secara langsung, ia menggunakan autoregressive RNN untuk memparameterkan likelihood pada setiap langkah, lalu menggambar predictive distribution penuh dengan menggulirkan likelihood itu ke depan.

Arsitektur. Pada setiap langkah $t$, network menerima observasi sebelumnya $z_{t-1}$ (diskalakan dengan faktor per-series), covariate $\mathbf{x}_t$, dan optional static feature. LSTM membawa state:

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

Dense head memetakan $\mathbf{h}_t$ ke likelihood parameter $\theta_t = \text{MLP}(\mathbf{h}_t)$ — $(\mu_t, \sigma_t)$ untuk Gaussian, $(\mu_t, \sigma_t, \nu_t)$ untuk Student-t. Training memaksimalkan $\sum_i \sum_t \log p(z_{i,t} \mid \theta_{i,t})$ pada semua series. Saat inference Anda sample dari $p(\cdot \mid \theta_t)$, memasukkan sample kembali sebagai input berikutnya, dan mengulang $S$ kali untuk memperoleh $S$ trajectory.

Dua konsekuensi penting untuk trading. Likelihood adalah pilihan model, jadi fat tail dipilih dan bukan diharapkan — Student-t untuk heavy tail, atau Gaussian mixture $\sum_k w_k \mathcal{N}(z_t; \mu_k(\mathbf{h}_t), \sigma_k^2(\mathbf{h}_t))$ untuk perilaku event bimodal. Forecast multi-step koheren: setiap sample path adalah trajectory masuk akal yang mempertahankan serial correlation, yang diperlukan untuk horizon lebih panjang dari satu langkah. (Selling point ketiga — satu global model pada banyak series mengalahkan $N$ model per-asset — adalah argumen data-efficiency yang sama dalam artikel TFT, dengan per-series scaling dan static covariate berperan seperti static encoder di sana.)

DeepAR dengan 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 mengontrol jumlah jalur Monte Carlo. Untuk penggunaan live, 100-200 biasanya cukup; untuk evaluasi offline gunakan 500-1000 karena CRPS yang diestimasi dari terlalu sedikit sample bias ke bawah.

Jalur Lain Menuju Predictive Distribution

Tiga alternatif layak diketahui tetapi tidak memerlukan bagian sendiri. MC Dropout (Gal dan Ghahramani, 2016) mempertahankan dropout aktif saat inference dan mengambil mean serta variance dari $S$ forward pass; ini adalah approximate variational inference terselubung dan cara tercepat menambahkan uncertainty pada model yang sudah dilatih. Deep ensembles (Lakshminarayanan et al., 2017) melatih $M$ salinan dengan seed berbeda dan memperlakukan predictive distribution sebagai mixture — kuat secara konsisten dalam benchmark, dengan biaya $M\times$, sehingga tidak cocok untuk horizon yang sensitif latency tetapi masih cocok untuk 4h atau harian. Normalizing flows (Rasul et al., 2021) mempelajari map invertible dari base density sederhana ke target arbitrer, menangkap multimodality dan asymmetry tanpa memilih keluarga parametrik. Ketiganya menghasilkan sample, jadi semua yang ada di bagian berikut berlaku tanpa perubahan.

CRPS: Metrik yang Tepat untuk Predictive Distribution

MSE dan MAE menilai point forecast. Keduanya tidak dapat memberi tahu apakah sebuah distribution bagus karena hanya melihat satu ringkasannya. Pengganti standar adalah Continuous Ranked Probability Score, dan ini adalah hal paling berguna dalam artikel ini.

CRPS adalah jarak kuadrat terintegrasi antara predicted CDF dan degenerate CDF yang menaruh seluruh massa pada apa yang benar-benar terjadi:

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

Tiga sifat membuatnya layak digunakan:

• Proper scoring rule (Gneiting dan Raftery, 2007): dalam expectation diminimalkan hanya ketika predicted distribution sama dengan true distribution. Overconfidence yang disengaja maupun hedging dengan distribution yang dibuat terlalu lebar tidak memperbaiki score. MSE-on-the-median tidak memiliki sifat ini, sehingga CRPS bukan pilihan opsional.
• Generalisasi MAE. Untuk point forecast degenerate, ia menyusut menjadi absolute error, sehingga CRPS berada dalam unit target — untuk forecast hourly log-return, CRPS 0.004 dapat langsung dibandingkan dengan mean absolute error 40 bps, bukan angka tanpa unit yang tidak dapat diperiksa kewajarannya.
• Memberi reward pada sharpness dengan syarat calibration. Di antara dua forecast yang sama-sama calibrated, forecast yang lebih sempit mendapat score lebih baik. Itu juga alasan CRPS saja tidak cukup: score buruk tidak memberi tahu kondisi mana yang gagal, dan itulah tugas bagian berikut.

Jembatan ke Quantile Loss

Inilah hubungan yang hilang dari bagian blog lainnya. Dengan grid quantile $\mathcal{T}$ alih-alih full CDF, CRPS didekati oleh pinball loss:

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

Baca secara literal: "quantile loss" yang dilaporkan artikel TFT dan CRPS di sini adalah kuantitas yang sama hingga faktor 2, dengan approximation yang makin ketat ketika grid makin padat. Keduanya bukan metric yang bersaing dan tidak ada alasan melaporkan keduanya.

Satu caveat unit karena sering menggigit. mean_wQuantileLoss GluonTS adalah pinball loss rata-rata yang sama, tetapi dinormalisasi dengan jumlah absolute target, sehingga dimensionless dan comparable antar-asset dengan skala berbeda. CRPS yang sebenarnya tidak dinormalisasi dan tetap dalam unit return. Jangan mencetak mean_wQuantileLoss dengan label "CRPS" — draft artikel ini pernah melakukan persis itu, dan itu cara mudah membandingkan dua angka dengan skala berbeda.

CRPS dari Sample

Dengan Monte Carlo sample $\{y^{(s)}\}_{s=1}^{S}$ (DeepAR, ensemble, flow), gunakan energy form:

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

Term pertama memberi reward pada accuracy, term kedua menghukum overdispersion. Ditulis naif, term kedua memiliki $O(S^2)$ dan menjadi bottleneck ketika ribuan forecast dinilai. Sort terlebih dahulu untuk mereduksinya: untuk order statistics naik $y_{(1)} \le \dots \le y_{(S)}$, double sum sama dengan $2\sum_k (2k - S - 1)\, y_{(k)}$, sehingga keseluruhan menjadi $O(S \log S)$ yang didominasi sort.

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)

Kedua bentuk seharusnya hampir sama pada predictive distribution yang sama; jika tidak, curigai quantile crossing atau sample terlalu sedikit. Dalam production, properscoring.crps_ensemble(observation, samples) adalah drop-in yang telah teruji.

Calibration: Apakah Uncertainty yang Dinyatakan Jujur?

CRPS yang baik tidak menjamin interval berarti seperti yang diklaim. Model dengan nominal interval 90% yang hanya mencakup 70% outcome terlalu percaya diri, dan dalam sizing rule yang membaca lebar interval ia akan menaikkan leverage tepat ketika seharusnya tidak. Artikel TFT menyampaikan peringatan dan meresepkan conformal prediction sebagai perbaikan; bagian berikut menjelaskan cara mendeteksi dan mendiagnosis kegagalan itu.

Histogram PIT

Untuk setiap observation $y_t$, hitung quantile-nya di bawah predicted CDF langkah tersebut:

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

Jika model calibrated, $\{u_t\}$ uniform pada $[0,1]$. Kekuatan diagnostic-nya adalah bentuk penyimpangan menamai failure mode:

• Berbentuk U: overconfident — terlalu banyak massa jatuh di tail, distribution terlalu sempit.
• Berbentuk bukit: underconfident — observation mengelompok di tengah, distribution terlalu lebar.
• Miring ke kiri: model secara sistematis overpredict.
• Miring ke kanan: model secara sistematis underpredict.

Catatan agar tidak tertukar dengan copula model untuk joint risk, yang juga menggunakan probability integral transform: di sana ia adalah marginal transform, langkah preprocessing yang mengubah marginal GARCH-EVT menjadi pseudo-uniform observation agar copula dapat di-fit. Di sini transform diterapkan pada out-of-sample forecast dan uniformity adalah hasil yang diuji, bukan input yang dibuat. Matematikanya sama, arah inferensinya berlawanan.

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

Dua caveat untuk KS test. Ia mengasumsikan PIT value independent, sementara overlapping multi-horizon forecast sangat autocorrelated — jadi perlakukan p-value sebagai rough flag dan bentuk histogram sebagai evidence sebenarnya. Dengan observation cukup banyak, test dapat menolak uniformity untuk miscalibration yang terlalu kecil untuk berarti; effect size yang harus menggerakkan keputusan.

Coverage Check

Diagnostic yang lebih kasar tetapi lebih langsung ditindaklanjuti: apakah interval $\alpha$% memuat $\alpha$% outcome?

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

Rule of thumb: gap di bawah 2 percentage point adalah noise pada sample size tipikal, sementara gap di atas 3-5 point adalah calibration problem nyata yang akan muncul dalam position sizing. Jalankan per horizon, bukan pooled — coverage hampir selalu memburuk ketika horizon memanjang dan pooling menyembunyikannya.

Saat Calibration Gagal

Tiga perbaikan standar, dengan jaminan yang makin kuat. Temperature scaling membagi predicted scale parameter dengan learned $T$ yang di-fit pada held-out data — satu parameter, sangat murah, memperbaiki over/underconfidence uniform tetapi tidak memperbaiki bentuk yang bergantung pada shape. Isotonic recalibration memetakan predicted quantile level ke observed frequency secara monoton, sehingga menangani shape distortion tetapi membutuhkan calibration set yang cukup besar. Conformal prediction membungkus model apa pun dan memberi coverage guarantee finite-sample; algoritma split-conformal lengkap, exact order-statistic rank, serta jebakan interpolation dan clamping yang diundang oleh ringkasan satu baris ada dalam conformal prediction untuk trading.

Pertimbangan Praktis

Non-stationarity. Calibration drift ketika volatility regime berubah, sehingga fixed calibration set menurun nilainya. Mekanisme yang dirancang tepat untuk ini adalah Adaptive Conformal Inference, yang memperbarui miscoverage level secara online dan mempertahankan long-run coverage guarantee bahkan pada adversarial sequence — lihat bagian ACI dari conformal prediction untuk trading.

Computational cost. DeepAR dengan 500 sample path kira-kira berbiaya $500\times$ single forward pass. Untuk pekerjaan intraday, pilih quantile regression (semua quantile dalam satu pass) atau parametric head (prediksi $\mu, \sigma$ sekali); simpan autoregressive sampling untuk horizon 4h dan harian.

Regime change. Pasangkan forecaster dengan regime detector dan simpan calibration parameter per-regime — regime detection dengan HMM memiliki implementasi dan backtest lengkap.

Forecast multivariat. Marginal distribution per asset tidak cukup untuk portfolio risk; joint tail-lah yang penting, dan copula model untuk joint risk mengkuantifikasi seberapa parah asumsi independence meremehkannya.

Downstream consumer. VaR dan expected shortfall langsung muncul dari sample-based forecast sebagai quantile dan conditional tail mean — definisi dan resep Monte Carlo ada dalam copula model untuk joint risk. Kelly sizing adalah satu hal yang tidak gratis: menurunkan Kelly fraction dari forecast interval membutuhkan asumsi tambahan tentang distribution di dalam interval itu, dan artikel conformal secara eksplisit menentang memasang interval ratio pada $f^*$. Lihat kriteria Kelly untuk strategy sizing untuk hal yang sebenarnya dibutuhkan fraction tersebut.

Kesimpulan

Evaluation stack untuk probabilistic forecast singkat dan tidak dapat ditawar: CRPS untuk score karena ia proper dan berada dalam unit target; histogram PIT untuk diagnosis karena menamai failure mode, bukan sekadar memberi flag; coverage per-horizon untuk keputusan karena itulah angka yang dipetakan ke position size. Apa pun yang dilaporkan sebagai "quantile loss" adalah CRPS hingga faktor 2, jadi hanya ada satu metric di sini, bukan dua.

Bagian yang tidak nyaman adalah semua machinery ini dapat menunjukkan bahwa model lebih buruk daripada harapan. Itulah tujuannya. Model yang secara jujur melebarkan interval ketika tidak tahu jauh lebih berguna daripada model yang tetap sempit dengan penuh percaya diri, dan satu-satunya cara membedakan keduanya adalah menilainya dengan benar.

Referensi

• Salinas, D., Flunkert, V., Gasthaus, J., & Januschowski, T. (2020). "DeepAR: Probabilistic forecasting with autoregressive recurrent networks." International Journal of Forecasting, 36(3), 1181-1191.
• Koenker, R. & Bassett, G. (1978). "Regression Quantiles." Econometrica, 46(1), 33-50.
• Gneiting, T. & Raftery, A. E. (2007). "Strictly Proper Scoring Rules, Prediction, and Estimation." Journal of the American Statistical Association, 102(477), 359-378.
• Gneiting, T., Balabdaoui, F., & Raftery, A. E. (2007). "Probabilistic forecasts, calibration and sharpness." Journal of the Royal Statistical Society: Series B, 69(2), 243-268.
• Gal, Y. & Ghahramani, Z. (2016). "Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning." ICML.
• Lakshminarayanan, B., Pritzel, A., & Blundell, C. (2017). "Simple and Scalable Predictive Uncertainty Estimation using Deep Ensembles." NeurIPS.
• Rasul, K., Sheikh, A.-S., Schuster, I., Bergmann, U., & Vollgraf, R. (2021). "Multivariate Probabilistic Time Series Forecasting via Conditioned Normalizing Flows." ICLR.
• GluonTS: Probabilistic Time Series Modeling in Python. https://ts.gluon.ai