Scoring Probabilistic Forecasts: CRPS, PIT Calibration, and 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 등, 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)은 추론 시 드롭아웃을 활성화한 채 $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$

CRPS가 중요한 자리를 차지하는 이유는 세 가지다.

• 적절한 점수 규칙이다(Gneiting과 Raftery, 2007). 예측 분포가 실제 분포와 같을 때에만 기대값이 최소가 된다. 의도적인 과신이나 인위적으로 넓은 분포를 이용한 회피는 점수를 높이지 못한다. 중위수에 대한 MSE에는 이 성질이 없으므로 CRPS는 선택 사항이 아니다.
• MAE를 일반화한다. 퇴화한 점 예측에서는 절대 오차로 줄어들므로 CRPS는 타깃의 단위를 사용한다. 시간별 로그 수익률을 예측할 때 CRPS 0.004는 40bp 평균 절대 오차와 직접 비교할 수 있으며, 감으로 점검할 수 없는 무차원 수가 아니다.
• 보정된 상태에서 예리함을 보상한다. 똑같이 보정된 두 예측 중 더 좁은 예측의 점수가 높다. 따라서 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자형: 과신——꼬리에 너무 많은 질량이 놓여 분포가 지나치게 좁다.
• 봉우리형: 과소신——관측값이 중심에 몰려 분포가 지나치게 넓다.
• 왼쪽으로 치우침: 모델이 체계적으로 과대예측한다.
• 오른쪽으로 치우침: 모델이 체계적으로 과소예측한다.

공동 위험의 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 절을 보라.

계산 비용. 500개 샘플 경로를 생성하는 DeepAR는 단일 순전파의 약 $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). "베이즈 근사로서의 드롭아웃: 딥러닝에서 모델 불확실성 표현하기." 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