Scoring Probabilistic Forecasts: CRPS, PIT Calibration, and DeepAR

Этот блог уже обосновал, почему стоит предсказывать распределения, а не отдельные точки, и уже разобрал, что делать с интервалом после его получения: конформное предсказание для позиционирования с учетом риска выводит интервалы без предположения о распределении и правило позиционирования, которое их использует, а статья о Temporal Fusion Transformer реализует выходной слой с несколькими квантилями. Но ни одна из них не рассматривает то, что определяет, заслуживает ли результат доверия: как оценивать прогнозное распределение и как проверять честность заявленной неопределенности.

Об этом и статья. Точнее, о трех вещах:

1. CRPS — правильное правило оценки вероятностных прогнозов и точная связь с pinball loss, уже показанным в статье о TFT.
2. Гистограмма PIT — диагностический инструмент калибровки, который показывает, как модель некалибрована, а не только сам факт этого.
3. DeepAR — семейство авторегрессионных моделей с выборкой, которое в этом блоге появляется лишь в примечании к benchmark и нигде не объясняется.

Сначала уточним одну рамку, поскольку она определяет нужный инструментарий. Точечный прогноз ошибается по-разному в разных режимах: при тренде с низкой волатильностью условное распределение узкое и почти симметричное, поэтому точка хорошо его резюмирует; перед запланированным событием оно бимодально, и условное среднее оказывается именно там, где цена наименее вероятна; в кризисе доминирует левый хвост, а среднее сильно занижает риск снижения. Полное распределение можно восстановить тремя путями — параметрическим (предсказывать параметры выбранного семейства: быстро, но с риском неправильной спецификации), квантильным (предсказывать фиксированную сетку: без допущений, но дискретно) и выборочным (пути Монте-Карло из авторегрессионной выборки, dropout или ансамблей: гибко, но дорого). В финансах последние два подхода особенно важны, потому что форма распределения меняется между этими режимами.

Квантильные прогнозы: кратко

Pinball loss и конкретная сетка квантилей уже приведены в выходном квантильном слое TFT, поэтому повторять формулу не будем. Важно усвоить ее асимметрию: при $\tau = 0.5$ переоценка и недооценка стоят одинаково, и loss сводится к MAE, но при $\tau = 0.95$ недооценка штрафуется в $19\times$ сильнее переоценки. Отношение $\tau/(1-\tau)$ и есть механизм: он вытягивает обученное значение к уровню, который должны превышать только 5% наблюдений.

Практическая тонкость, о которой опубликованная статья не говорит, — пересечение квантилей: ничто в loss для отдельных квантилей не запрещает $\hat{q}_{0.25} > \hat{q}_{0.75}$. При достаточном объеме данных и общем backbone это редко, но на малых выборках случается для крайних квантилей и незаметно ломает любой последующий расчет CRPS или coverage. Дешевое исправление — отсортировать вектор предсказанных квантилей после обучения; принципиальное — параметризовать выход монотонно (предсказывать $q_{\min}$ и неотрицательные приращения).

Минимальная реализация с нуля пригодится, если нужны квантильные выходы без перехода на forecasting-фреймворк:

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

Связь ширины интервала с размером позиции, включая фильтр no-trade по отношению edge, подробно выведена в конформном предсказании для позиционирования с учетом риска; эквивалентность $1/\sigma$ таргетированию волатильности выведена в таргетировании волатильности с прогнозами GARCH.

DeepAR: вероятностное авторегрессионное прогнозирование

DeepAR (Salinas и др., 2020) использует выборочный подход. Вместо прямого предсказания квантилей он применяет авторегрессионную RNN, параметризующую likelihood на каждом шаге, а затем получает полные прогнозные распределения, прокручивая этот likelihood вперед.

Архитектура. На каждом шаге $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$ в параметры likelihood $\theta_t = \text{MLP}(\mathbf{h}_t)$ — $(\mu_t, \sigma_t)$ для Gaussian, $(\mu_t, \sigma_t, \nu_t)$ для Student-t. Обучение максимизирует $\sum_i \sum_t \log p(z_{i,t} \mid \theta_{i,t})$ по всем рядам. На инференсе выполняется выборка из $p(\cdot \mid \theta_t)$, выбранное значение подается как следующий вход, и это повторяется $S$ раз для получения $S$ траекторий.

Для торговли важны два следствия. Likelihood — это выбор модели, поэтому толстые хвосты нужно выбрать, а не надеяться на них: Student-t для тяжелых хвостов или смесь Gaussian $\sum_k w_k \mathcal{N}(z_t; \mu_k(\mathbf{h}_t), \sigma_k^2(\mathbf{h}_t))$ для бимодального поведения во время событий. Многошаговые прогнозы согласованы: каждый выборочный путь является правдоподобной траекторией и сохраняет серийную корреляцию, необходимую для любого горизонта длиннее одного шага. (Третий стандартный аргумент — одна глобальная модель для многих рядов лучше $N$ моделей по активам — это тот же аргумент эффективности данных, что в статье о TFT, где масштабирование по рядам и статические ковариаты играют роль статических энкодеров.)

DeepAR с 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 управляет количеством траекторий Монте-Карло. Для live-использования обычно достаточно 100-200; для офлайн-оценки используйте 500-1000, поскольку CRPS, оцененный по слишком малому числу выборок, будет смещен вниз.

Другие пути к прогнозному распределению

