Scoring Probabilistic Forecasts: CRPS, PIT Calibration, and DeepAR

This blog has already made the case for predicting distributions instead of points, and has already worked out what to do with the interval once you have it: conformal prediction for risk-aware position sizing derives distribution-free intervals and the sizing rule that consumes them, and the Temporal Fusion Transformer article ships a multi-quantile output layer. What neither covers is the part that decides whether any of it is worth trusting: how you score a predictive distribution, and how you check that its stated uncertainty is honest.

That is what this article is about. Three things specifically:

1. CRPS — the proper scoring rule for probabilistic forecasts, and its exact relationship to the pinball loss the TFT article already reports.
2. The PIT histogram — a calibration diagnostic that reads off how a model is miscalibrated, not just whether it is.
3. DeepAR — the autoregressive-sampling model family, which appears on this blog only as a benchmark footnote and is never explained.

One piece of framing first, because it determines which machinery you need. A point forecast fails differently by regime: in a low-volatility trend the conditional distribution is narrow and near-symmetric and a point forecast is a fine summary; before a scheduled event it is bimodal, and the conditional mean sits exactly where the price is least likely to land; in a crisis the left tail dominates and the mean badly understates downside. Three routes recover the full distribution — parametric (predict the parameters of an assumed family: fast, misspecification risk), quantile-based (predict a fixed grid: assumption-free, discrete), and sample-based (Monte Carlo paths from autoregressive sampling, dropout, or ensembles: flexible, expensive). The latter two dominate in finance precisely because the distribution changes shape across those regimes.

Quantile Forecasts, Briefly

The pinball loss and a concrete quantile grid are already given in the TFT quantile output layer, so I will not restate the formula. The one thing worth internalizing is its asymmetry: at $\tau = 0.5$ over- and under-prediction cost the same and the loss reduces to MAE, but at $\tau = 0.95$ under-prediction is penalized $19\times$ more than over-prediction. That ratio $\tau/(1-\tau)$ is the mechanism — it is what drags the fitted value up to a level only 5% of observations should exceed.

The practical wrinkle the published article does not mention is quantile crossing: nothing in a per-quantile loss prevents $\hat{q}_{0.25} > \hat{q}_{0.75}$. With enough data and a shared backbone it is rare, but it does happen at the extreme quantiles on thin samples, and it silently breaks any downstream CRPS or coverage computation. The cheap fix is a post-hoc sort of the predicted quantile vector; the principled one is a monotonic output parameterization (predict $q_{\min}$ plus non-negative increments).

A minimal from-scratch implementation, useful when you want quantile outputs without adopting a 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 maps to position size, and we work that mapping out in full — including the edge-ratio no-trade filter — in conformal prediction for risk-aware position sizing; the $1/\sigma$ equivalence to volatility targeting is derived in volatility targeting with GARCH forecasts.

DeepAR: Autoregressive Probabilistic Forecasting

DeepAR (Salinas et al., 2020) takes the sample-based route. Rather than predicting quantiles directly, it uses an autoregressive RNN to parameterize a likelihood at each step, then draws full predictive distributions by rolling that likelihood forward.

Architecture. At each step $t$ the network consumes the previous observation $z_{t-1}$ (scaled by a per-series factor), covariates $\mathbf{x}_t$, and optional static features. An LSTM carries the state:

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

A dense head maps $\mathbf{h}_t$ to likelihood parameters $\theta_t = \text{MLP}(\mathbf{h}_t)$ — $(\mu_t, \sigma_t)$ for Gaussian, $(\mu_t, \sigma_t, \nu_t)$ for Student-t. Training maximizes $\sum_i \sum_t \log p(z_{i,t} \mid \theta_{i,t})$ across all series. At inference you sample from $p(\cdot \mid \theta_t)$, feed the sample back as the next input, and repeat $S$ times to get $S$ trajectories.

Two consequences matter for trading. The likelihood is a modeling choice, so fat tails are something you select rather than something you hope for — Student-t for heavy tails, or a Gaussian mixture $\sum_k w_k \mathcal{N}(z_t; \mu_k(\mathbf{h}_t), \sigma_k^2(\mathbf{h}_t))$ for bimodal event behavior. Multi-step forecasts are coherent: each sample path is a plausible trajectory that preserves serial correlation, which is what you need for any horizon longer than one step. (The third usual selling point — one global model across many series beats $N$ per-asset models — is the same data-efficiency argument made in the TFT article, with per-series scaling and static covariates playing the role of the static encoders there.)

