Multivariate Probabilistic Forecast Evaluation
Scoring, calibration, tail diagnostics, and dependence-aware inference for joint predictive distributions
QM015 · Forecast Evaluation · Advanced
Core idea. A multivariate probabilistic forecast should be judged as a distribution, not as a point estimate: evaluate location, dispersion, calibration, cross-variable dependence, and decision-relevant tails with complementary diagnostics.
Use it for. Comparing joint scenario generators, density forecasts, ensemble forecasts, and Monte Carlo predictive distributions for several assets or variables at once.
It does not establish. That the model is “good” because one score improves, that nominal interval coverage proves correct joint dependence, or that model-implied tail losses are guaranteed realized-loss forecasts.
The Question
A point forecast answers a question such as “What is the expected return?” A probabilistic forecast answers a larger question: what distribution of outcomes is plausible, and how are those outcomes expected to move together?
For a multivariate forecast, the object of interest may be a joint distribution
\[ F_t(\mathbf y)=\Pr(\mathbf Y_t\le \mathbf y\mid \mathcal I_t), \]
where \(\mathbf Y_t\in\mathbb R^d\) contains several future quantities and \(\mathcal I_t\) is the information set available when the forecast is made.
The evaluation problem is therefore not reducible to RMSE on a conditional mean. A model can have reasonable marginal forecasts but poor dependence, acceptable average scores but weak tail calibration, or good interval coverage with unnecessarily diffuse predictions. A robust evaluation design separates these dimensions.
Why It Matters
Probabilistic forecasting is intended to represent uncertainty. Proper scoring rules provide a decision-theoretic way to compare predictive distributions, while calibration diagnostics ask whether realized outcomes are statistically consistent with those distributions (Gneiting and Raftery 2007; Gneiting et al. 2007).
In multivariate applications, dependence matters because a portfolio, system, or joint decision depends on combinations of variables rather than on each marginal distribution independently. Scheuerer and Hamill show that the energy score can be comparatively insensitive to some dependence misspecification and motivate the variogram score as a complementary multivariate diagnostic (Scheuerer and Hamill 2015).
This leads to a general principle:
Use a scorecard, not a single winner metric.
A single number can summarize one aspect of predictive quality, but it cannot reveal every failure mode that matters in a joint scenario model.
Intuition
Suppose two models generate one thousand 20-day scenarios for eight assets.
- Model A gets each asset’s marginal volatility roughly right but makes the assets too independent.
- Model B reproduces cross-asset co-movement better but is slightly too diffuse in some marginals.
A univariate score computed asset by asset may prefer A. A dependence-sensitive diagnostic may prefer B. A portfolio-tail calibration check may reveal which misspecification matters more for the intended risk decision.
The correct conclusion is not that one diagnostic is wrong. The diagnostics answer different questions.
The Method
Point forecast versus probabilistic forecast
For a point forecast \(\hat{\mathbf y}_t\), evaluation typically applies a loss such as squared error:
\[ L_t^{\mathrm{MSE}}=\|\mathbf y_t-\hat{\mathbf y}_t\|_2^2. \]
For a probabilistic forecast \(F_t\), the input to the evaluator is the entire predictive distribution. In simulation-based work, \(F_t\) is often represented by an ensemble
\[ \mathbf X_t^{(1)},\ldots,\mathbf X_t^{(M)}\sim F_t. \]
The realized outcome \(\mathbf y_t\) is then compared with the predictive cloud.
Proper scoring rules
A scoring rule is proper when, in expectation, a forecaster does not improve the score by reporting a distribution different from the distribution it actually believes. Strict propriety strengthens this condition so that the truthful distribution is uniquely optimal under the score (Gneiting and Raftery 2007).
This article uses the common loss orientation:
lower score = better forecast.
Some literature writes scores as rewards to maximize. The orientation must be checked before comparing values across implementations.
Energy Score
For a \(d\)-dimensional predictive distribution \(F\) and realization \(\mathbf y\), the Energy Score is
\[ \operatorname{ES}(F,\mathbf y) = \mathbb E\|\mathbf X-\mathbf y\|_2 - \frac12\mathbb E\|\mathbf X-\mathbf X'\|_2, \]
where \(\mathbf X\) and \(\mathbf X'\) are independent draws from \(F\) (Gneiting and Raftery 2007).
For an ensemble of \(M\) draws, a direct empirical approximation is
\[ \widehat{\operatorname{ES}} = \frac1M\sum_{m=1}^M\|\mathbf X^{(m)}-\mathbf y\|_2 - \frac{1}{2M^2}\sum_{m=1}^M\sum_{\ell=1}^M \|\mathbf X^{(m)}-\mathbf X^{(\ell)}\|_2. \]
A finite-ensemble implementation may instead use a leave-diagonal-out or “fair” estimator. The chosen estimator should be stated rather than silently mixed across models.
The Energy Score rewards a predictive cloud that is close to the realization without being needlessly dispersed. It is a useful multivariate generalization of the CRPS, but it should not be treated as a complete dependence diagnostic (Scheuerer and Hamill 2015).
Variogram Score
For order \(p>0\), one form of the Variogram Score is
\[ \operatorname{VS}_p(F,\mathbf y) = \sum_{i=1}^d\sum_{j=1}^d w_{ij} \left( |y_i-y_j|^p - \mathbb E_F|X_i-X_j|^p \right)^2, \]
with nonnegative pair weights \(w_{ij}\) (Scheuerer and Hamill 2015).
The score targets pairwise structure directly. In an ensemble, the expectation is estimated by averaging \(|X_i^{(m)}-X_j^{(m)}|^p\) across draws. The choice of \(p\) and \(w_{ij}\) is part of the specification.
The Variogram Score is especially useful when dependence or relative movement across components is important, but it is not a universal substitute for marginal calibration or tail diagnostics.
CRPS and standardized CRPS
For a univariate predictive CDF \(F_j\) and realization \(y_j\), the continuous ranked probability score can be written as
\[ \operatorname{CRPS}(F_j,y_j) = \int_{-\infty}^{\infty} \left(F_j(z)-\mathbf 1\{y_j\le z\}\right)^2dz. \]
It is also a proper score and has the same units as the variable (Gneiting and Raftery 2007; Hersbach 2000).
In a multi-asset problem, raw CRPS values can be dominated by components with larger natural scales. A standardized CRPS can therefore be defined as
\[ \operatorname{sCRPS}_{j,t} = \frac{\operatorname{CRPS}_{j,t}}{s_j}, \]
where \(s_j>0\) is a pre-specified scale such as a training-sample standard deviation or another fixed robust scale.
The normalization scale must be fixed from admissible information and reported explicitly. Re-estimating the denominator with future data or choosing a denominator after inspecting the results changes the evaluation target.
Marginal CRPS is useful even in a multivariate model because it reveals which individual variables are poorly calibrated or excessively diffuse. It does not verify the cross-variable dependence structure.
PIT and calibration
For a continuous univariate predictive CDF \(F_t\), the probability integral transform is
\[ u_t=F_t(y_t). \]
Under a correctly specified conditional density forecast, PIT values have a Uniform\((0,1)\) benchmark, with stronger sequential conditions relevant to independence as well (Diebold et al. 1998; Gneiting et al. 2007).
A PIT histogram can therefore detect systematic U-shapes, humps, or asymmetry consistent with under-dispersion, over-dispersion, or location bias. But a nearly flat histogram does not prove correct multivariate dependence.
For ensemble forecasts, randomized or rank-based PIT variants may be required because the empirical predictive distribution is discrete.
Interval coverage
For a nominal \((1-\alpha)\) prediction interval \([L_t,U_t]\), empirical coverage is
\[ \widehat C_{1-\alpha} = \frac1T\sum_{t=1}^T \mathbf 1\{L_t\le y_t\le U_t\}. \]
Calibration suggests that empirical coverage should be close to nominal coverage over an appropriate evaluation sample. Coverage should be inspected at more than one nominal level, because a model may look reasonable at 80% but fail in the 95% tail.
Coverage alone is incomplete. Extremely wide intervals can attain high coverage while being uninformative. Proper scores help balance calibration and sharpness (Gneiting et al. 2007).
Portfolio-tail calibration
When scenarios are ultimately used for portfolio risk, the joint predictive cloud can be projected through a pre-specified portfolio map. Let \(R_{p,t}^{(m)}\) be the portfolio return implied by scenario \(m\) and let \(R_{p,t}^{\mathrm{real}}\) be the realized return.
Examples of portfolio-tail calibration diagnostics include:
- empirical exceedance frequency of a model-implied VaR threshold;
- realized severity conditional on tail exceedance;
- coverage of model-implied portfolio intervals; and
- calibration of path-based events such as maximum drawdown exceeding a threshold.
These checks test whether the joint model is useful after being mapped into a decision-relevant quantity. They are not substitutes for asset-level marginal and dependence diagnostics.
Dependence diagnostics
A multivariate forecast should also be challenged on the structure connecting its components. Useful diagnostics include:
- Variogram Score;
- forecast-versus-realized covariance or correlation patterns;
- rank correlation;
- co-exceedance frequencies in adverse tails;
- conditional cross-asset behavior during large moves; and
- portfolio projections that are sensitive to diversification breakdown.
Some of these are descriptive diagnostics rather than proper scores. That distinction should be stated.
Why One Metric Is Not Enough
Different metrics have different blind spots.
| Diagnostic | Primary target | Important limitation |
|---|---|---|
| Energy Score | Overall multivariate distribution | Can have limited sensitivity to some dependence errors |
| Variogram Score | Pairwise dependence geometry | Does not replace marginal calibration checks |
| CRPS / sCRPS | Marginal probabilistic accuracy | Ignores joint dependence when used componentwise |
| PIT | Marginal distributional calibration | Flatness alone does not establish joint correctness |
| Interval coverage | Quantile / interval calibration | Can reward overly wide intervals if used alone |
| Portfolio-tail calibration | Decision-relevant joint tail behavior | Depends on the chosen portfolio and risk map |
A credible evaluation therefore reports a joint scorecard and interprets disagreements rather than suppressing them into one rank.
Overlapping Forecast Horizons and Serial Dependence
Suppose a 20-day forecast is issued every trading day. The target for origin \(t\) uses returns from roughly \(t+1\) through \(t+20\), while the target for origin \(t+1\) uses \(t+2\) through \(t+21\).
Those realized target windows share 19 days. Consequently, score sequences, PIT-based statistics, coverage indicators, and model-loss differentials can be serially dependent even when the forecast procedure itself is unchanged.
This overlap matters for inference. A naive standard error that treats \(T\) forecast origins as independent can materially overstate the effective sample size.
Dependence-aware options include:
- heteroskedasticity-and-autocorrelation-consistent covariance estimators (Hansen and Hodrick 1980; Newey and West 1987);
- moving-block or other dependence-preserving bootstrap procedures; and
- sensitivity analysis over plausible lag or block-length choices.
QM006 — Block Bootstrap develops one practical route.
For an \(h\)-period target evaluated every period, overlap can generate dependence through at least the shared target window. Re-estimation, persistent state variables, and market dynamics can create additional dependence beyond that mechanical overlap.
How to Interpret the Result
A lower average Energy Score or CRPS means the model performed better under that scoring rule’s loss orientation over the evaluated origins. It does not imply that every calibration dimension improved.
A better Variogram Score supports the narrower conclusion that the evaluated pairwise structure was closer under the chosen \(p\) and weights.
A well-calibrated PIT or interval-coverage plot supports consistency with those marginal calibration checks. It does not establish accurate portfolio tails unless the portfolio transformation is also evaluated.
When model A beats model B by a small average score difference, the uncertainty of that difference matters—especially with overlapping horizons. Statistical uncertainty should be assessed with a dependence-aware procedure rather than by counting origins as independent draws.
Financial / Economic Example
Assume two scenario models forecast daily returns for a stock ETF and a bond ETF.
The realized one-day return is
\[ \mathbf y=(-2.0\%,\,+0.5\%). \]
Model A generates scenarios with realistic marginal volatility but almost zero stock-bond dependence. Model B produces a moderate negative dependence that is closer to recent data, but its bond marginal is slightly too wide.
A componentwise CRPS comparison can prefer Model A because its bond marginal is sharper. The Variogram Score can prefer Model B because the stock-bond spread geometry is more realistic. If the portfolio is 60/40, portfolio-tail calibration can then show which model better reproduces the relevant loss distribution.
This example explains why score disagreement is informative. It is not empirical evidence that one particular model family is superior.
Implementation
The included Python example computes an empirical Energy Score, Variogram Score, ensemble CRPS, PIT rank, and interval coverage from synthetic draws. A minimal Energy Score implementation is:
import numpy as np
def energy_score(samples, y):
x = np.asarray(samples, dtype=float)
y = np.asarray(y, dtype=float)
term1 = np.linalg.norm(x - y, axis=1).mean()
pairwise = np.linalg.norm(x[:, None, :] - x[None, :, :], axis=2)
term2 = 0.5 * pairwise.mean()
return term1 - term2For large \(M\), the full \(M\times M\) pairwise matrix can be expensive. Chunking, subsampling, or algebraically equivalent estimators may be needed, but the estimator used should remain comparable across models.
Common Mistakes
Treating a probabilistic forecast as a point forecast. Comparing only predictive means discards dispersion and tail information.
Calling one score a complete model-quality measure. Proper scores are principled, but each score still emphasizes particular aspects of the forecast.
Using marginal CRPS to certify joint dependence. Componentwise accuracy does not prove that assets co-move correctly.
Treating PIT uniformity as the whole calibration problem. PIT is powerful but principally marginal unless extended into a deliberate multivariate diagnostic framework.
Ignoring score orientation. Some sources maximize rewards; many implementations minimize losses. A sign error can reverse a model ranking.
Standardizing with future information. Scaling CRPS using the full sample leaks future scale information into historical evaluation.
Using iid standard errors for overlapping targets. Shared future windows mechanically induce dependence in evaluation sequences.
Interpreting model-implied tail metrics as guaranteed realized losses. They are conditional statements under the forecast distribution and its information set.
When Not to Use It
If the research question is purely about a conditional mean point forecast, simpler point-forecast evaluation may be sufficient. If the predictive object is univariate, multivariate scores add no value beyond a well-designed univariate probabilistic evaluation.
Do not use a multivariate score merely because the dataset has several columns. The forecast itself must represent a joint predictive object for the score to have its intended meaning.
Used in SlackQuant Research
- Multi-Asset Scenario Stress Lab — Technical White Paper (SSRN 7354238): uses joint probabilistic forecast evaluation as part of a broader assessment of multi-asset scenario quality and stress-decision usefulness.
Reproducibility
The included synthetic Python lab reproduces the score definitions and demonstrates a deliberately constructed case in which marginal and dependence diagnostics disagree. No empirical performance claim is attached to the synthetic outputs.