Model selection with proper scoring rules on data sets of time series Giorgio Corani∗, Stefano Damato, Dario Azzimonti, Lorenzo Zambon SUPSI, Istituto Dalle Molle di Studi sull’Intelligenza Artificiale (IDSIA), Lugano, Switzerland
arXiv:2606.24715v1 [stat.ML] 23 Jun 2026
Abstract We consider the problem of model selection between probabilistic models on data sets of time series. Chosen a proper scoring rule, we denote by the term score the average value of the scoring rule on the test of an individual time series. For model selection, we need aggregating the values of the scores across multiple time series. Three summary statistics are commonly used for model selection: mean score, median score, and mean rank. Results in previous papers show that these statistics can yield conflicting decisions; we show how the conflicting conclusions are due to the skewness of the distribution of scores. We also show that as the test set of each time series of the data set increases, the different model selection criteria progressively converge to the same conclusion. However, for short tests sets, only the mean score identifies the true model as the best. We illustrate these phenomena with an analysis on intermittent time series, including the data set of the M5 competition, where we underline the importance of having a large test set. In such experiments, we further notice that model selection based on mean ranks remains unchanged using different scaling factors.
1. Introduction Scoring rules (Gneiting & Raftery, 2007; Gneiting & Katzfuss, 2014) quantify how good a probabilistic forecast is. A univariate scoring rule S is a function S(P, y) which assigns a numerical score to a probability distribution P , given the actual observation y ∈ R. We consider negatively oriented scoring rules, i.e., better forecasts corresponds to lower scores. A scoring rule is proper if it is minimized in expectation by the true distribution, i.e., if Ey∼Q [S(Q, y)] ≤ Ey∼Q [S(P, y)] for any valid forecast distributions P, Q. The most common proper scoring rules are the continuous ranked probability score (CRPS), the quantile score (QS), the interval score (IS), and the energy score (ES). In this paper, the test set refers to a set of nte out-of-sample data points, related to different rolling origins and/or forecast horizons; the pointwise score is the scoring rule evaluated on a single data point, S(P, yi ) Pni.e. te for i ∈ {1, . . . , nte }, and the score is the average value of the scoring rule on the test set, i.e. n1te i=1 S(P, yi ). Model selection on a single time series is straightforward: we select the model with the lowest score, assuming nte to be large enough for the score to approximate the expected value. We can instead adopt different approaches for model selection on a data set containing n time series: the most common approaches are the mean scaled score and the mean rank. In the first case, we scale the n scores to make them comparable; then we average them and we select the model with the lowest mean scaled score. While the need for scaling is established (Bolin & Wallin, 2023), the choice of the scaling factor is not obvious; most scaling factors (Hyndman & Koehler, 2006; Makridakis, Spiliotis & Assimakopoulos, 2022a; Athanasopoulos & Kourentzes, 2023; Svetunkov & Boylan, 2023) are related to the scale of the time series or to the scale of the error of a benchmark model. Generally however no scaling factor is adequate for all the non-stationarities and non-normalities which might be present in the different time series of a data set (Hewamalage, Ackermann & Bergmeir, 2023). ∗ Corresponding author
Email addresses: [email protected] (Giorgio Corani ), [email protected] (Stefano Damato), [email protected] (Dario Azzimonti), [email protected] (Lorenzo Zambon)
The perceived risk is thus that the outcome of model selection might depend on the chosen scaling factor; however, to our knowledge, there is yet no study on this point. The mean rank approach is appealing as it avoids the choice of the scaling factor. It ranks the models according to their score on the test set; then it averages the ranks over time series and selects the model with lowest mean rank. The significance of the differences between the mean ranks is tested with a Friedman test (Demšar, 2006) with Nemenyi post-hoc; this is also referred to as multiple comparison with the best (MCB) (Koning, Franses, Hibon & Stekler, 2005) in the forecasting literature. Yet, sometimes the mean scaled score and the mean rank lead to conflicting conclusions. For instance, Spiliotis, Makridakis, Kaltsounis & Assimakopoulos (2021, Fig. 4) report an experiment on the ∼30,000 bottom time series of the M5 competition (Makridakis et al., 2022a). On the high quantiles, the Poisson static distribution is competitive according to the mean rank, but its mean scaled score is higher (i.e., worst) than most competitors (Spiliotis et al., 2021, Tab. 1). It remains unclear whether it is preferable to perform model selection based on the mean scaled score or the mean rank, and which criterion to trust when they lead to conflicting conclusions; in this paper we fill this gap. The paper is organized as follows. In Sect. 2, we show that the distribution of the pointwise scores is typically strongly skewed. In Sect. 3, we show that also the distribution of the scores, and of their differences, is generally skewed for common values of nte . We also show that, given a skewed distribution of the scores, a misspecified model can have a better mean rank than the true model. The same problem is found for other non-mean statistics such as the median score and the win rate. In Sect. 4, we perform model selection on real data sets of intermittent time series, choosing between the Poisson and the negative binomial distribution. We discuss the cases in which the mean rank and the mean scaled score lead to different conclusions, why this happens and which criterion is arguably more reliable. We also study the sensitivity of the decisions of the mean scaled score to the scaling factor and eventually we provide our recommendations. 2. The skew of the distribution of the pointwise scores We now introduce the most common univariate scoring rules, defined to be negatively oriented. The quantile score with probability p, QSp (y, ŷp ), is: ( QSp (y, ŷp ) =
2(1 − p) ŷp − y 2(p) y − ŷp
if y < ŷp (overestimation), if y ≥ ŷp (underestimation),
(1)
where ŷp is the quantile at level p of the forecast distribution and y is the observation. Let us discuss some properties of the pointwise scores of Eq. (1), assuming p > 0.5; the case p < 0.5 is symmetric. They aee skewed since they are lower-bounded by zero and without upper bound. Moreover, the skew increases with p as we explain in the following. The pointwise scores are drawn from a mixture with two components. The overestimation component yields many small pointwise scores, while the underestimation component yields few large pointwise scores. Given an underestimation and an overestimation error of the same size, underestimation results in a pointwise scores which is k = 1−p p times larger than overestimation. For p = {0.8, 0.9, 0.99}, we have k = {4, 9, 99}. We also notice that, for well calibrated forecasts, the underestimation errors are k times rarer than the overestimation ones. The interval score (ISp ) (Gneiting & Raftery, 2007) for a prediction interval with coverage p is: 2 2 ISp (l, u, y) = (u − l) + (l − y)1{y < l} + (y − u)1{y > u} 1−p 1−p i 2 h = QS 1−p (y, l) + QS 1+p (y, u) 2 2 1−p
(2) (3)
where l and u are the lower and upper endpoints, and 1 is the indicator function, which is equal to 1 if the argument is true, and 0 otherwise. We assume a symmetric prediction interval, for which l and u are the 1+p quantiles at level 1−p 2 and 2 . Given the relation between ISq and QSp , in general also the skew of the distribution of the pointwise scores of ISq increases with p.
2
The Continuous Ranked Probability Score (CRPS) scores the entire predictive distribution and corresponds to the average of the quantile scores over all values of p: Z Z 2 CRPS(F, y) = (F (x) − 1{y < x}) dx = 2 QSp (y, ŷp )dp, (4) where F is the cumulative distribution function of the forecast distribution. More details on the relation between QS, IS and CRPS can be found in Fakoor, Kim, Mueller, Smola & Tibshirani (2023); Tibshirani (2023). p
0.5
0.8
0.9
0.95
0.99
ISp
1.9
3.7
5.5
8.2
19.5
QSp
1.0
2.2
4.1
6.9
18.3
CRPS
1.9
Table 1: Skewness of the distribution of the pointwise scores of different scoring rules on 105 data points. The N (0, 1) is used both as true distribution and as forecasting distribution.
To make the point concretely, we drawing 105 actual values y from the standard Gaussian N (0, 1) and we score the true distribution N (0, 1) on each of them. Thus, we use the N (0, 1) both as true distribution 3/2 and forecasting distribution. We then measure the skewness of the pointwise scores as M3 / M2 where Pn Mj = i (xi − x̄)j / n. The skewness is a dimensionless measure, which we classify as strong and extreme when its absolute value respectively exceeds 1 and 5 (Bulmer, 1979). In our experiment (Tab. 1), the pointwise scores of CRPS, ISp and QSp with p = 0.5 or p = 0.8 are strongly skewed; they are instead extremely skewed for p = 0.95 and p = 0.99. At the same level p, ISp is more skewed than QSp , since ISp = QS 1−p + 2 QS 1+p (Eq 3). 2
3. Model selection on a data set of n time series In the following simulation we perform model selection between two static distributions: the true model N (0, 1) and the misspecified model N (0, 0.852 ). We sample n = 105 test sets of length nte from the true model; on each test set we measure the average QS0.9 and QS0.99 of both models. We consider nte ∈ {4, 28, 100}. In this experiment no scaling is necessary, so the mean scaled score is simply the mean of the n scores. To implement the mean rank, on each time series we assign rank 1 and 2 to respectively the model with lower and higher score (there are no ties); we average the ranks over time series and we select the model with the lowest mean rank. CRPS
nte = 4
nte = 28
nte = 100
QS0.9
QS0.99
Mean score
Mean rank
Mean score
Mean rank
Mean score
Mean rank
True
0.56
1.45
0.35
1.59
0.05
1.91
Misspec.
0.57
1.55
0.36
1.41
0.06
1.09
True
0.56
1.31
0.35
1.38
0.05
1.58
Misspec.
0.57
1.69
0.36
1.62
0.06
1.42
True
0.56
1.16
0.35
1.29
0.05
1.37
Misspec.
0.57
1.84
0.36
1.71
0.06
1.63
Model selection over 105 simulations.
Table 2: true model is selected, and red otherwise.
For the selected models we highlight the background; we use green when the
As reported in Tab.2, on the CRPS both criteria select the true model. Yet, the mean rank selects the misspecified model when assessing QS0.9 or QS0.99 with small nte How is this possible? To start, notice 3
that for any p > 0.5 the misspecified model has a smaller ŷp than the true model. Thus frequently the misspecified model has lower pointwise score than the true one, since the overestimation error (which are the common ones) are penalized proportionally to ŷp − y. If the elicited quantile is high and/or nte is small, many test sets contain no underestimation errors; thus most scores are optimistically biased compared to the actual expected value of the scoring rule. By ranking independently the models on each test set the mean rank overlooks the importance the underestimation errors, which are present only on few test sets. The problem progressively disappears with larger nte , which makes the distribution of the score more symmetric around the expected value. For instance, dealing with QS0.99 , the mean rank selects the misspecified model for nte = 4 and nte = 28; eventually, it selects the true model for nte = 100. Because of the central limit theorem the distribution of the difference of the scores becomes symmetric as nte increases. However the higher the skewness of the population, the slower the convergence of the sample means to a normal distribution (Durrett, 2019, Chap. 3.4). For the common values of nte , the distribution of the score might still be skewed; we show that this can lead the mean rank and the mean scaled score to conflicting decisions. In our experiment, the CRPS does not appear to be problematic, since its pointwise scores and thus also its (differences of) scores are less skewed compared to QS0.9 and QS0.99 . Discussion B Let us denote by sA i and si the score of models A and B on time series i, and as di their difference, i.e., A B di = si − si . Let us that there P are no ties and no scaling factors. The mean scaled P assume1for P simplicity 1 B score selects A if n1 i sA i < n i si , which implies n i di < 0, i.e., that the mean difference is negative. A The rank of model A on time series i is r = 1 (d > 0) + 1, where 1Pis the indicator function. The mean i i P rank criterion selects A if r̄A = n1 i riA < 1.5, or equivalently, if n1 i 1(di > 0) < 0.5. This means that more than half of the times di < 0, i.e., the median difference is negative. The mean scaled score and the mean rank are thus estimators of different parameters, which can lead to the different decisions. We now re-analyze the simulation of Sect. 3 focusing on QS0.9 . In Fig. 1, the right panel shows the distribution of the scores of the true and misspecified models, while the left panel shows the distribution of their difference (true − misspecified), for different values of nte . For nte = 4 (top), the mean difference is negative, but the median is positive: the mean score correctly selects the true model, while the mean rank selects the misspecified model (Table 2). Indeed, the distribution of the differences di is strongly skewed and multimodal. The right peak corresponds to the simulations where no actual value y exceeds the predicted quantiles; in these cases, the score of the misspecified model is lower, because its overestimation errors are smaller. Since the test set is short, this happens frequently, causing the median difference to be positive. Since the test set is short, this happens on many test sets; thus, the median difference is positive. The other density peaks correspond to cases where, among the nte = 4 actual values, one or more fall in the right tail of the distribution. These cases correspond to negative score differences, as the underestimation errors of the true model are smaller, and move the mean of the distribution of the score differences below zero. As nte grows, the median is pushed towards the mean: the central limit theorem makes the distribution smoother and more symmetric, as the peaks become more numerous and closer to each other. Eventually, for nte = 100, the distribution is practically Gaussian. The mean rank is not a reliable criterion when nte is not sufficiently large. Proper scoring rules are minimized by the true distribution only in expectation: the mean rank correctly identifies the best model only if the score, obtained as the mean of the nte pointwise scores, is an accurate approximation of this expected value. For the quantile score, this requires having observed an adequate number of underestimation losses, which become more rare for higher quantiles. Determining exactly when nte is “large enough” remains difficult, as it depends on both the chosen scoring rule and the true distribution of the actual values. Therefore, we offer a general recommendation: practitioners should exercise caution when relying on mean ranks, especially when evaluating the quantile score at high quantile levels. In scenarios where model selection criteria yield conflicting results, the mean scaled score should generally be regarded as the more reliable metric. The same recommendations apply to other model selection criteria. For instance, the win rate, used by Shchur, Ansari, Turkmen, Stella, Erickson, Guerron, Bohlke-Schneider & Wang (2026) to compare foundation forecasting models, is the fraction of time series on which each model ranks first. When comparing two models, it is equivalent to the mean rank and thus suffers from the same limitations. With more than two models, the win rate is no longer equivalent to the mean rank, but it remains a function of the ranks alone
4
Quantile Score (alpha = 0.9): true vs misspec. distribution nte = 4
Misspec. True
−0.15
−0.10
−0.05
0.00
0.05
0.2
0.3
0.4
0.5
0.6
−0.10
−0.05
0.00
0.05
0.2
0.3
0.4
0.5
0.6
−0.10
−0.05
0.00
0.05
0.2
0.3
0.4
0.5
0.6
nte = 28
−0.15
nte = 100
−0.15
Score difference
Score
Figure 1: Distribution of the average QS over h = 4 (top) and h = 28 (bottom) steps ahead. The true distribution N (0, 1) is compared against a misspecified distribution N (0, 0.852 ). While the mean of the QS (solid line) is consistently lower for the true distribution, the median (dashed line) is higher when the test set is short. This corresponds to a mean rank exceeding 1.5, leading to wrong model selection.
and is therefore subject to the same issues. Since it estimates how often a model achieves the lowest score rather than estimating the expected score, it can disagree with the mean scaled score whenever nte is too small. The median scaled score can also be problematic. In the top right of Fig.1, we show that the true model has lower mean but higher median scaled score than the misspecified model. As nte grows, the median scores move closer the mean scores, eventually yielding the same model selection. We refer to Svetunkov & Boylan (2023, Tabs. 2–4, E.4) for an empirical example on the M5 bottom-level series (nte = 28) where mean and median scaled scores lead to conflicting conclusions at the 99th quantile. For example, the negative binomial model is among the best models according to the mean, but among the worst according to the median. Another family of criteria is based on relative scores, obtained by dividing the score of each model by that of a reference model and then aggregating across series, typically via a geometric mean. The skill score, defined as
5
the relative improvement over the reference model, is essentially equivalent. These criteria have the practical advantage of not requiring scaling factors, since the ratios are already dimensionless. However, since they combine the scores nonlinearly, they do not estimate the expected values of the scores; they are thus reliable only when nte is large enough for the scores to closely approximate their expectations. Wheatcroft (2019) reaches a similar conclusion, showing that the skill score is biased due to the ratio in its definition. Table A.5 in Appendix A illustrates these points on the same simulation as Sect. 3, reporting both the median and the relative score. We show that both criteria can yield wrong model selection, especially for small nte and for the quantile score with high quantile levels. 4. Empirical model selection between the Poisson and the negative binomial distribution We now compare the mean scaled score and the mean rank in a model selection experiment on different data sets of intermittent demand. We also check the sensitivity of the decision of the mean scaled score on the scaling factor. Our candidate models are the Poisson and the negative binomial distribution. The negative binomial is generally recognized (Syntetos, Babai & Altay, 2012; Kolassa, 2022) as a more suitable model than the Poisson for intermittent demand since, thanks to overdispersion, it can better estimate the probability of demand spikes (Boylan & Syntetos, 2021, Chap. 5). Experimental design. We consider the data sets listed in Tab. 3; as scoring rules we evaluate CRPS and several quantile scores (QS0.75 , QS0.835 , QS0.975 , QS0.995 ). We assess the impact of nte by varying it from 1 up to the full available test set length (28 for M5; 6 for Auto and 12 for RAF). Dataset
# of ts
freq.
T
max nte
avg. ADI
M5 (bottom only)
30490
D
1941
28
6.2
Auto
3000
M
18
6
1.3
RAF
5000
M
72
12
9.8
Table 3: The considered datasets; “D" stands for daily and “M" for monthly; T denotes the length of the training and test set. Avg. ADI is the average inter-demand interval, averaged over all time series of the data set. We took the M5 data set from the Kaggle page of the competition and the Auto and RAF data sets from Damato et al. (2026b).
Scaling factor
Formula PT −1 1 i=1 |yi+1 − yi | T −1 PT 1 i=1 yi T P T 1 i=1 S(yi , ED) T
MAE naive in-sample mean in-sample score of ED
Table 4: The considered scaling factors. We denote the training set as y1 , . . . , yT . The in-sample score of ED is the mean score achieved on the training set by the empirical distribution of the data.
We consider the scaling factors of Tab.4. The MAE naive was introduced as scaling factor of the MASE (Hyndman & Koehler, 2006). It was used to scale scoring rules in the M5 competition (Makridakis, Spiliotis, Assimakopoulos, Chen, Gaba, Tsetlin & Winkler, 2022b); its main shortcoming is that it might inflate the influence of the most intermittent series. The in-sample mean is an appropriate scaling factor for stationary time series but it might be inadequate for data sets containing a mixture of stationary and trending series (Athanasopoulos & Kourentzes, 2023; Hewamalage et al., 2023). It has been used for instance by Svetunkov & Boylan (2023). The in-sample score of ED scaling factor is constituted by the average value of the scoring rule on the training data, achieved by a reference model. As in Damato, Azzimonti & Corani (2026a), our reference model is the empirical distribution (ED) of the data, which is a sound probabilistic baseline for intermittent time series (Spiliotis et al., 2021; Kolassa, 2022; Long, Bui, Oktavian, Schmidt, Bergmeir, Godahewa, Lee, Zhao & Condylis, 2025). When the scoring rule is QSp , the scaling factor corresponds to the in-sample score of the empirical quantiles. 6
M5 data set QS0.975
CRPS
QS0.995
Mean rank
2.1
1.6
1.75
1.8
1.5
1.50
1.5
1.4
1.25
1.2
0.9 0.11
Mean score (scaled by MAE naive)
0.20
0.09
0.70 0.18 0.07 0.16
0.65
0.05
Score
0.14
Mean score (scaled by in−sample score of ED)
1.4 1.20
2.0
1.3 1.15 1.2 1.10
1.6
1.1 1.2
1.05
1.0
0.125
Mean score (scaled by in−sample mean)
0.90 0.250
0.100
0.85 0.225
0.075 0.80 0.200 0.050 5
10
15
20
25
5
10
15
20
25
5
10
15
20
25
Test set Model
NB
Pois
Figure 2: Comparison of the Poisson and the negative binomial distribution on on the M5 data set, considering different scoring rules, scaling factors and test set lengths.
We now analyze the results on the bottom time series of the M5 dataset. In Fig. 2 each column corresponds to a different scoring rule and each row to a different model selection criterion, namely the mean rank or the mean scaled score with various scaling factors. Depending on the row, we select the model with either the lower rank or the lower mean scaled score. Each panel shows how the mean ranks and mean scaled scores vary with nte . To improve readability, Fig. 2 only shows CRPS, QS0.975 and QS0.995 . In Fig D.5 we show the full results. On CRPS (Fig. 2, col.1) mean rank and mean scaled score with any scaling factor select the negative binomial at any value of nte . They also agree also on QS0.835 and QS0.75 (only shown in Fig D.5), this time selecting the Poisson. Indeed, different models can be preferable for different scoring rules. For QS0.975 , the decision of the mean rank criterion is sensitive on nte (Fig. 2, col.2, row 1): it selects the Poisson if provided with a short test set and instead the negative binomial if provided with a long enough (nte >22) test set. Indeed, the mean rank of the negative binomial and the Poisson has respectively a strong downward and upward trend with nte . Instead, the mean scaled score selects the negative binomial model for all nte with any scaling factor (Fig. 2, col.2, rows 2–4). We assume the negative binomial to be the correct decision, being both supported by larger test sets and consistent with the literature. The behavior of the mean rank criterion can be explained as the simulation of Sect. 3: given a negative binomial and a Poisson distribution with the same mean, the Poisson has a lower p-th quantile for p > 0.5 (although the two quantile forecasts may coincide in some cases due to discreteness). Dealing with high quantile levels and small nte , many test set contain no underestimation error; thus the Poisson might has
7
lower score than the negative binomial on them, obtaining eventually a better mean rank. As nte increases, however, the decisions of the mean rank converges to those of the mean scaled score. Thus the mean rank can require substantially larger nte than the mean scaled score to stabilize its decision when the distribution of the difference of the scores is strongly skewed, as it is typically the case of scoring rules with high levels of p. On QS0.995 the issue is even more apparent. Also in this case the mean ranks of the two models show a strong trends with nte : upward for the Poisson, downward for the negative binomial. We assume that, if we could further increase nte , the mean rank of the negative binomial would eventually become lower than that of the Poisson. However, the available test set is too short too observe this; as a result, the mean rank criterion selects the Poisson model for all nte (Fig. 2, row 1, col.3). Instead, the mean scaled score selects the negative binomial for all values of nte , with any scaling factor. This conflicting conclusions is reported also by (Spiliotis et al., 2021, Tab. 1 and Fig. 4), who considered nte = 28 and the mae naive scaling factor, and whose experiments included the negative binomial and the Poisson among many competing models. Indeed, they acknowledge the conflicting conclusions achieved by the mean rank and the mean scaled score. We extended their analysis by studying the sensitivity to nte and to the scaling factor. We report in Fig. B.3 and Fig. C.4 the results for the Auto and RAF data sets, which have smaller test sets (6 and 12 respectively). On both data sets we observe the same patterns. Mean rank and mean scaled score with any scaling factor take the same decisions on CRPS, QS0.75 and QS0.835 (only CRPS is shown in the figures). On QS0.975 and QS0.995 , on both data sets, the mean rank criterion selects the Poisson distribution at any nte ; yet, the negative binomial has a strong downward trend with nte , similarly to the previous cases. We assume that, with a larger test set, the mean rank criterion would eventually switch its decision in favor of the negative binomial, on both quantiles and on both data sets. Instead, on both quantiles and on both data sets, the mean scaled score selects the negative binomial at any nte and with any scaling factor. RAF is the dataset with highest ADI (Tab. 3) and thus the highest percentage of zeros. In this case, assuming the variance to be equal to the mean is especially limiting . Indeed, the mean scaled scored selects the negative binomial for any quantile, value of nte and scaling factor. Instead, the mean ranks selects the negative binomial distribution on QS0.975 and QS0.995 . In both cases, there is a strong downward trend of the mean rank of the negative binomial model, which suggest that decision of the mean rank could eventually switch to the negative binomial, if a longer test was available. In general, the decision of the mean scaled score is consistent at any nte and with any scaling factor. The estimated percentage difference between models is also broadly consistent across scaling factors, although it has some sensitivity on nte (in our experiments, the percentage improvements in mean scaled provided by the negative binomial over the Poisson tends to increase with nte ). 5. Conclusions Model selection with scoring rules requires estimating the expected value of the score. We have shown that the distribution of the scores (and of the difference of the score between competing models) is generally skewed, especially when dealing with small nte or scoring rules involving high quantiles (QSp or ISp with high p). If the distribution of the (difference of the) scores is skewed, non-mean statistics such as the mean rank, the win rate or the median score can pick the wrong model even if the data set contains a huge number of time series. The mean scaled score is a more appropriate summary statistic as it is actually an estimator of the expected value. In our experiments the decision of the mean scaled score are robust to the choice of scaling factor, whose choice appears to be less critical than previously assumed. We anyway recommend validating the conclusions by checking results across two or more scaling factors. An open problem is how to check the significance of the differences between competing models. Indeed, in large data sets the Diebold–Mariano test tends to reject the null hypothesis even with small differences in mean scaled scores. The problem might be addressed by adopting as Bayesian tests with region of practical equivalence (Kruschke, 2013). References Athanasopoulos, G., & Kourentzes, N. (2023). On the evaluation of hierarchical forecasts. International Journal of Forecasting, 39 , 1502–1511.
8
Bolin, D., & Wallin, J. (2023). Local scale invariance and robustness of proper scoring rules. Statistical Science, 38 , 140–159. Boylan, J. E., & Syntetos, A. A. (2021). Intermittent demand forecasting: context, methods and applications. Wiley. Bulmer, M. G. (1979). Principles of statistics. Courier Corporation. Damato, S., Azzimonti, D., & Corani, G. (2026a). Forecasting intermittent time series with Gaussian Processes and Tweedie likelihood. International Journal of Forecasting, in press. doi:https://doi.org/ 10.1016/j.ijforecast.2025.10.001. Damato, S., Zambon, L., & Azzimonti, D. (2026b). fable.intermittent: Forecasting Models for Intermittent Time Series. URL: https://CRAN.R-project.org/package=fable.intermittent. doi:10.32614/CRAN. package.fable.intermittent r package version 0.1.0. Demšar, J. (2006). Statistical comparisons of classifiers over multiple data sets. Journal of Machine Learning Research, 7 , 1–30. Durrett, R. (2019). Probability: Theory and Examples. (5th ed.). Cambridge University Press. Fakoor, R., Kim, T., Mueller, J., Smola, A. J., & Tibshirani, R. J. (2023). Flexible model aggregation for quantile regression. Journal of Machine Learning Research, 24 , 1–45. Gneiting, T., & Katzfuss, M. (2014). Probabilistic forecasting. Annual Review of Statistics and Its Application, 1 , 125–151. Gneiting, T., & Raftery, A. E. (2007). Strictly proper scoring rules, prediction, and estimation. Journal of the American statistical Association, 102 , 359–378. Hewamalage, H., Ackermann, K., & Bergmeir, C. (2023). Forecast evaluation for data scientists: common pitfalls and best practices. Data Mining and Knowledge Discovery, 37 , 788–832. Hyndman, R. J., & Koehler, A. B. (2006). Another look at measures of forecast accuracy. International Journal of Forecasting, 22 , 679–688. Kolassa, S. (2022). Commentary on the M5 forecasting competition. International Journal of Forecasting, 38 , 1562–1568. Koning, A. J., Franses, P. H., Hibon, M., & Stekler, H. (2005). The M3 competition: Statistical tests of the results. International Journal of Forecasting, 21 , 397–409. Kruschke, J. K. (2013). Bayesian estimation supersedes the t test. Journal of experimental psychology: General , 142 , 573. Long, X., Bui, Q., Oktavian, G., Schmidt, D. F., Bergmeir, C., Godahewa, R., Lee, S. P., Zhao, K., & Condylis, P. (2025). Scalable probabilistic forecasting in retail with gradient boosted trees: A practitioner’s approach. International Journal of Production Economics, 279 , 109449. Makridakis, S., Spiliotis, E., & Assimakopoulos, V. (2022a). M5 accuracy competition: Results, findings, and conclusions. International Journal of Forecasting, 38 , 1346–1364. Makridakis, S., Spiliotis, E., Assimakopoulos, V., Chen, Z., Gaba, A., Tsetlin, I., & Winkler, R. L. (2022b). The M5 uncertainty competition: results, findings and conclusions. International Journal of Forecasting, 38 , 1365–1385. Shchur, O., Ansari, A. F., Turkmen, C., Stella, L., Erickson, N., Guerron, P., Bohlke-Schneider, M., & Wang, Y. (2026). fev-bench: A realistic benchmark for time series forecasting. URL: https://arxiv.org/abs/ 2509.26468. arXiv:2509.26468.
9
Spiliotis, E., Makridakis, S., Kaltsounis, A., & Assimakopoulos, V. (2021). Product sales probabilistic forecasting: An empirical evaluation using the M5 competition data. International Journal of Production Economics, 240 , 108237. Svetunkov, I., & Boylan, J. E. (2023). iETS: State space model for intermittent demand forecasting. International Journal of Production Economics, 265 , 109013. Syntetos, A. A., Babai, M. Z., & Altay, N. (2012). On the demand distributions of spare parts. International Journal of Production Research, 50 , 2101–2117. Tibshirani, R. J. (2023). Forecast scoring and calibration. Lecture notes. URL: https://www.stat. berkeley.edu/~ryantibs/statlearn-s23/lectures/calibration.pdf. Wheatcroft, E. (2019). Interpreting the skill score form of forecast performance metrics. International Journal of Forecasting, 35 , 573–579. Appendix A. Median and relative scores for the simulation of Sec.3
CRPS
nte = 4
nte = 28
nte = 100
QS0.9
QS0.99
Median score
Rel. score
Median score
Rel. score
Median score
Rel. score
True
0.530
1.000
0.316
1.000
0.047
1.000
Misspec.
0.533
0.995
0.305
0.988
0.042
0.934
True
0.560
1.000
0.342
1.000
0.048
1.000
Misspec.
0.563
1.005
0.348
1.014
0.045
1.024
True
0.563
1.000
0.349
1.000
0.050
1.000
Misspec.
0.567
1.006
0.356
1.018
0.054
1.061
Table A.5: Model selection over 105 simulations. For the selected models we highlight the background; we use green when the true model is selected, and red otherwise. For short test sets and/or high quantiles, both the median score and the relative score select the misspecified model. Similar to the mean rank, on QS0.99 even for nte = 28 the wrong model is selected. On the contrary, the true distribution has a better CRPS for any test set size. Indeed, eq. (4) shows that it takes into account all quantile levels, balancing the effect of the QSp where underforecasts are rare yet strongly penalised.
10
Appendix B. Results on Auto data sets Auto data set QS0.75
CRPS
QS0.835
QS0.975
QS0.995 1.9
1.60
1.7
1.55
1.6
1.7
1.75
1.50
1.5
1.5
1.50
1.45
1.4
1.3
1.25
1.40
1.3
Mean rank
1.7 1.6 1.5 1.4 1.3 1.1
0.65 0.34
Mean score (scaled by MAE naive)
0.32
0.60
0.30
0.10
0.04
0.28
0.08 2.6
0.02
2.4
6
2.2
5
0.35
1.3 1.25 1.3 1.2
1.20
2.0 1.2
1.15
4
1.8 3
1.1 1.6
1.10 1.1
2
1.4
Mean score (scaled by in−sample mean)
0.14 0.56
0.35
0.30
0.07 0.12
0.52
0.33
0.28
0.05 0.10
0.48
0.31
0.26
0.03
0.08
1
2
3
4
5
6
1
2
3
4
5
6
1
2
3
4
5
6
1
2
3
4
5
6
Test set Model
NB
Pois
Figure B.3: Mean ranks and mean scaled scores on the Auto data set.
11
1
2
3
4
5
6
Score
1.30
Mean score (scaled by in−sample score of ED)
0.06
0.12
0.37
0.55
0.08
0.14
0.39
Appendix C. Results on RAF data set RAF data set QS0.75
CRPS
QS0.835
QS0.975
QS0.995
Mean rank
1.9
1.8
1.7
1.5
1.5
1.2
1.3
1.75
1.50
1.25
1.1
2.00
2.1
1.75
1.8
1.50
1.5
1.25
1.2
1.00
0.9
0.60 0.55
Mean score (scaled by MAE naive)
0.7
0.40
0.25
0.6
0.20
0.40
0.45
0.35
0.40
0.32
0.15 0.10
1.7
1.4
5
1.3
1.4
1.6 1.3
1.2
4
1.5
1.2 1.2 1.1
1.4
1.1 1.0
1.0
1.0
1.0
1.2
2
1.2
1.1
1.3
Mean score (scaled by in−sample mean)
1.3
0.9
0.9
3
1.0
0.75
0.6
0.70
0.5
0.65
0.4
0.60
0.3
0.9 1.1 0.9
0.8
1.0 0.9
0.8
0.7
0.2
0.55 2
4
6
8
10
0.8 12
2
4
6
8
10
12
2
4
6
8
10
12
2
4
6
8
10
12
Test set Model
NB
Pois
Figure C.4: Mean ranks and mean scaled scores on the RAF data set.
12
2
4
6
8
10
12
Score
Mean score (scaled by in−sample score of ED)
0.36
0.50
0.45
0.5
0.30
0.55
0.50
Appendix D. Full results on M5 data set M5 data set QS0.75
CRPS
QS0.835
QS0.975
QS0.995
1.58
2.1
1.52
Mean rank
1.6
1.51
1.5
1.50
1.49
1.4
1.54
1.75
1.8
1.50
1.50
1.5
1.46
1.25
1.2
1.48 0.9 0.11
1.42
0.57
0.50 0.20
0.09
Mean score (scaled by MAE naive)
0.54 0.48
0.70
0.18 0.51
0.07
0.46 0.16
0.65
0.05
0.44
0.48
Score
0.14
Mean score (scaled by in−sample score of ED)
0.42 1.4
1.20 1.20
1.15
2.0
1.3
1.15 1.15
1.10
1.2
1.6
1.10 1.10
1.1
1.05
1.2
1.05 1.05
1.0 1.00 1.00 0.125
Mean score (scaled by in−sample mean)
0.90 0.69
0.250
0.63
0.100
0.85
0.66 0.60
0.225
0.57
0.200
0.075
0.63
0.80
0.60 5
10
15
20
25
5
10
15
20
25
0.050 5
10
15
20
25
5
10
15
20
25
Test set Model
NB
Pois
Figure D.5: Mean ranks and mean scaled scores on the M5 data set.
13
5
10
15
20
25