DeepAR with 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 controls the Monte Carlo path count. For live use 100-200 usually suffices; for offline evaluation use 500-1000, since CRPS estimated from too few samples is biased low.

Other Routes to a Predictive Distribution

Three alternatives are worth knowing but do not need their own sections. MC Dropout (Gal and Ghahramani, 2016) keeps dropout active at inference and takes the mean and variance across $S$ forward passes; it is approximate variational inference in disguise, and it is the fastest way to bolt uncertainty onto a model you already trained. Deep ensembles (Lakshminarayanan et al., 2017) train $M$ copies from different seeds and treat the predictive distribution as a mixture — consistently strong in benchmarks, at $M\times$ the cost, which rules it out for latency-sensitive horizons but not for 4h or daily. Normalizing flows (Rasul et al., 2021) learn an invertible map from a simple base density to an arbitrary target, capturing multimodality and asymmetry without picking a parametric family. All three produce samples, so everything in the next section applies to them unchanged.

CRPS: The Right Metric for a Predictive Distribution

MSE and MAE score point forecasts. They cannot tell you whether a distribution was good, because they only ever look at one summary of it. The standard replacement is the Continuous Ranked Probability Score, and it is the single most useful thing in this article.

CRPS is the integrated squared distance between the predicted CDF and the degenerate CDF that puts all mass on what actually happened:

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

Three properties earn it its place:

• It is a proper scoring rule (Gneiting and Raftery, 2007): minimized in expectation only when the predicted distribution equals the true one. Neither deliberate overconfidence nor hedging with an artificially wide distribution improves your score. MSE-on-the-median has no such property, which is why CRPS is not optional.
• It generalizes MAE. For a degenerate point forecast it collapses to absolute error, so CRPS lives in the units of the target — forecasting hourly log returns, a CRPS of 0.004 is directly comparable to a 40 bps mean absolute error, rather than a unitless number you cannot sanity-check.
• It rewards sharpness subject to calibration. Between two equally calibrated forecasts the narrower scores better. That is also why CRPS alone is insufficient: a bad score does not tell you which of the two conditions failed, which is what the next section is for.

The Bridge to Quantile Loss

This is the connection the rest of the blog is missing. Given a quantile grid $\mathcal{T}$ rather than a full CDF, CRPS is approximated by the pinball loss:

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

Read that literally: the "quantile loss" reported in the TFT article and the CRPS discussed here are the same quantity up to a factor of 2, with the approximation tightening as the grid densifies. They are not competing metrics and there is no reason to report both.

One unit caveat, because it bites. GluonTS's mean_wQuantileLoss is that same averaged pinball loss normalized by the sum of absolute target values, which makes it dimensionless and comparable across assets of different scale. CRPS proper is not normalized and stays in return units. Do not print mean_wQuantileLoss under a label that says "CRPS" — the draft version of this article did exactly that, and it is an easy way to compare two numbers that are not on the same scale.

CRPS from Samples

With Monte Carlo samples $\{y^{(s)}\}_{s=1}^{S}$ (DeepAR, ensembles, flows), use the 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')}|$

The first term rewards accuracy, the second penalizes overdispersion. Written naively the second term is $O(S^2)$ and becomes the bottleneck when you score thousands of forecasts. Sorting first collapses it: for ascending order statistics $y_{(1)} \le \dots \le y_{(S)}$, the double sum equals $2\sum_k (2k - S - 1)\, y_{(k)}$, so the whole thing is $O(S \log S)$ dominated by the 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)

Both forms should agree closely on the same predictive distribution; if they do not, suspect quantile crossing or too few samples. In production, properscoring.crps_ensemble(observation, samples) is a well-tested drop-in.

Calibration: Is the Stated Uncertainty Honest?

A good CRPS does not guarantee that the intervals mean what they claim. A model whose nominal 90% interval captures 70% of outcomes is overconfident, and in a sizing rule that reads interval width it will lever you up precisely when it should not. The TFT article states the warning and prescribes conformal prediction as the fix; what follows is how you actually detect and diagnose the failure.

