cultivars.forecast.scoring#
Proper scores for density forecasts: what a predictive sample earned.
A density forecast is a claim about a distribution, and only a strictly proper scoring rule – one a forecaster cannot improve by hedging away from their honest belief – deserves to grade it (Gneiting & Raftery 2007). That is the admission criterion here. For a predictive distribution \(F\) and an outcome \(y\), a score \(S(F, y)\) is strictly proper when \(\mathbb{E}_{Y \sim G}\, S(F, Y)\) is minimized over \(F\) uniquely at \(F = G\); each score below is computed from the predictive sample \(x^{(1)}, \dots, x^{(S)}\) that stands in for \(F\). The CRPS grades each series exactly from the sample,
the energy score is its multivariate generalization with the Euclidean norm in place of the absolute value, grading the joint claim; the log score \(-\log f(y)\) is the classical alternative, and the one number in this module that is estimated (by a kernel density on the draws) rather than exact, which its docstring says plainly. The point losses ride along as degenerate cases because backtests report them, not because they suffice: a point forecast is a density forecast that has thrown its uncertainty away. The pinball loss is the proper score for the one point forecast that keeps a distributional claim, a quantile.
Two commitments shape the surface. First, everything scores raw arrays
– forecast_paths output, or any rival forecaster’s simulations –
so the comparison is never restricted to this package’s own models.
This module carries a load the design placed deliberately: the Gibbs
BVAR family reports no marginal likelihood, so out-of-sample density
scores are the honest instrument that ranks a conjugate fit against a
horseshoe fit against anything else that produces predictive paths.
Second, all scores are negatively oriented (smaller is better), and
everything is reported per horizon and per series before any averaging,
because aggregation is a modelling choice the user should make in
daylight; the summary averages for display and says so.
Layout. DensityScore scores one origin and returns a
DensityScoreResult with the per-cell arrays; pinball_loss()
scores a quantile forecast elementwise. The numerics live in _core:
_crps_from_draws is the exact sample CRPS, _energy_score the exact
sample energy score by chunked pairwise norms, _kernel_log_score the
Gaussian-kernel log density with Silverman’s bandwidth, and
_pinball_loss the check function. Across many origins the same
scores are collected by
BacktestResult.
References
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.
Gneiting, T. (2011). Making and evaluating point forecasts. Journal of the American Statistical Association, 106(494), 746-762.
Example
A conjugate and a horseshoe BVAR fitted to the same 120 observations and scored on the four that follow, the comparison the module exists for:
>>> import numpy as np
>>> from cultivars.bayes.priors import HorseshoePrior
>>> from cultivars.multivariate.large_dim.bayesian import BVAR
>>> from cultivars.multivariate.large_dim.gibbs import GibbsBVAR
>>> y = np.random.default_rng(0).standard_normal((124, 2))
>>> conjugate = BVAR(y[:120], order=1).fit(n_draws=200, seed=0)
>>> sampler = GibbsBVAR(y[:120], order=1, prior=HorseshoePrior())
>>> horseshoe = sampler.fit(n_draws=400, n_burn=200, seed=0)
>>> first = DensityScore(conjugate.forecast_paths(4, seed=0), y[120:]).compute()
>>> second = DensityScore(horseshoe.forecast_paths(4, seed=0), y[120:]).compute()
>>> first.crps.shape, second.crps.shape
((4, 2), (4, 2))
>>> bool(np.all(first.crps > 0.0)) and bool(np.all(second.energy > 0.0))
True
Classes#
Score one origin's density forecast against what happened. |
|
One forecast origin's scores, per horizon and per series. |