Стоит знать еще три альтернативы, но отдельные разделы им не нужны. MC Dropout (Gal и Ghahramani, 2016) оставляет dropout активным на инференсе и берет среднее и дисперсию по $S$ прямым проходам; это приближенный вариационный вывод, замаскированный под другой подход, и самый быстрый способ добавить неопределенность уже обученной модели. Глубокие ансамбли (Lakshminarayanan и др., 2017) обучают $M$ копий с разными seed и трактуют прогнозное распределение как смесь — они стабильно сильны в benchmark, но стоят $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 измеряется в единицах цели: при прогнозировании часовых лог-доходностей CRPS 0.004 напрямую сопоставим со средней абсолютной ошибкой 40 bps, а не с безразмерным числом, которое невозможно проверить на здравый смысл.
• Он вознаграждает остроту при условии калибровки. Из двух одинаково калиброванных прогнозов лучше оценивается более узкий. Именно поэтому одного CRPS недостаточно: плохая оценка не говорит, какое из двух условий нарушено, а этому посвящен следующий раздел.

Связь с квантильным loss

Это связь, которой не хватает остальному блогу. Имея сетку квантилей $\mathcal{T}$ вместо полной CDF, CRPS можно приблизить pinball loss:

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

Прочитайте это буквально: указанный в статье о TFT quantile loss и обсуждаемый здесь CRPS — одна и та же величина с множителем 2, а точность приближения растет при уплотнении сетки. Это не конкурирующие метрики, и нет причин отчитываться по обеим.

Важно уточнить единицы. mean_wQuantileLoss в GluonTS — тот же усредненный pinball loss, нормированный на сумму абсолютных значений цели, поэтому он безразмерен и сопоставим между активами разного масштаба. Сам 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)

Обе формы должны близко совпадать для одного и того же прогнозного распределения; если это не так, проверьте пересечение квантилей или слишком малое число выборок. В production можно использовать хорошо протестированную замену 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-моделями совместного риска, где тоже используется probability integral transform. Там это маргинальное преобразование, шаг предобработки, превращающий маргинальные 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-value как грубый сигнал, а форму гистограммы — как основное свидетельство. Кроме того, при большом числе наблюдений тест отвергает равномерность даже при слишком малой для практики некалибровке; решать следует по размеру эффекта.

Проверка покрытия

Более грубая, но непосредственно применимая диагностика: содержат ли интервалы $\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 пунктов — реальная проблема калибровки, которая проявится в позиционировании. Запускайте расчет отдельно по горизонтам, а не объединяйте их: покрытие почти всегда ухудшается с ростом горизонта, а объединение это скрывает.

Когда калибровка не работает

Три стандартных исправления в порядке усиления гарантий. Temperature scaling делит предсказанный scale-параметр на обученный $T$, подогнанный на отложенных данных: один параметр, почти бесплатное решение, исправляет равномерную излишнюю или недостаточную уверенность, но не зависимые от формы искажения. Изотоническая рекалибровка монотонно отображает предсказанные уровни квантилей в наблюдаемые частоты; она работает с искажением формы, но требует достаточно большой калибровочной выборки. Конформное предсказание оборачивает любую модель и дает гарантии покрытия на конечной выборке; полный алгоритм split-conformal, точный ранг порядковой статистики, а также ловушки интерполяции и clamping, которые скрывает однострочное резюме, разобраны в конформном предсказании для торговли.

Практические соображения

Нестационарность. Калибровка дрейфует при смене режимов волатильности, поэтому фиксированная калибровочная выборка устаревает. Механизм, созданный именно для этого, — Adaptive Conformal Inference: он обновляет онлайн-уровень непокрытия и сохраняет гарантию долгосрочного покрытия даже для неблагоприятных последовательностей. См. раздел ACI в конформном предсказании для торговли.

Вычислительная стоимость. DeepAR с 500 траекториями стоит примерно $500\times$ одного прямого прохода. Для внутридневной работы предпочитайте квантильную регрессию (все квантили за один проход) или параметрическую голову (однократно предсказать $\mu, \sigma$); авторегрессионную выборку оставляйте для 4-часовых и дневных горизонтов.

Смена режимов. Соединяйте прогнозист с детектором режима и храните параметры калибровки по каждому режиму — в детектировании режимов с HMM есть готовая реализация и backtest.

Многомерные прогнозы. Маргинальных распределений по активам недостаточно для риска портфеля; важен совместный хвост, а copula-модели совместного риска показывают, насколько сильно предположение о независимости его занижает.

Потребители на downstream. VaR и expected shortfall напрямую получаются из выборочного прогноза как квантиль и условное среднее хвоста — определения и рецепт Монте-Карло приведены в copula-моделях совместного риска. Kelly sizing не дается бесплатно: вывод доли Kelly из прогнозного интервала требует дополнительного предположения о распределении внутри интервала, а статья о конформном подходе прямо возражает против присоединения отношения интервала к $f^*$. См. критерий Келли для стратегий, чтобы понять, что действительно требуется для этой доли.

Заключение

Стек оценки вероятностного прогноза короток и не подлежит торгу: CRPS для оценки, потому что он является правильным правилом и выражен в единицах цели; гистограмма PIT для диагностики, потому что называет тип ошибки, а не просто сигнализирует о ней; покрытие по каждому горизонту для принятия решения, потому что именно это число отображается в размер позиции. Все, что называется quantile loss, — это CRPS с множителем 2, поэтому метрика здесь одна, а не две.

Неприятная часть в том, что весь этот инструментарий может показать: модель хуже, чем хотелось бы. В этом и смысл. Модель, которая честно расширяет интервалы, когда не знает ответа, строго полезнее модели, сохраняющей уверенно узкие интервалы, и отличить их можно только корректной оценкой.

Ссылки

• 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