The PIT Histogram

For each observation $y_t$, compute its quantile under that step's predicted CDF:

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

If the model is calibrated, $\{u_t\}$ is uniform on $[0,1]$. The diagnostic power is that the shape of the deviation names the failure mode:

• U-shaped: overconfident — too much mass lands in the tails, the distribution is too narrow.
• Hump-shaped: underconfident — observations cluster near the center, the distribution is too wide.
• Left-skewed: the model systematically overpredicts.
• Right-skewed: the model systematically underpredicts.

A note to avoid confusion with copula models for joint risk, which also uses the probability integral transform: there it is a marginal transform, a preprocessing step that converts GARCH-EVT marginals into pseudo-uniform observations so a copula can be fitted to them. Here the transform is applied to out-of-sample forecasts and the uniformity is the result being tested, not an input being manufactured. Same mathematics, opposite direction of inference.

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

Two caveats on the KS test. It assumes independent PIT values, and overlapping multi-horizon forecasts are strongly autocorrelated — so treat the p-value as a rough flag and the histogram shape as the real evidence. And with enough observations the test rejects uniformity for miscalibration too small to matter; the effect size is what should drive a decision.

Coverage Check

A coarser but more directly actionable diagnostic: do the $\alpha$% intervals contain $\alpha$% of outcomes?

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

The rule of thumb: gaps under 2 percentage points are noise at typical sample sizes, and gaps beyond 3-5 points are a real calibration problem that will show up in position sizing. Run it per horizon, not pooled — coverage almost always degrades as the horizon extends, and pooling hides it.

When Calibration Fails

Three standard repairs, in increasing order of guarantee. Temperature scaling divides the predicted scale parameter by a learned $T$ fitted on held-out data — one parameter, trivially cheap, fixes uniform over/underconfidence but nothing shape-dependent. Isotonic recalibration maps predicted quantile levels to observed frequencies monotonically, which handles shape distortion but needs a reasonably large calibration set. Conformal prediction wraps any model and delivers finite-sample coverage guarantees; the full split-conformal algorithm, the exact order-statistic rank, and the interpolation and clamping traps that a one-line summary invites are all in conformal prediction for trading.

Practical Considerations

Non-stationarity. Calibration drifts as volatility regimes change, so a fixed calibration set decays. The mechanism designed for exactly this is Adaptive Conformal Inference, which updates an online miscoverage level and carries a long-run coverage guarantee even under adversarial sequences — see the ACI section of conformal prediction for trading.

Computational cost. DeepAR at 500 sample paths costs roughly $500\times$ a single forward pass. For intraday work prefer quantile regression (all quantiles in one pass) or a parametric head (predict $\mu, \sigma$ once); reserve autoregressive sampling for 4h and daily horizons.

Regime changes. Pair the forecaster with a regime detector and keep per-regime calibration parameters — regime detection with HMM has the worked implementation and backtest.

Multivariate forecasts. Marginal distributions per asset are not enough for portfolio risk; the joint tail is what matters, and copula models for joint risk quantifies how badly independence assumptions understate it.

Downstream consumers. VaR and expected shortfall fall straight out of a sample-based forecast as a quantile and a conditional tail mean — the definitions and the Monte Carlo recipe are in copula models for joint risk. Kelly sizing is the one thing that does not come for free: deriving a Kelly fraction from a forecast interval requires an extra assumption about the distribution within that interval, and the conformal article argues explicitly against bolting an interval ratio onto $f^*$. See the Kelly criterion for strategies for what the fraction actually requires.

Conclusion

The evaluation stack for a probabilistic forecast is short and non-negotiable: CRPS for the score, because it is proper and lives in the units of your target; the PIT histogram for the diagnosis, because it names the failure mode rather than just flagging one; per-horizon coverage for the decision, because that is the number that maps to position size. Anything reported as "quantile loss" is CRPS up to a factor of 2, so there is one metric here, not two.

The uncomfortable part is that all of this is machinery for finding out that a model is worse than you hoped. That is the point. A model that honestly widens its intervals when it does not know is strictly more useful than one that stays confidently narrow, and the only way to tell them apart is to score them properly.

References

• 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