Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
90,054 characters · 17 sections · 74 citation commands
Uncertainty Quantification in Forecast Comparisons
\noindentKeywords: Confidence Bands, Forecast Evaluation, Predictive Performance, Probabilistic Forecasting
Forecasts are the basis for sound decision-making. Consequently, forecasting is an important endeavor in many disciplines such as medicine, meteorology, economics, or finance. Traditionally, forecasts for the central tendency (mean or median) of a variable of interest play a key role. In recent years, other types like quantile, interval or probabilistic forecasts, that is, forecasts for the full probability distribution, have become increasingly popular as they give a more detailed picture, and, in particular, inform about uncertainty gneiting2014.
To ensure that forecasts are of high quality and to improve future forecasts, suitable evaluation methods are needed. The standard approach to (relative) forecast evaluation is to compare the performance of competing forecasts via suitable loss functions tailored to the type of forecast at hand: consistent scoring functions gneiting2011 or proper scoring rules gneiting2007proper. Average scores over an evaluation sample of forecasts and observations provide a ranking of different forecasting methods in terms of their accuracy and allow for a principled choice of forecasting approaches.
As the magnitude of average scores is usually not interpretable, it is common to compute relative average scores, that is, to divide the average score of a forecasting method of interest by the average score of a benchmark method. In this way, the relative improvement in forecast accuracy over the benchmark can be easily read off. Often, in particular in the meteorological literature, this relative improvement is computed directly, and called skill score gneiting2007proper, wilks2011, thorarinsdottir2018. Hence, relative average scores, or skill scores, not only allow for a ranking in terms of forecast accuracy, but also for an assessment of the magnitude of the difference in forecasting accuracy and thus of its relevance. For example, a skill score of 10% of a new method compared to a simple benchmark, which means an increase in forecast accuracy by 10%, might lead to a decision to implement that method. In contrast, with a skill score of only 1% the additional effort might not be worthwhile, depending on the application.
Average scores, relative average scores and empirical skill scores are of course estimates of the respective population quantities (expected scores, relative expected scores and population skill scores) based on an evaluation sample of forecasts and observations. Thus, they should be accompanied by appropriate measures of sampling uncertainty. The classical and most natural way to quantify and communicate sampling uncertainty are confidence intervals. However, such confidence intervals are usually not reported in the forecast evaluation literature. Instead, p-values from Diebold-Mariano tests of forecast accuracy diebold1995 are very popular. One reason why test results and not the corresponding confidence intervals are reported is certainly that this test is concerned with differences in expected scores and due to their lack of interpretability, the corresponding confidence intervals would also lack interpretability. However, confidence intervals for relative expected scores or skill scores would yield the same conclusions about the null hypothesis of equal predictive accuracy and would additionally provide interpretable statements about sampling uncertainty by reporting the percentage improvements in forecasting performance supported by a certain confidence level.
In practice, usually not only a single comparison between two forecasting methods in terms of their forecast accuracy is executed, but many: often multiple forecast horizons, variables, forecasting methods and/or locations are of interest. This leads to large tables with many p-values of pairwise forecast accuracy tests, which are, of course, plagued by serious multiple testing problems, often rendering the conclusions from those tests essentially useless. Analogously, pointwise confidence intervals, which are intended for a single comparison, are invalid and usually too narrow in such settings. Existing approaches addressing specific instances of this problem from a hypothesis testing perspective are model confidence sets hansen2011 for the case of multiple forecasting methods, and a multi-horizon Diebold-Mariano test quaedvlieg2021.
A natural solution is to replace pointwise intervals with simultaneous confidence bands for a vector of population skill scores (or expected scores), which contain this vector with a prespecified probability of $1-\alpha$. Such bands carry the advantages of confidence intervals to the multivariate setting. They give a valid statement of sampling uncertainty, report the skill score vectors supported by the data at a given confidence level, and, contrary to other types of confidence regions such as ellipsoids, they can be easily visualized regardless of dimension. As a byproduct, they further automatically provide test results for relevant null hypotheses. For example, to test the hypothesis of equal predictive accuracy of two methods over a set of forecast horizons at significance level $\alpha$, one just needs to check if zero lies inside the corresponding confidence band of skill scores of level $1-\alpha$ for all those horizons.
The construction of confidence intervals or confidence bands for skill scores has been considered in the meteorological literature. However, the discussion of this topic is rather limited and often subject to deficiencies. First, the focus is almost exclusively on pointwise confidence intervals stephenson2000, ahrens2008, bradley2008. Second, observations are usually assumed to be independent, making the approaches invalid for time series data. Noteworthy exceptions regarding the second point are baran2023 and baran2024, who report pointwise confidence bands constructed via a block bootstrap, thus accounting for temporal dependence.
To address these limitations, we introduce two types of simultaneous confidence bands, sup-t and Bonferroni, for vectors of skill scores (as well as relative expected scores and expected scores themselves). We combine recent results on confidence bands in a nonlinear setting and novel bootstrap implementations by montielolea2019 with assumptions that can be seen as multivariate versions of the classical Diebold-Mariano assumptions diebold1995,diebold2015 to establish the validity of our bands: Both types of bands have asymptotic coverage of at least the nominal coverage. While the sup-t bands even have exact asymptotic coverage, the Bonferroni bands are mildly conservative. Our proposed implementation involves resampling of scores via a moving block bootstrap kunsch1989 for time series data and via a classical bootstrap for iid data. Our approaches are easy to implement and to apply, but at the same time very versatile: They allow to construct confidence bands over multiple forecast horizons, variables, forecasting methods, locations or combinations thereof. The bands can be easily represented graphically even if they span several of those dimensions. Thus, they are a beneficial addition to the usual graphs showing average scores or skill scores over, for example, forecast horizons for multiple methods. Furthermore, they allow for testing many hypotheses of interest directly while preventing multiple testing problems. The proposed approach can be used for basically all types of forecasts (more precisely for forecasts for all elicitable functionals, that is, types of forecasts for which a consistent scoring function or proper scoring rule exists), for example, univariate or multivariate mean, quantile or probabilistic forecasts.
In a simulation study, we examine the finite-sample performance of our confidence bands. In the absence of temporal dependence, both simultaneous confidence bands (sup-t and Bonferroni) are very close to the nominal coverage even when the vector of skill scores (or expected scores) is large. With increasing temporal dependence, the sup-t bands tend to achieve undercoverage for large vectors of scores. This is due to well-known issues with variance estimation under temporal dependence, which is a notoriously difficult problem. Variance estimators in time series settings usually suffer from a downward bias, which is larger the stronger the temporal dependence and the smaller the sample size is. This is a well-documented and -researched phenomenon for heteroskedasticity and autocorrelation consistent (HAC) variance estimators (for an overview of this literature see lazarus2018), but also for the moving block bootstrap fitzenberger1997. To partially mitigate this undercoverage we propose to use Bonferroni bands in time series settings, whose coverage rates are closer to the nominal level for larger samples due to their asymptotic conservativeness. Both simultaneous bands show huge improvements over the above-mentioned competitors currently used in practice. Pointwise confidence bands not accounting for temporal dependence expectedly show very low coverage rates, which decrease with the strength of temporal dependence and the length of the skill score vector, while pointwise bands accounting for temporal dependence are definitely an improvement over them, but still show drastic undercoverage even for small skill score vectors.
We demonstrate the practical utility of the approach in two case studies. First, we quantify the benefits of time-varying parameters for point and probabilistic macroeconomic forecasts. Macroeconomic forecasts have to adapt to rapidly changing conditions, shocks to the economy and gradual structural change. Thus, time-varying parameter models have been found to be beneficial in this setting (clark2011, dagostino2013, clark2015, knueppel2022). Our confidence bands enable us to quantify those benefits and assess their sampling uncertainty and statistical significance. We compare the forecasting performance of a Bayesian Vector Autoregression (BVAR) with time-varying parameters and stochastic volatility primiceri2005, a state-of-the art model in macroeconomic forecasting, for real GDP growth, inflation and the federal funds rate to the same model with constant parameters. Figure (ref) shows the estimated skill score of multivariate probabilistic forecasts for all three variables, evaluated via the energy score, over eight forecast horizons. The dashed lines depict the four types of confidence bands at 90% confidence level. The skill score is positive for all horizons, starting with an impressive improvement of 15% before gradually declining. The pointwise band under the iid assumption is very narrow, while accounting for temporal dependence roughly doubles its width. Still, the simultaneous bands are considerably wider and show that there is substantial sampling uncertainty in the empirical skill scores, which is dramatically underestimated by the pointwise bands. The simultaneous bands indicate a statistically significant improvement of the time-varying parameter model over the constant-parameter model up to the third quarter. Aggregated over horizons and variables, the estimated skill scores suggest an improvement in forecast accuracy by approximately 1.8% for mean forecasts and by about 6.2% for density forecasts. However, skill and estimation uncertainty vary strongly across type of forecast, forecast horizon and target variable. The uncertainty is considerably higher for mean forecasts, longer forecast horizons and the federal funds rate, illustrating that point estimates should be interpreted with caution.
Our second case study compares weather forecasts from physics-based numerical weather prediction models to data-driven artificial intelligence-based models. The former have been used successfully for decades in operational weather forecasting, while the latter have entered the stage in recent years and shown potential to outperform the former in terms of point forecasting performance (Pathak2022, GraphCast, BiEtAl2023). This leads to the question whether these purely data-driven models can also be leveraged to create better probabilistic forecasts. As a step in this direction, buelte2026 propose and compare different uncertainty quantification (UQ) methods to generate probabilistic weather forecasts from the deterministic Pangu-Weather model BiEtAl2023. We build on their study, examining forecasts from Pangu-Weather combined either with a Gaussian noise perturbation approach (GNP) or with the EasyUQ method introduced by EasyUQ. The benchmark is a traditional physics-based numerical weather prediction model, namely the ensemble forecast of the European Centre for Medium-Range Weather Forecasts (ECMWF). We focus on forecasts of temperature at 2m for the year 2022 over a grid of 35200 points covering Europe and for 31 forecast horizons in 6-hourly steps, and examine the underlying sampling uncertainty. Overall, Pangu-Weather + EasyUQ provides the best forecasts. It outperforms the physics-based ECMWF for lead times up to about 100h significantly at 10%, while the latter is more accurate for horizons greater than 150h. Moreover, there is substantial variation across locations. Both data-driven models achieve the largest improvements in forecast accuracy over the ECMWF in coastal and mountainous areas. At a lead time of 72h, estimation uncertainty is considerably higher over sea than over land grid points. These results demonstrate the need to assess sampling uncertainty in such large-scale forecast comparisons and that our confidence bands provide a valid yet simple approach for this.
The remainder of the paper is organized as follows. In Section (ref), we review consistent scoring functions and proper scoring rules, discuss skill scores, and introduce our multivariate setup. In Section (ref), we introduce our confidence bands, discuss underlying assumptions, establish their validity, and discuss their properties. Section (ref) contains the simulation study, Section (ref) presents the case studies and Section (ref) concludes. The appendix contains proofs and further details regarding assumptions, properties of the confidence bands, simulations, and applications. An R package implementing the simultaneous confidence bands is available at \url{https://github.com/TanjaZahn/UQforecasts}. Replication material is available at \url{https://github.com/TanjaZahn/UQforecasts_replication}.
We denote by $F_W$ the cumulative distribution function (CDF) of a random variable $W$, by $f_W$ its density, by $\mu(W)$ its mean and by $q_{\tau}(W)$ its $\tau$-quantile, $\tau \in (0,1)$. Let $\xrightarrow[]{d}$ and $\xrightarrow[]{p}$ denote convergence in distribution and in probability, respectively. We let $\mathcal{N}(.)$ stand for a (multivariate) normal distribution and let $z_{\tau}$ denote the $\tau$-quantile of a univariate standard normal distribution. Further, we denote vectors and multivariate functions by bold letters, and use uppercase letters for random variables and lowercase letters for their realizations.
Consider a forecaster who wants to forecast a variable of interest $Y_{t}$ from forecast origin $t-h$, where $h$ is the forecast horizon. The forecaster's information set $\mathcal{F}_{t-h}$ is generated by the past of a vector-valued stochastic process $\{\mathbf{Z}_t\}$, which includes $\{Y_t\}$, i.e., $\mathcal{F}_{t-h}=\sigma(\{\mathbf{Z}_s\}_{s \leq t-h})$. The forecaster's goal is to issue a forecast for a certain property of the conditional distribution of $Y_{t}$ given the information set $\mathcal{F}_{t-h}$. Put formally, the forecaster wants to forecast a statistical functional $T$ of the distribution $F_{Y_{t}|\mathcal{F}_{t-h}}$, that is, $T(Y_{t}|\mathcal{F}_{t-h}):=T(F_{Y_{t}|\mathcal{F}_{t-h}})$. Popular choices of the functional $T$ include the mean, $\mu(Y_{t}|\mathcal{F}_{t-h})$, or a certain quantile, $q_{\tau}(Y_{t}|\mathcal{F}_{t-h})$. For probabilistic forecasts, $T(Y_{t}|\mathcal{F}_{t-h})$ becomes $f_{Y_{t}|\mathcal{F}_{t-h}}$ or $F_{Y_{t}|\mathcal{F}_{t-h}}$ itself. We denote the forecast for $T(Y_{t}|\mathcal{F}_{t-h})$ by $X_{t,h}$. Often, we want to compare this forecast to a benchmark forecast denoted by $X^{ben}_{t,h}$. For simplicity, we omit the indices $t$ and $h$ in the following exposition of scoring functions and skill scores.
The key tool underlying relative forecast evaluation are suitable loss or scoring functions $s$, which map forecast-observation pairs $(x,y)$ to the real line. They are called consistent scoring functions in the case of point-valued forecasts and proper scoring rules in the case of probabilistic forecasts gneiting2011,gneiting2007proper. We call loss functions for arbitrary statistical functionals consistent scoring functions throughout the paper, regarding proper scoring rules as a special case pohle2020, fissler2021multi. Consistency is a fundamental property of scoring functions, which ensures that forecasters are incentivized to issue their true beliefs about the respective functional of the variable of interest; for details see gneiting2011. The standard approach to relative evaluation is to rank the forecasts according to their expected score, $\mathrm{E}[s(X,Y)]$, which is the classical measure of forecast accuracy. We define scoring functions such that they are negatively oriented, meaning that lower expected scores indicate more accurate forecasts.
We now briefly review some of the most widely-used scoring functions; for more details and further scoring functions see gneiting2011 and gneiting2007proper. For the sake of this recap, we index the forecast $x_T$ by its target functional $T$. The squared error,
is the classical consistent scoring function for the mean. The Brier score arises as a special case when the variable of interest $y$ is binary and the target functional is the probability $p$ that $y$ is equal to one. The quantile score, $$QS_{\tau} (x_{q_{\tau}},y) = \rho_{\tau} ( y - x_{q_{\tau}} ), \text{ where } \rho_{\tau}(u) = u(\tau - \mathds{1}_{\{u<0\}}),$$ is the most popular consistent scoring function for the $\tau$-quantile. A widely-used scoring function for probabilistic forecasts is the continuous ranked probability score
where $x_F$ is a CDF.
For multivariate forecasts, the target variable is a $D$-dimensional vector $\mathbf{y}$. A consistent scoring function for the mean $\boldsymbol{\mu}$, which is a vector of the same dimension, is the multivariate squared error
where $\lVert \cdot \rVert_2$ denotes the Euclidean norm. It generalizes the univariate squared error and is just the squared Euclidean distance between forecast and observation vector, or equivalently, the sum of all the univariate squared errors. For multivariate probabilistic forecasts, where $\mathbf{x}_{\boldsymbol{F}}$ takes the form of a multivariate CDF, the energy score, a multivariate generalization of the CRPS, is a proper scoring rule:
where $\mathbf{D}_{\mathbf{x}_{\boldsymbol{F}}}$ and $\mathbf{D}^*_{\mathbf{x}_{\boldsymbol{F}}}$ are two independent draws from the forecast distribution $\mathbf{x}_{\boldsymbol{F}}$.
In many forecasting problems we want to compare forecasts for multiple horizons, variables and/or multiple locations as laid out in Section (ref). In these settings, it may be useful to aggregate scores in one or several dimensions to get a summarized assessment, before studying more disaggregated results. For example, our goal is to measure the aggregated accuracy of forecasts $x_{i,T}$ for a functional $T$ for each element $y_i$ from a vector $\mathbf{y}=(y_1, y_2,...,y_I)$, where the index $i$ possibly runs through multiple horizons, variables and/or locations. If the scoring function $s$ is consistent for $T$, any linear combination of the single scores is a consistent scoring function for the $I$-dimensional vector of the functionals $\boldsymbol{T} = (T,...,T)^\prime$, see, e.g., dawid2014 and pic2025. In our case studies, we use an equal-weighted aggregation, that is, the sum of the individual scores,
A skill score measures the relative improvement in terms of the expected score of a forecast $X$ over a benchmark forecast $X^{ben}$. We make the following assumption to ensure the existence of skill scores.
This assumption is usually mild and is satisfied for all scoring functions discussed in Section (ref) as long as the forecaster does not have perfect foresight; see Appendix (ref) for a more detailed discussion on the assumptions underlying the usage of skill scores. It excludes the use of the popular log score for probabilistic forecasts, for which it does not make sense to compute skill scores. When our confidence bands are constructed for expected scores or their differences, Assumption (ref) is not required. Thus, in those cases our bands provide uncertainty quantification for the log score as well.
Thus, it is simply the relative change in accuracy as measured by the expected score over the benchmark. A positive skill score indicates that the forecast improves over the benchmark forecast, while a negative skill indicates a worse performance than the benchmark. The skill is zero if both forecasts are equally accurate. The skill score is at most one, but the upper bound is usually not attained, see Appendix (ref) for more details.
Skill scores lead to the same ranking as the underlying expected scores and thus incentivize the forecaster correctly as well. What they add compared to expected scores or their differences, $\mathrm{E} [ s(X^{ben},Y) ] - \mathrm{E} \left[ s \left( X,Y \right) \right],$ is interpretability: The value of an expected score or the difference in expected scores is not interpretable per se, whereas a skill score directly expresses the relative improvement in forecast accuracy. Whether the improvement is practically relevant or negligible, of course, might depend strongly on the application at hand. The desire for better interpretability is certainly one of the reasons why relative mean squared errors are often used in the evaluation of mean forecasts. The analogous quantity can also be computed for other types of forecasts and this is what we call relative accuracy. Skill score and relative accuracy carry the same information and it is just a matter of taste, which one is used. We use the skill score throughout the rest of the paper, but our approach to uncertainty quantification works in exactly the same way with relative accuracy.
Now, we make the multivariate setup explicit. We consider $M$ forecasting methods indexed by $m = 1,\ldots,M$, $H$ forecast horizons indexed by $h = 1,\ldots,H$, $D$ target variables indexed by $d = 1,\ldots,D$, and a generic multi-index $(i_1,\ldots,i_K)$ with $i_k = 1,\ldots,I_k$ for $k=1,\ldots,K$, which may represent, for example, spatial locations on a two-dimensional grid. The forecasts $X_{t, i_1,...,i_K, d, h, m}$ are then indexed by time, the generic index, the variables, the forecast horizons, and the methods. Thus, at every point in time $t$, forecasts and observations are arrays,
and
leading to an array of scores at time $t$,
where $S_{t, i_1,...,i_K, d, h, m} := s(X_{t, i_1,...,i_K, d, h, m},Y_{t, i_1,...,i_K, d})$.
The expected scores come in this form as well:
In the shorthand notation of this section, we write the skill score of method $m_{1}$ relative to benchmark method $m_2$ as
Consider the set $\mathcal{M} \subset \{ 1,...,M \} \times \{ 1,...,M \}$ containing all combinations of methods $m_1$ and $m_2$, for which we want to consider skill scores in our problem at hand and let the index $u=1,...,U=|\mathcal{M}|$ run through all those combinations. Often, all methods are compared to a single benchmark method from $ \{ 1,...,M \} $ such that $|\mathcal{M}|=M-1$. The corresponding array of skill scores,
can then be computed from the array of expected scores via (ref).
For our theoretical discussion of confidence bands, the array structure of scores, expected scores and skill scores does not play a role as we treat all dimensions in the same way. For this purpose, we vectorize the arrays. Letting an index $p=1,...,P$ with $P=I_1\cdot...\cdot I_k D H M$ run through $i_1,...,i_K$, $d$, $h$ and $m$, the arrays of scores and expected scores from (ref) and (ref) become vectors of length $P$,
The array of skill scores from (ref) becomes a vector of length $J= I_1\cdot...\cdot I_K D H U$,
which arises from $\mathrm{E}[\mathbf{S}_t]$ via a differentiable function $\mathbf{g}: \mathbb{R}^{P} \rightarrow \mathbb{R}^{J}$ defined via the formula for skill scores in terms of expected scores (ref) and the choice of the set $\mathcal{M}$, i.e.,
Our goal is to construct confidence bands for this vector of skill scores.
In practice, we observe an evaluation sample of size $N$, where at every time point $t=1,...,N$ the forecast and observation arrays from (ref) and (ref) are observed. Thus, the evaluation sample comes in the form of arrays as well by adding the time index,
and
from which the array of scores can be computed:
Averaging over time leads to an array of average scores, which is the empirical analog of the array of expected scores from (ref),
with $$ \quad \overline{S}_{i_1,...,i_K,d, h,m} = \frac{1}{N} \sum_{t=1}^{N} S_{t,i_1,...,i_K,d,h,m} .$$
From this, the empirical skill scores can be computed in the same way as the population skill scores from the expected scores via
which leads to an array of empirical skill scores, the empirical analog of (ref),
Again, we consider the vectorized versions. The sample of scores from (ref) becomes a matrix,
which we rather write as
in the following to emphasize that it is a sample of score vectors. The arrays of average scores from (ref) and the empirical skill scores from (ref) become vectors
which are the empirical analogs of (ref) and (ref). Similar to the theoretical analogs, see (ref), the empirical skill scores arise from the average scores via the differentiable function $g$, $\mathbf{\widehat{SS}} = \mathbf{g} \left( \mathbf{\overline{S}}\right)$.
Under our Assumption (ref) below, a law of large numbers holds for $\{\mathbf{S}_t\}$, that is, the average scores $\overline{\mathbf{S}}$ are consistent estimators for the expected scores, $\mathbf{\overline{S}} \overset{p}{\to} \mathrm{E}[\mathbf{S}_t]$. Under consistency of the average scores, the continuous mapping theorem implies because of (ref) that the vector of empirical skill scores is a consistent estimator of its population counterpart, $\mathbf{\widehat{SS}} \overset{p}{\to} \mathbf{SS}$.
Our construction of confidence bands builds on the scores $\mathbf{S}_t$ from (ref) themselves. To construct the bands, we need the sample of scores $\{\mathbf{S}_t\}_{t=1}^N$ from (ref). we make our assumptions directly on the stochastic process of scores $\{\mathbf{S}_t\}_{t \in \mathbb{Z}}$ and not on the processes of forecasts and observations. This has the advantage that all types of forecasts, e.g., mean, quantile and probabilistic forecasts, can be handled within a single, versatile framework. In this respect, we follow diebold1995, who essentially assume that the score differences of the two forecasting methods they want to compare follow a univariate central limit theorem; see also diebold2015 for an insightful discussion of the classical Diebold-Mariano assumptions. We aim for a multivariate version of those assumptions. Our key assumption is that a multivariate central limit theorem holds for the scores $\{\mathbf{S}_t\}_{t \in \mathbb{Z}}$.
We further assume that there exists a consistent estimator for $\mathbf{\Omega}$.
Assumptions (ref) and (ref) represent a suitable multivariate extension of the classical Diebold-Mariano assumptions. Without limiting normality and the possibility to estimate the corresponding covariance matrix, inference on expected scores and skill scores will hardly be possible. Thus, we also regard them as minimal assumptions that are needed in this context. While we state them at a high level to remain agnostic about the precise dependence structure, Appendix (ref) discusses classical sufficient conditions on the dependence structure, homogeneity and moments of the process $\{\mathbf{S}_t\}_{t \in \mathbb{Z}}$. Depending on the problem at hand, different sets of sufficient conditions may be more suitable and stating the high-level assumptions makes it possible to choose such a set of conditions tailored to it. One could also discuss sufficient conditions on the observations and forecasts themselves that imply our assumptions on the scores, see again diebold2015 for a discussion of this in the univariate setting.
As skill scores arise from expected scores via a differentiable function $\boldsymbol{g}$, see (ref), invoking the delta method, the normality assumption on the average scores (Assumption (ref)) implies that the empirical skill scores are asymptotically normal as well,
with
where $\nabla \boldsymbol{g}$ denotes the Jacobian of $\boldsymbol{g}$. We lastly assume that the diagonal elements of $\boldsymbol{\Sigma}$ are positive, which essentially just means that there is sampling variability in the skill scores. Otherwise, inference on the skill scores would not be of interest anyway.
Consider an asymptotic confidence interval of level $1-\alpha$ for a single skill score,
where $\widehat{\sigma}_{j}$ is a consistent estimator of the standard deviation of $\widehat{SS}_{j}$, $\sigma_{j} = \sqrt{\boldsymbol{\Sigma}_{jj}/{N}}$. By construction, this interval fulfills the defining condition for a confidence interval of level $1-\alpha$, $\lim_{N \to \infty} P \left( {SS}_{j} \in \widehat{CI}^{1-\alpha} (SS_{j}) \right) \geq 1 - \alpha$, with equality. If such a condition is fulfilled, we say that a confidence interval or band has correct asymptotic coverage. If the condition is fulfilled with equality, we say that it has exact asymptotic coverage.
When we are interested in a vector (or array) of skill scores $\mathbf{SS}$ and want to quantify the sampling uncertainty of its estimator $\mathbf{\widehat{SS}}$, it is natural to consider simultaneous confidence bands, which contain the whole parameter vector with a probability of at least $1-\alpha$.
The advantage of such a confidence band compared to other confidence regions is that it is rectangular and thus can be easily visualized no matter what the length of the parameter vector $J$ is, whereas, e.g., confidence ellipsoids cannot be visualized in dimensions higher than 2.
Just using a pointwise confidence band, i.e., the Cartesian product of the confidence intervals from (ref), which amounts to setting $\widehat{c}:=z_{1-\alpha/2}$, leads to a coverage that is unknown and usually much smaller than $1-\alpha$ for $J>1$, especially for large $J$; see Figure (ref) and the discussion in Appendix (ref) and the simulations in Section (ref).
The simultaneous confidence band replaces the standard normal quantiles $z_{1-\alpha/2}$ of the pointwise band with a larger scaling factor $\widehat{c}$. There are different approaches to constructing simultaneous confidence bands, that is, to determine $\widehat{c}$, which montielolea2019 review and compare. We focus on Bonferroni and sup-t bands.
To understand the construction of those bands (and where the name sup-t bands comes from), we can express the joint coverage probability from (ref) in terms of the distribution of the maximum absolute values of individual t-statistics:
Combining this representation with our assumptions (ref) to (ref), we have that (see Lemma (ref) in the appendix)
where again $\bs V = \left( V_1, \cdots , V_J \right)' \sim \mathcal{N}(\bs 0,\bs \Sigma)$ and $c$ is defined as the probability limit of $\widehat{c}$, $\widehat{c} \xrightarrow[]{p} c$. This explains the choice of $\widehat{c}$ in the sup-t band as the $1-\alpha$ quantile of this distribution of the maximum of the absolute value of correlated standard normal random variables, $q_{ \bs \Sigma, 1-\alpha}$, or its estimator $\widehat{q}_{\bs \Sigma,1-\alpha}$, respectively, and shows that the sup-t band has exact asymptotic coverage, i.e., (ref) holds with equality. $q_{\bs \Sigma, 1-\alpha}$ is also called the $1-\alpha$ equicoordinate quantile genz2025 of the multivariate normal distribution with covariance matrix $\bs \Sigma$ as it fulfills $ P \left( | \bs \Sigma_{11}^{-1/2} V_1| \leq q_{\bs \Sigma, 1-\alpha}, \ \dots \ , | \bs \Sigma_{JJ}^{-1/2} V_J| \leq q_{\bs \Sigma, 1-\alpha} \right) = 1 - \alpha \, $.
For Bonferroni bands we simply choose the scaling factor $\widehat{c}$ as yet another quantile of the standard normal distribution, namely $z_{1-\alpha/(2J)}$. Here the scaling factor is a constant and requires no estimation of the covariance matrix. It suffices to estimate the diagonal elements $\boldsymbol{\Sigma}_{jj}/N$ of the asymptotic covariance matrix. Since $z_{1-\alpha/(2J)} \geq q_{\bs \Sigma, 1-\alpha}$ (see, e.g., Lemma 1 in hassler2025), of which $\widehat{q}_{\bs \Sigma,1-\alpha}$ is a consistent estimator, Bonferroni bands are asymptotically wider than sup-t bands. Thus, they have correct asymptotic coverage (ref), but tend to be conservative.
In Appendix (ref), we compare the width and coverage properties of pointwise, sup-t, and Bonferroni bands from a large-sample perspective, i.e.\ assuming that $N \to \infty$, which enables us to use the population covariance matrix to calculate the equicoordinate quantiles and the formula for the asymptotic coverage probability from (ref). We illustrate there, in particular in figures (ref) and (ref), that the width of the simultaneous bands only grows slowly with $J$. Thus, our simultaneous bands are not prohibitively wide and still informative for large vectors of skill scores, while the pointwise bands suffer from serious undercoverage already for small $J$ and are thus virtually useless. Concerning the relation of Bonferroni and sup-t bands, we find that the differences in width and coverage are negligible when the elements of the empirical skill score vector are independent or mildly correlated, while they become more pronounced when the correlation gets stronger. From an asymptotic perspective, the Bonferroni band is thus less desirable than the sup-t band.
These findings on the relative width of Bonferroni and sup-t bands are consistent with the finite-sample performance in the simulation study from Section (ref). However, in our asymptotic comparisons, we set aside the issue of variance estimation. It is a well-known problem that variance estimation under temporal dependence is plagued by a downward bias, which gets stronger with the degree of temporal dependence and only disappears in large samples, leading to oversized tests and undercoverage of confidence intervals. This issue is documented in the large literature on heteroscedasticity and autocorrelation consistent (HAC) variance estimation (see lazarus2018 for an overview), but also for the moving block bootstrap fitzenberger1997, which we use (for details see Section (ref)). Thus, our confidence bands are too narrow as well under temporal dependence. As a partial remedy, we recommend the Bonferroni bands as a default choice in time series settings. They show less severe undercoverage in our simulations because their conservativeness under strong dependence works against the undercoverage that comes from the variance estimation problem.
We use a bootstrap algorithm to construct our simultaneous confidence bands, that is, to estimate $\sigma_{j}$ and for the sup-t bands also $q_{\bs \Sigma, 1-\alpha}$. The initial step of the algorithm (step 2 in Algorithm (ref) below) is generating bootstrap resamples of scores $\{\boldsymbol{S}_t\}_{t=1}^N$. For this step, any type of bootstrap can be chosen that is suitable for the type of data at hand. This suitability is formulated in the following assumption, which guarantees that the bootstrap reproduces the asymptotic distribution of $\boldsymbol{\overline{S}}$.
Usually, the existence of a central limit theorem already essentially guarantees that there is a valid bootstrap algorithm lahiri2003, that is, Assumption (ref) is not really a strong additional assumption on top of Assumption (ref). The possibility of flexibly choosing the type of bootstrap underlying the algorithm makes our bootstrap bands very versatile as it can be adapted by the user to different types of data. As we are concerned with time series data in our applications, we use the moving block bootstrap discussed below but other valid bootstrap procedures can be used analogously. For example, for data, for which the iid assumption is reasonable, a classical iid bootstrap can be used.
The full bootstrap algorithm is as follows.
After generating $B$ bootstrap samples of scores, Algorithm 1 estimates the empirical standard deviations of the resulting bootstrap skill scores. For the Bonferri bands, this is the only quantity that needs to be estimated, so in this case the algorithm amounts to a bootstrap variance estimator since the only thing left is to set the scaling factor to $z_{1-\alpha/(2J)}$ and compute the bands. To construct sup-t bands, the scaling factor is determined as the $1-\alpha$ quantile of the bootstrap sample of the maximum of the absolute t-statistics as proposed by montielolea2019. We establish the validity of the bands in Proposition (ref).
Of course, as an initial step in the algorithm, the arrays of scores from (ref) have to be vectorized to arrive at $\left\{ \boldsymbol{S}_t \right\}_{t=1}^N$. As a final step, the process has to be reversed to arrive at an array of bands, more precisely, one array of upper and one of lower bounds of the same nature as the array of empirical skill scores from (ref). Those three together is the final output of the algorithm:
We use the moving block bootstrap by kunsch1989 in step 2 of Algorithm (ref), which replicates the temporal dependence from the original series of scores by resampling blocks of length $l$ and stringing them together. The block length $l$ has to grow with the sample size, $l \rightarrow \infty$ as $N \rightarrow \infty$, such that asymptotically the dependence structure of the original series is fully replicated. At the same time, it has to grow slower than the sample size itself, $\frac l N \rightarrow 0$, such that the number of blocks goes to infinity as well. Thus, the block length needs to fulfil $\frac l N + \frac 1 l \rightarrow 0$. Apart from that, the moving block bootstrap requires essentially that a central limit theorem holds to be valid. We spell out sufficient conditions in Appendix (ref). To generate bootstrap resamples of size $N$, the moving block bootstrap draws $\left \lceil \frac{N}{l} \right \rceil$ blocks with replacement from the set of $N-l+1$ blocks $\{\mathcal{B}_1,...,\mathcal{B}_{N-l+1} \}$, where $\mathcal{B}_i=(\boldsymbol{S}_i,...,\boldsymbol{S}_{i+l-1})$, and strings them together (and finally discards the last $\left \lceil \frac{N}{l} \right \rceil - N$ observations) to generate a bootstrap sample $\left\{ \boldsymbol{S}^*_t \right\}_{t=1}^N$. To replicate the cross-sectional dependence, we draw the entire vectors $\boldsymbol{S}_i$, keeping the original structure of the elements. For an in-depth discussion of the moving block bootstrap see lahiri2003. The moving block bootstrap requires the choice of a block length $l$. The theoretically optimal $l$ (in the sense of minimizing the mean squared error) is of order $N^{\frac 1 4}$ lahiri2003. Data-driven choices of the block length are difficult as the optimal block length depends on the distribution of parameters of a higher order than the one which is bootstrapped. Therefore, one usually needs to apply a second bootstrap after the first one, for which one wants to determine the block length in the first place. Consequently, $l$ is usually chosen ad hoc in practice. Nevertheless, in the univariate case such procedures exist lahiri2003, while we are not aware of any in a multivariate setting. Based on our simulation results and the surrounding discussion in Section (ref), we recommend choosing a block length in the order of magnitude of $l = 3 \lfloor N^{1/4} \rfloor $ and checking robustness against some alternative choices. When $l=1$, the moving block bootstrap reduces to a classical iid bootstrap.
The applicability of our simultaneous confidence bands is not restricted to skill scores $\boldsymbol{SS}$, but can be used in the same way for expected scores $\mathrm{E} [\boldsymbol{S}_t]$ or relative accuracy, $\boldsymbol{RA} = (1-SS_1,...,1-SS_J)^\prime$ (see Definition (ref)), by replacing empirical skill scores $\mathbf{\widehat{SS}}$ and their bootstrapped versions $\widehat{\boldsymbol{SS}}^{*,b}$ with $\mathbf{\overline{S}}$ and $\mathbf{\overline{S}}^{*,b}$ or $\mathbf{\widehat{RA}}$ and $\widehat{\boldsymbol{RA}}^{*,b}$, respectively.
To assess the finite-sample performance of our confidence bands, we simulate directly from a score process, following hansen2011 and quaedvlieg2021. To be able to control the strength of temporal dependence as well as the cross-correlation between the scores separately and by a single parameter each, we simulate the $P$-dimensional vector of scores from (ref) by a simple vector autoregressive model of order 1 (VAR(1) model) similar to knueppel2022:
where $\mathbf{I}$ denotes the identity matrix and $\mathbf{A} = a \mathbf{I} $ is a diagonal matrix with diagonal element $a$. The error term follows a multivariate normal distribution $\bs \varepsilon_t \sim \mathcal{N} ( \mathbf{0} , \mathbf{V} )$, where the covariance matrix has an equicorrelation structure. The variances are 1 and the covariances/correlations are all equal to $v$, $\mathbf{V} = v \mathbf{J} + (1-v) \mathbf{I}$, where $\mathbf{J}$ denotes a matrix of ones. Thus, the parameter $a$ governs the strength of temporal dependence, while $v$ controls the cross-correlation. We choose the following parameter values: $a=0,0.3,0.6$, $v=0,0.3,0.6$, and $P=2,5,25, 100, 400$. We set all elements of $\mathrm{E}[ \mathbf{S}_t ] $ to 10. The value of the expected scores should not matter as we are rather interested in the uncertainty surrounding it. However, the expected score of the benchmark method should be sufficiently far away from 0 so that the average score is not 0 and empirical skill scores can be calculated. For each parameter combination, we simulate 1000 $P$-dimensional time series of scores of length $N=100,400$, calculate the $J=P-1$ empirical skill scores relative to the $P$th score and construct $90 \%$ confidence bands employing a moving block bootstrap with block length $l=q \lfloor N^{1/4} \rfloor $ for $q=1,2,3$. We also compare results to an iid bootstrap, where $l=1$, which is appropriate for iid data ($a=0$). We also include pointwise bands (with both variants of the bootstrap) as an important competitor widely used in practice. We report the fraction of times the true skill score vector (which is $\mathbf{0}$) falls into the bands. All tabulated simulation results can be found in Appendix (ref).
First, we focus on the results for small to medium dimensions of the score vector, i.e., $P=2,5,25$, and compare the iid bootstrap with block length $l = 1$ to a block bootstrap with $ l = 3 \lfloor N^{1/4} \rfloor $. Table (ref) displays the empirical coverage of the confidence bands for the skill score vector when expected scores are independent over time, i.e., $a = 0$. As expected, the pointwise bands show a severe undercoverage for a skill score vector containing more than one element ($ P > 2$) for both choices of the block length, which becomes more pronounced for larger values of $P$ and smaller sample sizes $N$. In contrast, both simultaneous bands using the iid bootstrap are very close to the nominal coverage of 0.9. For large values of $P$, the sup-t bands show a slight undercoverage for $N = 100$ while the Bonferroni bands show a slight overcoverage for $N = 400$. Using a block length of $ l = 3 \lfloor N^{1/4} \rfloor $, which is quite large when there is no temporal dependence, decreases the coverage of the bands only to some degree for $N = 100$, while it barely affects results for $N = 400$. The level of cross-sectional dependence, which is controlled by $v$, does not alter the results considerably, suggesting that the bootstrap method correctly accounts for it.
Table (ref) displays the empirical coverage rates when the scores are autocorrelated ($a > 0$). Unsurprisingly, using a block bootstrap to capture the persistence yields better results than using the iid bootstrap, especially for high levels of $a$. Again, both simultaneous bands outperform the pointwise bands by a large margin for $P > 2$. Although the sup-t bands have exact asymptotic coverage by construction, the finite sample results reveal that their empirical coverage is lower than the nominal level of 0.9, which is more pronounced for smaller sample sizes $N$, larger dimensions $P$ of the score vector, and stronger levels of persistence $a$. This is as expected given the downward bias in variance estimators under temporal dependence discussed in Section (ref), which also plagues the moving block bootstrap fitzenberger1997. This issue naturally exponentiates with longer skill score vectors, i.e., with higher $P$ or $J$, respectively. As illustrated in Appendix (ref), the Bonferroni bands are very similar to the sup-t bands in terms of asymptotic width and coverage for weak levels of dependence. With increasing dependence, the Bonferroni bands suffer from overcoverage asymptotically. We therefore recommend the Bonferroni bands to mitigate the notorious undercoverage in finite samples under temporal dependence. In our simulations, the empirical coverage of the Bonferroni bands is better than the coverage of sup-t for $ P \geq 5 $ and it is close to 0.9 when $N = 400$. As before, the cross-sectional dependence seems to be captured well by the bootstrap as it does not influence the coverage rates considerably. The results are qualitatively similar when analyzing confidence bands for the expected scores themselves (see Tables (ref) and (ref)).
In Tables (ref) and (ref), we consider different choices of the block length $l = q \lfloor N^{1/4} \rfloor $ with $q=1,2,3$. In general, the coverage rates are quite robust to the choice of $l$, but, of course, the optimal block length tends to increase with the degree of temporal dependence. We suggest choosing $q = 3$ as a default for the following reasons. First, it is reasonable to assume that scores are temporally dependent to some degree in most applications. In our simulation study, selecting $q = 3$ frequently achieves the best coverage for $ N = 400$ when $ a > 0$. Second, the differences in coverage rates due to the block length are more severe for stronger levels of temporal dependence while they are less relevant for small $a$, which supports choosing a larger block length as default.
For our asymptotics, we consider the length of the score vector fixed while the sample size goes off to infinity. In order to check whether the empirical results are still decent for very large vectors of scores, we also report results for $P = 100$ and $P = 400$ in Tables (ref) and (ref). Even when $ P = 400$, the Bonferroni bands work very well for skill scores when the sample size $N$ is 400, especially for $a = 0$ and $a = 0.3$. Although the coverage becomes a bit worse for $a = 0.6$, it still about 80% even for the very large skill score vectors. In the small sample, results are satisfying for $a = 0$, depending on the block length. However, with increasing temporal dependence, the coverage rates become considerably smaller than 0.9 in the high-dimensional setting. Thus, with very large vectors of skill scores under strong temporal dependence one should be aware that even the simultaneous bands may suffer from serious undercoverage.
Bayesian Vector Autoregressions (BVAR) are workhorse models in macroeconomic forecasting. In this application, we study the benefits of including time-varying parameters and stochastic volatility in those models with respect to their forecasting performance. This question has been addressed, for example, by clark2011, dagostino2013, clark2015, and knueppel2022. We are able to quantify those benefits using skill scores, characterize the surrounding uncertainty, and assess their statistical significance via our confidence bands.
Using data from FRED-QD mccracken2020, we generate mean and probabilistic forecasts for three quarterly macroeconomic variables: real GDP growth, inflation, and the federal funds rate. We consider forecast horizons up to eight quarters ahead, i.e., $h=1,...,8$. Our evaluation sample runs from 1991:Q1 to 2019:Q4 and consists of $N=116$ observations. We compare the performance of two forecasting methods, that is, we analyze the skill score of method $m_1$ relative to the benchmark method $m_2$. The forecasting method of interest, $m_1$, is a BVAR that allows for time-varying parameters and stochastic volatility. We choose the model specification introduced by primiceri2005 with two lags:
where $\mathbf{Y}_t = (Y_{t,1}, Y_{t,2}, Y_{t,3})' $ denotes the vector of variables at time $t$ consisting of real GDP growth, inflation, and the federal funds rate. The vector of intercepts $\mathbf{c}_t$, the coefficient matrices, $\mathbf{B}_{1,t}$ and $\mathbf{B}_{2,t}$, as well as the error covariance matrix $\mathbf{\Psi}_t$ are allowed to vary over time. In short, they evolve as (geometric) random walks. Priors are chosen analogously to primiceri2005. For a more in-depth discussion of the model see primiceri2005 and delnegro2015. Our benchmark method $m_2$ is the same model as in (ref), but the parameters $\mathbf{c}$, $\mathbf{B}_{1}$, $\mathbf{B}_{2}$ and $\mathbf{\Psi}$ are assumed to be constant over time. Altering the respective priors `switches off' the time-varying elements. We estimate both models on rolling windows of 120 observations, from which we use 48 observations to determine prior parameters via least squares (LS). The main Bayesian estimation is performed on the remaining 72 observations. The number of draws in the MCMC algorithm is 50,000, from which we keep every tenth draw. We use the R package bvarsv krueger2015 for implementation.
Forecast evaluation is carried out via univariate and multivariate scoring functions. Thus, variables are treated separately or jointly. For mean forecasts, we use the univariate squared error and its multivariate version as defined in equations (ref) and (ref). For density forecasts, we use the CRPS and the energy score, see equations (ref) and (ref). In both applications we choose a confidence level of 90% and use Bonferroni bands with a block length of $ l = 3 \lfloor N^{1/4} \rfloor $ in the main analysis.
First, we aggregate the skill score over all dimensions, i.e., over forecast horizons and over variables. This allows us to assess the gains of time-varying parameters based on a single number for each type of forecast. The aggregation is as follows. First, we use the multivariate scoring functions to evaluate the forecasts for the variables jointly. Then, we apply the equal-weighting scheme from (ref) over forecast horizons. Figure (ref) depicts the results for mean forecasts in the upper panel and for density forecasts in the lower panel. The blue point is the estimated skill score and the blue vertical line depicts the respective confidence interval. The estimated skill score is positive in both panels with values of 1.79% and 6.2%. These figures suggest that time-varying parameters improve both types of forecasts, i.e., mean and density forecasts. However, the estimate for mean forecasts is subject to larger sampling uncertainty as the width of the confidence bands is 0.325. While the upper bound suggests an improvement in forecasting performance by 18% due to time-varying parameters, the lower bound indicates a decline in performance by 14.4%. Of course, the skill score is not significantly different from 0 at the 10% level since the bands include the null. A traditional hypothesis test alone would only indicate that the null of equal predictive accuracy is not rejected at the 10% level; the confidence band additionally reveals that the data are consistent with improvements as large as 18% or deteriorations as large as 14.4%, illustrating the practical value of uncertainty quantification beyond binary test decisions. For the density forecasts, the confidence bands are much narrower with a width of 0.17. Thus, estimation uncertainty is much lower. As before, the hypothesized value of 0 is included in the confidence band, i.e., the traditional test of equal predictive ability is not rejected at 10% significance level. However, the confidence bands show that the lower bound is at -2.28%, suggesting rather small losses, while the upper bound is at 14.7%.
Although the fully aggregated skill score allows us to summarize information, reduce dimensionality and is therefore a good starting point for every forecast comparison, it might conceal important patterns in the forecasting performance across forecast horizons and variables. Thus, we next disaggregate skill scores over one dimension while keeping the other dimension aggregated. On the left side of Figure (ref), we disaggregate over forecast horizons. The solid and dashed lines depict the estimated skill score and its confidence bands, respectively. The latter are simultaneous with respect to all eight forecast horizons. The estimated skill score is positive at short forecast horizons in both panels. For mean forecasts, estimation uncertainty is again quite high. For example, at a forecast horizon of one quarter, they span from 0.34% to 26.5% around a point estimate of 13.44%. With increasing forecast horizon, the estimated skill score tends to depreciate while uncertainty rises further. However, the lower bound remains above zero up to a forecast horizon of two quarters. The confidence bands for the density forecasts, which we already took a look at in Figure (ref) in the introduction, are much more concentrated. At the first two horizons the lower bounds are approximately 5.56% and 4.85%, while the upper bounds reach 23.8% and 24.2%, respectively, with point estimates of approximately 15% suggesting considerable gains in the forecasting performance due to the inclusion of time-varying parameters in the model. The right side of Figure (ref) displays the skill score disaggregated over variables by using the univariate scoring functions instead of the multivariate analogs. The confidence bands are simultaneous with respect to all three variables. First, they reveal a strong variation in uncertainty across variables despite the forecasts being generated by the same model. In both panels, the confidence bands for the federal funds rate are more than twice as wide as the bands for real GDP growth. Second, the upper bound is always larger in magnitude than the lower bound in all cases. This observation is more pronounced for the density forecasts. For example, consider the point estimate of 4.44% for real GDP growth. The lower bound for this estimate is only at -1.25% while the upper bound is at 10.1%. Thus, gains due to time-varying parameters could potentially be quite large while potential losses seem limited. As before, uncertainty is higher for mean forecasts than for density forecasts.
Finally, we analyze forecast performance disaggregated over both dimensions in Figure (ref). Now, the confidence bands are not only simultaneous within each graph but across all three graphs contained in the upper and lower panel, respectively. The mean forecasts for real GDP growth and the federal funds rate have a positive estimated skill score up to a forecast horizon of seven and five quarters. The estimate for inflation is also positive except for $h=2,3$. However, the confidence bands are rather wide in all cases. For example, they range from -12% to 22.7% around a point estimate of 5.4% for real GDP growth at $h=1$. As before, the bands are widest for the federal funds rate. However, they are above zero for short forecast horizons suggesting gains in the forecasting performance due to time-varying parameters. At $h=1$, the lower and upper bounds are 11.1% and 53.2%, respectively, while the point estimate is 32.2%. A similar trajectory of the estimated skill score emerges for density forecasts but estimation uncertainty is much lower. For example, the bands for real GDP growth at $h=1$ span from -3.93% to 13.7% around a point estimate of 4.86%. Similar to before, the confidence bands for the federal funds are well above zero up to a forecast horizon of two quarters. At $h=1$, they range from 16.3% to 50.2% with an estimated skill score of 33.2% suggesting a considerable improvement in the forecasting performance.
In Appendix (ref) we report additional details and results for this case study. In Table (ref) and the surrounding text we analyze the average width of the Bonferroni bands presented in Figures (ref) to (ref), which confirms the impression that sampling uncertainty for mean forecast skill is considerably higher here than for density forecast skill. While we focus on the Bonferroni band in the case studies since it is our recommended choice, we also present the competitors in Figure (ref) in the introduction and in Figures (ref) to (ref). Again, the pointwise bands assuming independent observations are much narrower than the pointwise bands constructed via the block bootstrap, which are in turn substantially narrower than the simultaneous bands. The sup-t bands are slightly narrower than the Bonferroni bands, which is also confirmed by the average widths in Tables (ref) and (ref).
Summarizing our main findings, the time-varying parameter model can lead to substantial gains in forecasting performance. Performance gains are higher and sampling uncertainty is lower for density forecasts than for mean forecasts and for the federal funds rate compared to inflation and GDP growth. Furthermore, sampling uncertainty tends to grow with the forecast horizon.
Physics-based numerical weather prediction (NWP) models have long been the standard for probabilistic weather forecasting, generating ensembles by running simulations with perturbed initial conditions. Recently, purely data-driven artificial intelligence (AI)-based models for weather forecasting have emerged as competitive alternatives ECMWF2023rise, with noteworthy examples including FourCastNet Pathak2022, Pangu-Weather BiEtAl2023, and GraphCast GraphCast. Unlike NWP systems, data-driven models focus on predicting future weather states from initial conditions, leveraging statistical relationships learned from historical data. In addition to improved forecasts, key advantages of these models include significantly reduced computational costs, lower energy consumption, and faster forecast generation once the model has been trained. Systematic and case study-based comparisons of physics- and AI-based weather models have been a recent focus of research interest, including the development of tailored benchmarking frameworks such as WeatherBench 2 WB2 and novel evaluation metrics gneiting_etal_2026_probabilistic.
A major limitation of many current data-driven weather models is that they solely provide point forecasts. To address this limitation, buelte2026 propose and compare different uncertainty quantification (UQ) methods to generate probabilistic weather forecasts from the deterministic Pangu-Weather model. Here, we follow their setup and compare two of the approaches. In the Gaussian noise perturbation (GNP) approach, an ensemble of initial conditions is generated by adding independent Gaussian noise at every grid point, and the deterministic Pangu-Weather model is started from those initial conditions. Further, we consider a post-hoc UQ method, where a statistical model is learned based on a training dataset of past forecasts and observations to turn the point forecasts into probabilistic ones by supplementing them with uncertainty information. Specifically, we utilize the EasyUQ method proposed by EasyUQ, which is based on isotonic distributional regression IDR. EasyUQ offers appealing theoretical properties and does not require any choices of tuning parameters. It is a straightforward and readily applicable UQ method to generate probabilistic forecasts. EasyUQ is applied separately for every grid point and location, see BuelteEtAl2024 for details.
In the following, we focus on probabilistic forecasts of temperature at 2m for the calendar year 2022 and choose the CRPS from (ref) as the score. Probabilistic forecasts, and thus score values, for all models are available on a two-dimensional grid over Europe with a resolution of 0.25$^\circ$, which corresponds to a total of 35200 grid point locations, and for 31 forecast horizons in 6-hourly steps (i.e., up to a maximum lead time of 186 hours). We use the operation ensemble forecasts of the European Centre for Medium-Range Weather Forecasts (ECMWF) as a physics-based benchmark model when computing skill scores. The ECMWF ensemble prediction system is based on the ECMWF Integrated Forecasting System and is a standard benchmark in weather forecasting research.
Figure (ref) shows the average score, i.e., the estimated expected CRPS, and the estimated skill score, which we call CRPSS, over the test set, along with 90% Bonferroni confidence bands. The Pangu-Weather + GNP approach produces the least skillful forecasts overall, and is only competitive with the ECMWF ensemble (without being significantly better) within the first 24 hours of lead time. The best overall forecasts are provided by the Pangu-Weather + EasyUQ method. Disaggregating over forecast horizons demonstrates that the data-driven model provides better forecasts for lead times of up to around 100h, whereas the physics-based ECMWF ensemble provides better forecasts for lead times exceeding around 150h, with no significant differences in between. We feel that the statement about the significance requires a word of caution regarding the interpretation of overlapping confidence bands: There is a common misconception that if two confidence intervals or bands of level $1-\alpha$ overlap, the hypothesis that the two corresponding parameter vectors are equal cannot be rejected at significance level $\alpha$. This is not true schenker2001. What is true is that if the bands do not overlap, we can reject the null that the two parameter vectors are equal at level $\alpha$. However, this test is conservative, that is, the actual significance level is usually even smaller. Thus, if we wanted to compare both Pangu-Weather models with each other instead of the ECMWF ensemble, we usually would need to compute the respective skill scores and the corresponding confidence band. However, since the confidence bands do not overlap here, we can directly conclude that the forecast accuracies of both Pangu-Weather models are significantly different from each other.
Figure (ref) in Appendix (ref) shows maps of the estimated CRPSS together with lower and upper bounds of 90% Bonferroni bands, as well as maps of the width of the corresponding bands. With 35,200 grid points the score vector far exceeds the dimensions considered in our simulations, and the theoretical guarantees of Section (ref) do not formally apply. We therefore treat Figure (ref) as an exploratory illustration of spatial uncertainty patterns rather than as inferentially valid bands. For both methods, the largest improvements over the ECMWF ensemble predictions can be observed in coastal areas and over mountainous regions. There is substantial variability across locations. For example, even though the Pangu-Weather + GNP forecasts showed significantly worse performance than the ECMWF model, there are locations where the upper bound of the CRPSS bands is positive, even when aggregated over all forecast horizons. At a forecast horizon of 72h, the upper bound is positive for most land grid points. Similarly, the Pangu-Weather + EasyUQ forecasts are not significantly better than the ECMWF ensemble everywhere, with a substantial fraction of grid cells showing negative values for the lower bounds of the CRPSS bands. For both methods, the confidence bands are notably wider over sea grid points than over land. In summary, the results presented here nicely complement those from buelte2026 by quantifying not only where data-driven models outperform the physics-based benchmark, but also where sampling uncertainty remains too large to draw conclusions.
We develop simultaneous confidence bands for vectors of forecast accuracy measures in order to quantify and communicate sampling uncertainty in forecast comparisons. These bands are versatile, allowing for a variety of multidimensional settings, and may be applied to any type of elicitable forecast. Our approach allows for joint inference, for example, over multiple forecast horizons, variables, methods, and/or locations, while preventing multiple comparison problems. Additionally, our confidence bands improve interpretability since they are not limited to (differences in) expected scores, but can be applied to skill scores (or relative accuracy), which measure the relative improvement in forecast accuracy. The bands are straightforward to implement via a bootstrap algorithm. We consider sup-t bands with exact asymptotic coverage and Bonferroni bands, which are asymptotically conservative, but offer better finite-sample reliability under temporal dependence and are, therefore, our recommended default.
We believe that our confidence bands provide a substantial step forward for quantifying and communicating sampling uncertainty in forecast comparisons. They add simultaneity, interpretability, and visualizability compared to conventional pairwise hypotheses tests of predictive accuracy. They improve over pointwise confidence bands using an iid bootstrap by simultaneity and validity under temporal dependence. Nevertheless, there is room for future research and improvement. First, variance estimation under temporal dependence is a crucial topic in general and specifically in the forecast comparison setting. Plug-in confidence bands based on a HAC variance estimator would be an alternative to our bootstrap algorithm. However, as discussed in the paper, they are not expected to cure the downward bias leading to undercoverage fitzenberger1997. Second, the theoretical guarantees established here treat the dimension $J$ of the skill score vector as fixed while the sample size $N$ grows. Extending the asymptotic theory to allow $J \to \infty$ with $N$ would be of interest for applications involving large spatial grids, such as the weather forecasting case study, where the number of locations exceeds the sample size. Third, data-driven block length selection for the moving block bootstrap remains an open problem in the multivariate setting. While our simulation results support $l = 3\lfloor N^{1/4} \rfloor$ as a practical default, a principled automatic selection procedure would strengthen the implementation. To facilitate the application of our confidence bands in practice, we provide a software implementation in the form of an R package (see the link at the end of the introduction). A Python implementation will follow.
We thank Daniel Gutknecht, Fabian Krüger, Thomas Muschinski and conference participants at the MathSEE Symposium 2023 at Karlsruhe Institute of Technology, the 17th International Conference on Computational and Financial Econometrics 2023 at HTW Berlin, the IWH Workshop on Forecasting in Times of Structural Change and Uncertainty 2024 at Halle Institute for Economic Research, the 44th International Symposium on Forecasting 2024 in Dijon, the 2024 Annual Conference of the German Economic Association at TU Berlin and the DAGStat conference 2025 at HU Berlin as well as seminar participants at Karlsruhe Institute of Technology, Goethe University Frankfurt and Bielefeld University for helpful comments. Sebastian Lerch and Marc Pohle gratefully acknowledge support by the Vector Stiftung through the Young Investigator Group “Artificial Intelligence for Probabilistic Weather Forecasting”. Marc-Oliver Pohle and Sebastian Lerch thank the Klaus Tschira Foundation for infrastructural support at the Heidelberg Institute for Theoretical Studies (HITS).
\addcontentsline{toc}{section}{\refname}