# Sequential model confidence sets

Sebastian Arnold , Georgios Gavrilopoulos11footnotemark: 1 22footnotemark: 2,
Benedikt Schulz, Johanna Ziegel22footnotemark: 2 Both authors contributed
equally to this work.Eidgenössische Technische Hochschule Zürich (ETHZ),
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Institute of Technology (KIT),
<EMAIL_ADDRESS>
###### Abstract

In most prediction and estimation situations, scientists consider various
statistical models for the same problem, and naturally want to select amongst
the best. Hansen et al. (2011) provide a powerful solution to this problem by
the so-called model confidence set, a subset of the original set of available
models that contains the best models with a given level of confidence.
Importantly, model confidence sets respect the underlying selection
uncertainty by being flexible in size. However, they presuppose a fixed sample
size which stands in contrast to the fact that model selection and forecast
evaluation are inherently sequential tasks where we successively collect new
data and where the decision to continue or conclude a study may depend on the
previous outcomes. In this article, we extend model confidence sets
sequentially over time by relying on sequential testing methods. Recently,
e-processes and confidence sequences have been introduced as new, safe methods
for assessing statistical evidence. Sequential model confidence sets allow to
continuously monitor the models’ performances and come with time-uniform,
nonasymptotic coverage guarantees.

_Keywords:_ Model confidence set, forecast evaluation, model selection,
sequential inference, multiple testing.

##  1 Introduction

In science as well as in our daily life, we frequently encounter situations in
which multiple statistical models or forecasts are available for the same
problem and where we have to decide which model(s) or forecast(s) we want to
trust. We might think for example of multiple inflation regression models with
respect to different sets of covariates or of multiple weather institutions
which issue precipitation predictions for the next day. In such situations, we
naturally want to select the best model(s) or forecast(s), where the term
“best” is defined in terms of a user-specified criterion.

The _model confidence set (MCS)_ proposed by Hansen et al. (2011) provides a
promising solution to this problem. It departs from standard practice where
just a single model is selected according to some appropriate loss and nothing
is said about the uncertainty associated with this selection. The MCS takes
this uncertainty into account by reducing the original set of available models
to a smaller set of flexible size that contains the best models with a given
level of confidence. Model confidence sets are of great importance in
applications without an obvious benchmark, may be easily constructed by the
MCS algorithm, and are widely applied in the econometrics community; see,
e.g., Weron (2014), Corbet et al. (2019) or Masini et al. (2023) amongst
others. However, even though the MCS provides an appealing solution to a
highly important problem, it comes with limitations. In particular, it is
assumed that the sequences of observed loss differences are stationary, an
assumption which may be questionable in practise since we may expect that the
models evolve over time, and correct for errors and systematic biases by using
past information. Apart from relying on strong stationarity assumptions, the
MCS requires a sample of some fixed size chosen independently of the data. In
other words, we assume that we observe losses of the models over some
prespecified evaluation period in order to compute the MCS only once at the
end of the period. This procedure highly contrasts our natural urge to assess
the models on a regular basis by successively including past observations to
the available dataset. Consider for example different weather institutions
that predict the accumulated precipitation for the following day on a daily-
basis. Then, firstly, it is not given that the initially better models remain
superior until the end of our study. Secondly, we want to assess the
forecasters and monitor their performance sequentially, say, at the end of
each day, week or month, since the forecasters issue their forecasts
sequentially as well and may behave nonstationary. If we gather evidence that
one institution is worse than the others but the evidence is not yet
statistically significant, then we want to continue to collect evidence for
this hypothesis without loosing the information of the earlier observations.
Crucially, we want to decide at the end of each day, week or month whether to
continue or conclude the study, depending on the previous outcomes.

In summary, model selection and forecast evaluation are inherently sequential
tasks and we should seek for methods that allow to sequentially choose amongst
the best models while respecting the underlying selection uncertainty. In this
article, we address this need and contribute with sequential model confidence
sets. Sequential model confidence sets generalize the core idea of Hansen et
al. (2011) sequentially over time and enjoy time-uniform coverage guarantees
of the theoretically superior objects. Our methods rely on e-processes and
time-uniform confidence sequences, which lay the foundation for safe anytime-
valid inference, a field which met a surge of new contributions over the last
years by various authors; see, e.g., Shafer (2021), Vovk and Wang (2021),
Ramdas et al. (2023) and Grünwald et al. (2019) amongst others.

Sequential testing methods rely on minimal distributional assumptions and have
been used before for sequential forecast evaluation: Henzi and Ziegel (2022)
provide e-processes to sequentially assess two probability forecasts for the
strong hypothesis that one forecaster is better than the other at all time
points. In contrast, Choe and Ramdas (2023) contribute with e-processes for
the weaker hypothesis that one forecaster is uniformly better than the other
on average. Both studies assess the performance of some forecaster at time
$t\in\mathbb{N}$ with respect to the (average) expected loss difference given
all information until the previous time step, which goes back to the test of
equal conditional predictive ability proposed by Giacomini and White (2006).
In contrast to the unconditional approach, most prominently in the seminal
work of Diebold and Mariano (1995), the conditional approach uses available
information to infer which forecaster is more accurate at a specific date. In
this article, we follow the conditional approach and study the strong as well
as the weak hypothesis which we further divide into two subhypotheses. In
contrast to all aforementioned contributions on forecast comparison, our
methods allow to simultaneously draw inference on the performance of $m\geq 2$
different forecasters rather than just being able to compare two different
forecasters.

Although the article mainly focuses on forecast evaluation, many of the
presented results directly extend to model selection as long as the used
information criterion satisfies the imposed assumptions. For this reason, we
introduce the problem in a general framework and illustrate the connections to
forecast comparison and model selection afterwards. In order to construct
sequential model confidence sets, we make use of the e-processes provided by
Choe and Ramdas (2023) under the assumption that the observed loss differences
are conditionally bounded, which does not only hold for bounded losses but
also for prominent choices of scoring rules such as the continuous ranked
probability score (CRPS; Matheson and Winkler, 1976) or the quantile score
(Gneiting, 2011).

The remainder of the paper is structured as follows. Section 2 presents the
problem and introduces sequential model confidence sets. In Section 3, we
review important definitions from sequential testing theory and provide
methods to construct sequential model confidence sets with respect to three
particular notions of superior objects. A simulation study is conducted in
Section 4 before we apply the proposed methods in case studies on Covid-19
related deaths and wind gust predictions in Section 5. The main part of the
paper closes with a discussion in Section 6. A review on forecast evaluation,
technical comments, and additional information for the case studies are
available in Appendices A through G.

##  2 Sequential model confidence sets

We consider a set $\mathcal{M}_{0}=\\{1,\dots,m\\}$ indexing $m\geq 2$
different statistical models or forecasters, simply referred to as _models_ ,
and let $L$ be a loss function which measures the quality of the models. We
let $L$ be negatively oriented, that is, lower scores correspond to better
performance, assume that we collect data in discrete time and denote the loss
of model $i\in\mathcal{M}_{0}$ at time $t\in\mathbb{N}$ by $L_{i,t}$. We
assume that all random quantities are defined on some underlying measurable
space $(\Omega,\mathcal{F})$ equipped with some filtration
$(\mathcal{F}_{t})_{t\in\mathbb{N}}$. We denote the family of all probability
measures on $(\Omega,\mathcal{F})$ by $\mathfrak{B}(\Omega)$ and write
$\mathbb{Q}$ for a generic element in $\mathfrak{B}(\Omega)$.

######  Example 2.1 (Sequential forecasters).

Let $\\{1,\dots,m\\}$ denote different forecasters which all sequentially
issue predictive distributions or point forecasts
$(f_{1,t})_{t\in\mathbb{N}},\dots,(f_{m,t})_{t\in\mathbb{N}}$ for some unknown
quantity $(Y_{t})_{t\in\mathbb{N}}\subseteq\mathbb{R}$. For each $i=1,\dots,m$
and $t\in\mathbb{N}$, the forecast $f_{i,t}$ is predictable and refers to the
observation $Y_{t}$. To assess the forecasters, one should employ _proper
scoring rules_ or _consistent scoring functions_ , see Appendix A. Then, the
time series $(L_{i,t})_{t\in\mathbb{N}}$ emerges by applying a proper scoring
rule (consistent scoring function) $L$ to the forecast-observation pairs
$(f_{i,t},Y_{t})_{t\in\mathbb{N}}$, that is $L_{i,t}=L(f_{i,t},Y_{t})$,
$i=1,\dots,m,t\in\mathbb{N}$.

######  Example 2.2 (Sequential model selection).

Let $\\{1,\dots,m\\}$ consist of nested and non-nested regression models for
some sequentially revealed response $(Y_{t})_{t\in\mathbb{N}}$. At each
$t\in\mathbb{N}$, model $i$ observes the covariate(s)
$\boldsymbol{X}_{i,t}\subseteq\mathbb{R}^{p_{i}}$, $p_{i}\geq 1,$ and
proclaims
$Y_{t}=\boldsymbol{\beta}_{i,t}^{\top}\boldsymbol{X}_{i,t}+\varepsilon_{i,t}$
for some $\boldsymbol{\beta}_{i,t}\subseteq\mathbb{R}^{p_{i}}$ and random
error $\varepsilon_{i,t}$. The quality of the models is typically evaluated
with respect to some likelihood-based _information criterion_ (Ding et al.,
2018). For example, for $\hat{\ell}_{i,t}$ the empirical log-likelihood of
model $i$ given all observations until $t$,
$L_{i,t}=-2\hat{\ell}_{i,t}+2p_{i}$ or $L_{i,t}=-2\hat{\ell}_{i,t}+p_{i}\log
t$ reveal the well-known Akaike information criterion (AIC) or Bayesian
information criterion (BIC), respectively.

Next, we clarify the notion of “best models at time $t\in\mathbb{N}$”. For
$i,j\in\mathcal{M}_{0}$ and $t\in\mathbb{N}$, consider

$d_{ij,t}=L_{i,t}-L_{j,t},\quad\mu_{ij,t}=\mathbb{E}(d_{ij,t}\mid\mathcal{F}_{t-1})$

and their averages

$\hat{\Delta}_{ij,t}=\frac{1}{t}\sum_{s=1}^{t}d_{ij,s},\quad\Delta_{ij,t}=\frac{1}{t}\sum_{s=1}^{t}\mu_{ij,s},$

where we omit the dependence of $\mu_{ij,t}$ and $\Delta_{ij,t}$ on the choice
of the probability measure $\mathbb{Q}$ on $(\Omega,\mathcal{F})$ for the sake
of clarity. Whereas Diebold and Mariano (1995) define model (forecast)
superiority with respect to the unconditional expected loss differences, we
follow Giacomini and White (2006), Lai et al. (2011), Henzi and Ziegel (2022)
and Choe and Ramdas (2023) by defining the superior models in terms of the
(average) conditional expected loss differences. In particular, we consider
the _strongly superior objects_

$\mathcal{M}^{\textrm{s},\star}_{t}=\left\\{i\in\mathcal{M}_{0}\mid\mu_{ij,r}\leq
0\textrm{ for all }r\leq t,\textrm{ for all
}j\in\mathcal{M}_{0}\right\\},\quad t\in\mathbb{N},$ (1)

the _uniformly weakly superior objects_

$\mathcal{M}^{\textrm{uw},\star}_{t}=\left\\{i\in\mathcal{M}_{0}\mid\Delta_{ij,r}\leq
0\textrm{ for all }j\in\mathcal{M}_{0},\textrm{ for all }r\leq
t\right\\},\quad t\in\mathbb{N},$ (2)

and the _weakly superior objects_

$\mathcal{M}^{\textrm{w},\star}_{t}=\left\\{i\in\mathcal{M}_{0}\mid\Delta_{ij,t}\leq
0\textrm{ for all }j\in\mathcal{M}_{0}\right\\},\quad t\in\mathbb{N}.$ (3)

For each $t\in\mathbb{N}$,
$\mathcal{M}^{\textrm{s},\star}_{t},\mathcal{M}^{\textrm{uw},\star}_{t}$ and
$\mathcal{M}^{\textrm{w},\star}_{t}$ are $\mathcal{F}_{t-1}$-measurable random
sets with
$\mathcal{M}^{\textrm{s},\star}_{t}\subseteq\mathcal{M}^{\textrm{uw},\star}_{t}\subseteq\mathcal{M}^{\textrm{w},\star}_{t}$.
Importantly, $\mathcal{M}^{\textrm{w},\star}_{t}\neq\emptyset$ for all
$t\in\mathbb{N}$, whereas $\mathcal{M}^{\textrm{s},\star}_{t}$ and
$\mathcal{M}^{\textrm{uw},\star}_{t}$ may eventually become empty for all
$t\geq t_{0}$ for some $t_{0}\in\mathbb{N}$.

Let $(\mathcal{M}^{\star}_{t})_{t\in\mathbb{N}}\subseteq\mathcal{M}_{0}$ be
the targeted sequence of superior objects. Of particular interest for this
article are the sequences given at (1), (2) and (3), but one could also think
of other possible notions of best models. For a given confidence level
$\alpha\in(0,1)$, we call
$(\widehat{\mathcal{M}}_{t})_{t\in\mathbb{N}}\subseteq\mathcal{M}_{0}$ a
_sequence of model confidence sets_ or _sequential model confidence sets_ ,
for short SMCS(s), for $(\mathcal{M}^{\star}_{t})_{t\in\mathbb{N}}$ at level
$\alpha$ if, for any $\mathbb{Q}\in\mathfrak{B}(\Omega)$,

$\mathbb{Q}\left(\exists t\geq
1:\mathcal{M}^{\star}_{t}\nsubseteq\widehat{\mathcal{M}}_{t}\right)\leq\alpha,$
(4)

or equivalently, if, for any $\mathbb{Q}\in\mathfrak{B}(\Omega)$,

$\mathbb{Q}\left(\forall t\geq
1:\mathcal{M}^{\star}_{t}\subseteq\widehat{\mathcal{M}}_{t}\right)\geq
1-\alpha.$ (5)

We refer to the equivalent conditions (4) and (5) as _time-uniform coverage
guarantees_ of the sequential model confidence sets and highlight the fact
that the time quantifiers are inside of the probability, which is a much
stronger requirement than assuming
$\mathbb{Q}(\mathcal{M}^{\star}_{t}\subseteq\widehat{\mathcal{M}}_{t})\geq
1-\alpha$ for all $t\in\mathbb{N}$ (or assuming that there exists a
$t\in\mathbb{N}$ with
$\mathbb{Q}(\mathcal{M}^{\star}_{t}\nsubseteq\widehat{\mathcal{M}}_{t})\leq\alpha$);
see, e.g., Howard et al. (2021, Lemma 3).

###  2.1 Discussion of the different sets of superior models

Studying different notions of superior models is important, since, depending
on the particular situation, one might interpret the term “best” differently.
In this section, we discuss the hypotheses which correspond to the sequences
of superior models given at (1), (2), and (3), respectively, in more detail,
and give concrete examples for each of them.

If it is reasonable to assume that some models outperform all other models at
all time points (in terms of the conditional expected score differences), then
we should try to find estimators for
$(\mathcal{M}^{\textrm{s},\star}_{t})_{t\in\mathbb{N}}$. We refer to this
assumption as the _strong hypothesis_. The strong hypothesis particularly
applies if we assume independent and identically distributed scores. Then, one
typically studies the unconditional expected score differences
$\mathbb{E}(d_{ij,t})$ independent of $t\in\mathbb{N}$, as it is done by
Diebold and Mariano (1995) or Hansen et al. (2011), who assume stationary
performances of the models. Another important instance where the strong
hypothesis applies is if the models have nested information sets, as it is
e.g. the case in Example 2.2 if the different sets of covariates are nested.
Then, the regression models with many covariates typically estimate the given
data more accurately, and the AIC or BIC take the corresponding complexities
into account by a penalty on the regression dimension. In Section 5.1, we
construct SMCSs for different Covid-19 related deaths where the assumption of
nested information sets seems plausible.

On the other hand, if we assume that some models have a lower conditional
expected score on average rather than at all single time points, then we
should target $(\mathcal{M}^{\textrm{uw},\star}_{t})_{t\in\mathbb{N}}$ under
the _uniformly weak hypothesis_. Choe and Ramdas (2023, Section 4.4.) argue
convincingly that there are many situations where testing for the strong
hypothesis may be misleading. Inspired by their arguments, Example 2.4 below
considers a forecaster which is ideal on most of the days and only slightly
worse than the other forecasters on Sundays. Any powerful method to monitor
the strongly superior models over time will exclude this forecaster from the
SMCS since she is indeed not strongly superior. Nevertheless, the forecaster
is superior under the uniformly weak hypothesis as she is generally better
than all other forecasters and only slightly worse than some of them on
Sundays.

Finally, we say that the _weak hypothesis_ applies if we expect that also the
average performances might change over time. The weak hypothesis is of
particular interest if we assume that the models might evolve over time and
improve their outputs individually by using past observations. In Section 5.2,
we study post-processing methods for wind gust predictions. The results show
that in this application it is indeed beneficial to monitor the weakly
superior objects as some models which are excluded from the SMCS at earlier
time points are included again at later stages, either due to systematic
changes in the underlying meteorological model, or due to the different
adaptive behaviour of the methods. If we had conducted the study under any of
the two stronger hypothesis, we could not have observed that some methods
become again competitive towards the end of the test period, since, under
these hypotheses, a model can not be included again in the SMCS once it is
excluded.

Since the weak hypothesis is the most general, one might wonder, why we do not
only consider this particular hypothesis as it is implied by the other two
hypotheses anyhow? The answer to this question is given in terms of power:
Whereas we impose more assumptions on the data-generating process for the
former hypotheses, our methods for the latter are generally less powerful if
we apply them to construct SMCSs for the strongly (uniformly weakly) superior
models. However, if we construct SMCSs for the strongly (uniformly weakly)
superior models even though the corresponding hypotheses do not apply, then we
may expect misleading results as argued before.

To conclude, we introduce the simulation settings from Section 4 to give
concrete examples for each of the discussed hypotheses.

######  Example 2.3 (Simulation 1).

For $n=1000$, we sample $(y_{t})_{t=1}^{n}$, where
$y_{t}\sim\mathcal{N}(y_{t-1},1),$ for $t=1,\dots,n$ and $y_{0}=0$. We
consider $m=49$ different forecasters $\\{1,\dots,m\\}$ which all sequentially
issue predictive distributions for $(y_{t})_{t=1}^{n}$ given by

$f_{i,t}=\mathcal{N}\left(y_{t-1}+\varepsilon_{i},1+\delta_{i}\right),\quad
t=1,\dots,n,$

for $(\varepsilon_{i},\delta_{i})\in\\{-0.6,-0.4,-0.2,0,0.2,0.4,0.6\\}^{2}$.
For these forecasters with different biases and dispersion errors, there is
exactly one ideal forecaster ${i_{0}}$ with
$\varepsilon_{i_{0}}=\delta_{i_{0}}=0$. We assess the different forecasters
with respect to the _continuous ranked probability score_ (CRPS; Matheson and
Winkler (1976)), see Appendix A. Since the CRPS is a proper scoring rule, it
follows that $\mathcal{M}^{\textrm{s},\star}_{t}=\\{i_{0}\\}$ for all
$t\in\mathbb{N}$.

######  Example 2.4 (Simulation 2).

We consider the same data-generating mechanism and the same forecasters as in
Example 2.3 with the only exception that the forecaster ${i_{0}}$ now issues
the forecasts
$f_{i_{0},t}=\mathcal{N}\left(y_{t-1}+\varepsilon_{i_{0},t},1+\delta_{i_{0},t}\right)$,
$t=1,\dots,n$, for

$\varepsilon_{i_{0},t}=\delta_{i_{0},t}=\begin{cases}0.3&,t\in 7\mathbb{N}\\\
0&,\textrm{else}.\end{cases}$

That is, forecaster $f_{i_{0}}$ is still ideal on, say, weekdays and
Saturdays, however, on Sundays, some forecasters have a smaller bias and
dispersion error, and hence $\mathcal{M}^{\textrm{s},\star}_{t}=\emptyset$ for
all $t\geq 7$. However, ${i_{0}}$ is uniformly weakly superior, that is
$\mathcal{M}^{\textrm{uw},\star}_{t}=\\{{i_{0}}\\}$ for all $t\in\mathbb{N}$.

######  Example 2.5 (Simulation 3).

We sample i.i.d. standard normally distributed observations
$(y_{t})_{t=1}^{n}$, for $n=800$, and compare $m=3$ different forecasters
$\\{1,2,3\\}$ which issue median predictions
$m_{i,t}=y_{t}+\varepsilon_{i,t}$, $i=1,2,3$, for

$\varepsilon_{1,t}=\beta,\quad\varepsilon_{2,t}=\gamma^{t},\quad\varepsilon_{3,t}=\delta
t,\quad t=1,\dots,n,$ (6)

for $\beta,\delta>0$ and $0<\gamma<1$. We call forecaster $1$ _constantly
biased_ , forecaster $2$ _improving_ and forecaster $3$ _disimproving_. We
assess them by the scoring function $L(m,y)=0.5(\Phi(m)-\Phi(y))$, for
$m,y\in\mathbb{R}$ and $\Phi$ the cdf of the standard normal distribution,
which is a consistent scoring function for the median, see Appendix A. If we
choose the parameters as given in Section 4, we have
$\mathcal{M}^{\textrm{w},\star}_{t}=\\{3\\}$ for $t\leq 153$,
$\mathcal{M}^{\textrm{w},\star}_{t}=\\{1\\}$ for $153<t<550$, and
$\mathcal{M}^{\textrm{w},\star}_{t}=\\{2\\}$ for $t\geq 550$. That is, whereas
the disproving forecaster is initially superior, the constantly biased
forecaster catches up after some time until the improving forecaster becomes
the best in the end.

##  3 Construction of sequential model confidence sets

In this section, we provide methods to construct SMCSs for the superior models
given at (1), (2) and (3). Our constructions build on e-processes, confidence
sequences and sequential multiple testing methods as introduced in the
following.

###  3.1 Sequential testing methods and multiple testing

Let $(\Omega,\mathcal{F})$ be a measurable space and $\mathfrak{B}(\Omega)$ be
the family of all probability measures on $(\Omega,\mathcal{F})$. A
_statistical (null) hypothesis_ $\mathcal{H}\subseteq\mathfrak{B}(\Omega)$ is
a set of probability measures that are potential candidates for the true
probability measure $\mathbb{P}$ governing the data generating process. We let
$(\mathcal{F}_{t})_{t\in\mathbb{N}}\subseteq\mathcal{F}$ be a filtration and
interpret $\mathcal{F}_{t}$ as the available information at time
$t\in\mathbb{N}$. We call a process $(M_{t})_{t\in\mathbb{N}}$ a _test
(super-)martingale_ for some hypothesis
$\mathcal{H}\subseteq\mathfrak{B}(\Omega)$ if it is a nonnegative
(super-)martingale with $M_{0}\leq 1$ for each $\mathbb{Q}\in\mathcal{H}$.
Test martingales allow for a financial interpretation in terms of betting as
the accumulated wealth of a gambler after having bet a number of times against
the null hypothesis (Shafer, 2021). If the null hypothesis is true, then one
cannot expect to gain money over time, and thus, large values of the test
martingale give evidence against the null hypothesis. More rigorously, any
test supermartingale $(M_{t})_{t\in\mathbb{N}}$ for
$\mathcal{H}\subseteq\mathfrak{B}(\Omega)$ satisfies

$\sup_{\mathbb{Q}\in\mathcal{H}}\mathbb{Q}\left(\exists
t\in\mathbb{N}:M_{t}\geq\frac{1}{\alpha}\right)\leq\alpha,\quad\textrm{for all
}\alpha\in(0,1).$ (7)

The inequality at (7), referred to as _Ville’s inequality_ (Ville, 1939),
follows from the optional stopping theorem for martingales and forms a time-
uniform extension of Markov’s inequality. It lays the foundation for safe
anytime-valid inference by enabling the construction of time-uniform
confidence sequences, sequential tests, or anytime-valid p-values amongst
others (Ramdas et al., 2023).

For some problems, it is theoretically impossible to construct nontrivial test
martingales (Ramdas et al., 2022). However, often one still can draw
sequential inference by using so-called e-processes: An _e-process_ for some
hypothesis $\mathcal{H}\subseteq\mathfrak{B}(\Omega)$ is a nonnegative adapted
stochastic process $(E_{t})_{t\in\mathbb{N}}$ with

$\mathbb{E}_{\mathbb{Q}}(E_{\tau})\leq 1$

for all $\mathbb{Q}\in\mathcal{H}$ and all (possibly infinite) stopping times
$\tau$. In other words, for any random stopping time $\tau$, $E_{\tau}$ is an
_e-variable_ (Vovk and Wang, 2021). Equivalently, an adapted nonnegative
process $(E_{t})_{t\in\mathbb{N}}$ is an e-process for $\mathcal{H}$ if it is
upper bounded by a test supermartingale for each $\mathbb{Q}\in\mathcal{H}$
and hence, any test supermartingale is also an e-process and Ville’s
inequality continues to hold for e-processes (Ramdas et al., 2020).

Test martingales and e-processes allow for sequential testing. A _sequential
testing procedure_ or _sequential test_ for some hypothesis
$\mathcal{H}\subseteq\mathfrak{B}(\Omega)$ is an adapted stochastic process
$(\Psi_{t})_{t\in\mathbb{N}}\subseteq\\{0,1\\}$, where $\Psi_{t}=1$ indicates
that the hypothesis $\mathcal{H}$ is rejected at $t\in\mathbb{N}$. A
sequential test $(\Psi_{t})_{t\in\mathbb{N}}$ for
$\mathcal{H}\subseteq\mathfrak{B}(\Omega)$ is said to _uniformly control the
type-I error at level $\alpha\in(0,1)$_ if

$\sup_{\mathbb{Q}\in\mathcal{H}}\mathbb{Q}(\exists
t\in\mathbb{N}:\Psi_{t}=1)\leq\alpha.$

By Ville’s inequality at (7), any e-process $(E_{t})_{t\in\mathbb{N}}$
directly leads to a sequential testing procedure $(\Psi_{t})_{t\in\mathbb{N}}$
with uniform type-I error control at level $\alpha\in(0,1)$ by defining

$\Psi_{t}=\mathbbm{1}\\{E_{t}\geq 1/\alpha\\},\quad t\in\mathbb{N}.$ (8)

If we test more than one hypothesis simultaneously, then we have to correct
for multiple testing. In the spirit of Vovk and Wang (2021), we call $m\geq 2$
e-processes $(E_{1,t})_{t\in\mathbb{N}},\dots,(E_{m,t})_{t\in\mathbb{N}}$ for
the hypotheses $\mathcal{H}_{1},\dots,\mathcal{H}_{m}$ _family-wise valid
(FWV)_ if there exists a _dominating family of e-processes_ , that is, if
there exists a family
$\\{(E^{\mathbb{Q}}_{t})_{t\in\mathbb{N}}\mid\mathbb{Q}\in\mathfrak{B}(\Omega)\\}$
such that $(E^{\mathbb{Q}}_{t})_{t\in\mathbb{N}}$ is an e-process for
$\\{\mathbb{Q}\\}$ for all $\mathbb{Q}\in\mathfrak{B}(\Omega)$, and for any
$i=1,\dots,m$ with $\mathbb{Q}\in\mathcal{H}_{i}$,
$(E^{\mathbb{Q}}_{t})_{t\in\mathbb{N}}$ upper bounds
$(E_{i,t})_{t\in\mathbb{N}}$. In Appendix C, we show that for family-wise
valid e-processes, the corresponding sequential tests given at (8) time-
uniformly control the family-wise error rate (FWER).

To derive SMCS for the sequence of weakly superior objects, we will make use
of confidence sequences with time-uniform coverage guarantees (Howard et al.,
2021; Ramdas et al., 2023). For $\alpha\in(0,1)$, a _(time-uniform)
$(1-\alpha)$-confidence sequence_ for a stochastic process
$(\theta_{t})_{t\in\mathbb{N}}\subseteq\Theta$ in some space $\Theta$ is a
sequence $(C_{t})_{t\in\mathbb{N}}\subseteq 2^{\Theta}$ such that

$\mathbb{Q}\left(\forall t\in\mathbb{N}:\theta_{t}\in C_{t}\right)\geq
1-\alpha,\quad\textrm{for all }\mathbb{Q}\in\mathfrak{B}(\Omega).$

###  3.2 SMCSs for the strongly superior models

In this section, we construct SMCSs for the sequence of strongly superior
models $(\mathcal{M}^{\textrm{s},\star}_{t})_{t\in\mathbb{N}}$ given at (1) by
a sequential testing procedure. Throughout this section, we assume that for
any $i,j\in\mathcal{M}_{0}$, the sequence of loss differences
$(d_{ij,t})_{t\in\mathbb{N}}$ is conditionally bounded with respect to some
predictable sequence $(c_{ij,t})_{t\in\mathbb{N}}\subseteq\mathbb{R}_{+}$,
that is $|d_{ij,t}|\leq c_{ij,t}/2$ for all $t\in\mathbb{N}$. For
$i,j\in\mathcal{M}$, we consider the strong hypotheses

$\mathcal{H}^{\textrm{s}}_{ij}=\\{\mathbb{Q}\in\mathfrak{B}(\Omega)\mid\mu_{ij,t}\leq
0,\forall t\in\mathbb{N}\\}$

and the strong intersection hypotheses

$\displaystyle\mathcal{H}^{\textrm{s}}_{i\cdot}=\bigcap_{j\in\mathcal{M}_{0}}\mathcal{H}^{\textrm{s}}_{ij}=\\{\mathbb{Q}\in\mathfrak{B}(\Omega)\mid\mu_{ij,t}\leq
0,\forall t\in\mathbb{N},\forall j\in\mathcal{M}_{0}\\}.$

Then, for any predictable $(\lambda_{ij,t})_{t\in\mathbb{N}}$ with
$0\leq\lambda_{ij,t}\leq 1/c_{ij,t}$, we have that

$M_{ij,t}=\prod_{r=1}^{t}(1+\lambda_{ij,r}(d_{ij,r}-\mu_{ij,r})),\quad
t\in\mathbb{N}$

is a nonnegative martingale. Importantly, for any
$\mathbb{Q}\in\mathcal{H}^{\textrm{s}}_{ij}$, $M_{ij,t}$ upper bounds

$E_{ij,t}=\prod_{r=1}^{t}(1+\lambda_{ij,r}d_{ij,r}),\quad t\in\mathbb{N},$

Hence, $(E_{ij,t})_{t\in\mathbb{N}}$ is an e-process for
$\mathcal{H}^{\textrm{s}}_{ij}$, and thus also for
$\mathcal{H}^{\textrm{s}}_{i\cdot}$. For $i=1,\dots,m$, we let

$E_{i\cdot,t}=\frac{1}{m-1}\sum_{j\neq i}E_{ij,t},\quad t\in\mathbb{N},$

which is also an e-process for $\mathcal{H}^{\textrm{s}}_{i\cdot}$ as the
arithmetic mean of $m-1$ e-processes. Finally, we adjust the e-processes

$(E_{1\cdot t})_{t\in\mathbb{N}},\dots,(E_{m\cdot,t})_{t\in\mathbb{N}},$ (9)

symmetrically by the closure principle in the spirit of Vovk and Wang (2021),
as explained in detail in Appendix C, to obtain family-wise valid e-processes
$(E_{1\cdot,t}^{\star})_{t\in\mathbb{N}},\dots,(E_{m\cdot,t}^{\star})_{t\in\mathbb{N}}$
for the hypotheses
$\mathcal{H}^{\textrm{s}}_{1\cdot},\dots,\mathcal{H}^{\textrm{s}}_{m\cdot}$.
We use the arithmetic mean as the $e$-merging function as it essentially
dominates all symmetric e-merging functions (Vovk and Wang, 2021, Proposition
3.1.). For $\alpha\in(0,1)$, we apply the sequential tests given at (8) and
let

$\widehat{\mathcal{M}}_{t}=\left\\{i\in\mathcal{M}_{0}\mid
E_{i\cdot,t}^{\star}\leq 1/\alpha\right\\}.$ (10)

That is, at each $t\in\mathbb{N}$, we exclude model $i\in\mathcal{M}_{0}$ from
the SMCS if we may reject the hypothesis $\mathcal{H}^{\textrm{s}}_{i\cdot}$
according to the sequential testing procedure, and include it otherwise.

######  Lemma 3.1.

For any $\alpha\in(0,1)$, the sequence
$(\widehat{\mathcal{M}}_{t})_{t\in\mathbb{N}}$ defined at (10) is a sequence
of model confidence sets at level $\alpha$ for the sequence of strongly
superior objects $(\mathcal{M}^{\textrm{s},\star}_{t})_{t\in\mathbb{N}}$, and
so is its running intersection $\widetilde{\mathcal{M}}_{t}=\bigcap_{r\leq
t}\widehat{\mathcal{M}}_{r},t\in\mathbb{N}$.

###### Proof.

For any $\mathbb{Q}\in\mathfrak{B}(\Omega)$, we have

$\displaystyle\mathbb{Q}(\exists
t\in\mathbb{N}:\mathcal{M}^{\textrm{s},\star}_{t}\not\subseteq\widehat{\mathcal{M}}_{t})=\mathbb{Q}\left(\exists
t\in\mathbb{N}:E_{i\cdot,t}^{\star}>1/\alpha\textrm{ for some
}i\in\mathcal{M}^{\textrm{s},\star}_{t}\right).$

Moreover, $i\in\mathcal{M}^{\textrm{s},\star}_{t}$ if and only if
$\mathbb{Q}\in\mathcal{H}^{\textrm{s}}_{i\cdot,t}=\\{\mathbb{Q}\in\mathfrak{B}(\Omega)\mid\mu_{ij,r}\leq
0,\forall r\leq t,\forall j\in\mathcal{M}_{0}\\}$. Since
$(E_{1\cdot,t}^{\star})_{t\in\mathbb{N}},\dots,(E_{m\cdot,t}^{\star})_{t\in\mathbb{N}}$
are FWV, there exists a dominating family of e-processes
$\\{(E_{t}^{\mathbb{Q}})_{t\in\mathbb{N}}\mid\mathbb{Q}\in\mathfrak{B}(\Omega)\\}$
such that $E_{r}^{\mathbb{Q}}\geq E_{i\cdot,r}^{\star}$ for all $r\leq t$ and
all $i\in\mathcal{M}^{\textrm{s},\star}_{t}$. So, we may conclude

$\displaystyle\mathbb{Q}(\exists
t\in\mathbb{N}:\mathcal{M}^{\textrm{s},\star}_{t}\not\subseteq\widehat{\mathcal{M}}_{t})$
$\displaystyle=\mathbb{Q}\left(\exists
t\in\mathbb{N}:E_{i\cdot,t}^{\star}>1/\alpha\textrm{ for some
}i\in\mathcal{M}^{\textrm{s},\star}_{t}\right)$
$\displaystyle\leq\mathbb{Q}\left(\exists
t\in\mathbb{N}:E_{t}^{\mathbb{Q}}>1/\alpha\right)\leq\alpha.$

It follows directly, that the running intersection
$\widetilde{\mathcal{M}}_{t}=\bigcap_{r\leq
t}\widehat{\mathcal{M}}_{r},t\in\mathbb{N}$, is a SMCS as well, since
$i\notin\mathcal{M}^{\star}_{t_{0}}$ implies $i\notin\mathcal{M}^{\star}_{t}$
for all $t\geq t_{0}.$ ∎

######  Remark 1.

Alternatively, one could adjust the $m\cdot(m-1)$ e-processes
$(E_{ij,t})_{t\in\mathbb{N}}$ directly to obtain FWV e-processes
$(E^{\star}_{ij,t})_{t\in\mathbb{N}}$ for the hypotheses
$\mathcal{H}^{\textrm{s}}_{ij}$ and exclude model $i$ from the SMCS at
$t\in\mathbb{N}$ if $E^{\star}_{ij,t}>1/\alpha$ for some $j\neq i$. However,
the corresponding tests are less powerful, since for any $i\in\mathcal{M}_{0}$
and $j\neq i$,

$\displaystyle E^{\star}_{ij,t}$ $\displaystyle=\min\left\\{\frac{1}{\lvert
I\rvert}\sum_{(k,l)\in
I}E_{kl,t}\;\Big{|}\;I\subseteq\\{1,\dots,m\\}^{2},(i,j)\in I\right\\}$
$\displaystyle\leq\min\left\\{\frac{1}{\lvert I\rvert(m-1)}\sum_{k\in
I}\sum_{l=1}^{m}E_{kl,t}\;\Big{|}\;I\subseteq\\{1,\dots,m\\},i\in
I\right\\}=E_{i\cdot}^{\star}.$

######  Remark 2.

Apart from the FWER, the _false discovery rate (FDR)_ is arguably the most
commonly used criterion in multiple testing (Benjamini and Hochberg, 1995;
Benjamini and Yekutieli, 2001). Recently, Wang and Ramdas (2022) studied FDR
control for e-values, which turn out to be particularly suitable for this
problem due to their behaviour under taking expectations. In Appendix D, we
demonstrate how one could alternatively construct SMCSs with a time-uniform
control of the FDR, which is generally less restrictive than FWER.

###  3.3 SMCSs for the uniformly weakly superior objects

We construct SMCSs for the sequence of uniformly weakly superior objects
$(\mathcal{M}^{\textrm{uw},\star}_{t})_{t\in\mathbb{N}}$ defined at (2)
analogously as in the previous section. For simplicity, we assume first that
for any $i,j\in\mathcal{M}$, there exists a constant $c_{ij}>0$ with $\lvert
d_{ij,t}\rvert\leq c_{ij}/2$ for all $t\in\mathbb{N}$. In Section 3.5, we
explain how we can weaken this assumption to conditionally bounded loss
differences.

For $i,j\in\mathcal{M}$, consider the weak hypotheses

$\mathcal{H}^{\textrm{uw}}_{ij}=\\{\mathbb{Q}\in\mathfrak{B}(\Omega)\mid\Delta_{ij,t}\leq
0,\forall t\in\mathbb{N}\\}$

and the weak intersection hypotheses

$\displaystyle\mathcal{H}^{\textrm{uw}}_{i\cdot}=\bigcap_{j\in\mathcal{M}_{0}}\mathcal{H}^{\textrm{uw}}_{ij}=\\{\mathbb{Q}\in\mathfrak{B}(\Omega)\mid\Delta_{ij,t}\leq
0,\forall t\in\mathbb{N},\forall j\in\mathcal{M}_{0}\\}.$

Choe and Ramdas (2023) show that, for any $0\leq\lambda_{ij}<1/c_{ij}$,

$M_{ij,t}=\exp\left\\{\lambda_{ij}t(\hat{\Delta}_{ij,t}-\Delta_{ij,t})-\psi_{E,c_{ij}}(\lambda_{ij})V_{ij,t}\right\\},\quad
t\in\mathbb{N}$

is a test supermartingale under any $\mathbb{Q}\in\mathfrak{B}(\Omega)$, where

$\psi_{E,c}(\lambda)=(-\log(1-c\lambda)-c\lambda)/c^{2}$ (11)

for $c>0$ and $0\leq\lambda<1/c$ and

$V_{ij,t}=\sum_{r=1}^{t}(d_{ij,r}-\gamma_{ij,r})^{2}$ (12)

for some predictable sequence $(\gamma_{ij,t})_{t\in\mathbb{N}}$ with
$\lvert\gamma_{ij,t}\lvert\leq c_{ij}/2$. Importantly, for any
$\mathbb{Q}\in\mathcal{H}^{\textrm{uw}}_{ij}$, $M_{ij,t}$ upper bounds

$E_{ij,t}=\exp\left\\{\lambda_{ij}t\hat{\Delta}_{ij,t}-\psi_{E,c_{ij}}(\lambda_{ij})V_{ij,t}\right\\}$
(13)

for all $t\in\mathbb{N}$. That is, $(E_{ij,t})_{t\in\mathbb{N}}$ is an
e-process for $\mathcal{H}^{\textrm{uw}}_{ij}$ and hence also for
$\mathcal{H}^{\textrm{uw}}_{i\cdot}$. Analogously as before, we adjust the
averages $E_{i\cdot,t}=1/(m-1)\sum_{j\neq i}E_{ij,t}$, $i=1,\dots,m$,
symmetrically by the closure principle to construct SMCSs for
$(\mathcal{M}^{\textrm{uw},\star}_{t})_{t\in\mathbb{N}}$ by the FWV
e-processes
$(E_{1\cdot,t}^{\star})_{t\in\mathbb{N}},\dots,(E_{m\cdot,t}^{\star})_{t\in\mathbb{N}}$
as in (10), compare Appendix C.

###  3.4 SMCSs for the weakly superior objects

Our SMCS construction for the weakly superior objects defined at (3) builds on
time-uniform confidence regions. For any $i,j\in\mathcal{M}_{0}$, we again
assume that there exists $c_{ij}>0$ with $\lvert d_{ij,t}\rvert\leq c_{ij}/2$
for all $t\in\mathbb{N}$ and discuss a possible weakening of this assumption
in Section 3.5. For $t\in\mathbb{N}$ and $0\leq\lambda_{ij}\leq 1/c_{ij}$, let

$M_{ij,t}(x)=\exp\left\\{\lambda_{ij}t\hat{\Delta}_{ij}-\lambda_{ij}tx-\psi_{E,c_{ij}}(\lambda_{ij})V_{ij,t}\right\\},\quad
x\in\mathbb{R},t\in\mathbb{N},$

where $\psi_{E,c_{ij}}$ and $V_{ij,t}$ are given in (11) and (12),
respectively. The following lemma shows the one-to-one correspondence between
confidence sequences and sequential tests, and is found e.g. in Howard et al.
(2021, Lemma 1) or Waudby-Smith and Ramdas (2024, Theorem 1).

######  Lemma 3.2.

For any $i,j\in\mathcal{M}_{0}$ and $\alpha\in(0,1)$, the sequence

$C_{t}=\left\\{x\in[-c_{ij}/2,c_{ij}/2]\mid M_{ij,t}(x)\leq
1/\alpha\right\\},\quad t\in\mathbb{N},$

is a $(1-\alpha)$-confidence sequence for $(\Delta_{ij,t})_{t\in\mathbb{N}}$.

The claim follows directly by Ville’s inequality observing that
$(M_{ij,t}(\Delta_{ij,t}))_{t\in\mathbb{N}}$ is a nonnegative supermartingale
for any $\mathbb{Q}\in\mathfrak{B}(\Omega)$. Note that $M_{ij,t}(x)\leq
1/\alpha$ if and only if
$x\geq\hat{\Delta}_{ij,t}-(\psi_{E,c_{ij}}(\lambda_{ij})V_{ij,t}+\log(1/\alpha))/(\lambda_{ij}t)$.
That is, Lemma 3.2 actually yields time-uniform lower confidence bounds for
$(\Delta_{ij,t})_{t\in\mathbb{N}}$. Next, we show how to combine these lower
confidence bounds simultaneously for all models $i\in\mathcal{M}_{0}$, and how
we can use them to construct SMCSs for
$(\mathcal{M}^{\textrm{w},\star}_{t})_{t\in\mathbb{N}}$.

We let $\mathbb{R}_{0}^{m\times m}$ be the family of all $m\times m$-matrices
with diagonal entries equal zero and let

$\boldsymbol{c}=\begin{pmatrix}0&c_{12}&\cdots&c_{1m}\\\
c_{21}&0&\cdots&c_{2m}\\\ \vdots&\vdots&\ddots&\vdots\\\
c_{m1}&c_{m2}&\dots&0\\\ \end{pmatrix}\in\mathbb{R}_{0}^{m\times
m},\quad\boldsymbol{\Delta}_{t}=\begin{pmatrix}0&\Delta_{12,t}&\cdots&\Delta_{1m,t}\\\
\Delta_{21,t}&0&\cdots&\Delta_{2m,t}\\\ \vdots&\vdots&\ddots&\vdots\\\
\Delta_{m1,t}&\Delta_{m2,t}&\cdots&0\\\
\end{pmatrix}\in\mathbb{R}_{0}^{m\times m},\quad t\in\mathbb{N}.$

For $\boldsymbol{A},\boldsymbol{B},\boldsymbol{C}\in\mathbb{R}_{0}^{m\times
m}$ and $C\in\mathbb{R}$, we understand $\boldsymbol{A}\leq\boldsymbol{B}$ and
$C\cdot\boldsymbol{A}$ componentwise, and write
$\boldsymbol{B}\leq\boldsymbol{A}\leq\boldsymbol{C}$ equivalently as
$\boldsymbol{A}\in[\boldsymbol{B},\boldsymbol{C}]$. For
$\boldsymbol{X}\in\mathbb{R}_{0}^{m\times m}$, we let

$M_{t}(\boldsymbol{X})=\frac{1}{m(m-1)}\sum_{\begin{subarray}{c}i,j=1,\dots,m\\\
i\neq j\end{subarray}}M_{ij,t}(X_{ij}),\quad t\in\mathbb{N}.$

######  Lemma 3.3.

For any $\alpha\in(0,1)$, the sequence

$C_{t,1-\alpha}=\left\\{\boldsymbol{X}\in[-\boldsymbol{c}/2,\boldsymbol{c}/2]\mid
M_{t}(\boldsymbol{X})\leq 1/\alpha\right\\}$

is a time-uniform $(1-\alpha)$-confidence sequence for the sequence
$(\boldsymbol{\Delta}_{t})_{t\in\mathbb{N}}$.

###### Proof.

For any $\mathbb{Q}\in\mathfrak{B}(\Omega)$, we have

$\displaystyle\mathbb{Q}(\exists t\in\mathbb{N}:\boldsymbol{\Delta}_{t}\notin
C_{t,1-\alpha})$ $\displaystyle=\mathbb{Q}(\exists
t\in\mathbb{N}:M_{t}(\boldsymbol{\Delta}_{t})>1/\alpha)\leq\alpha,$

by Ville’s inequality since
$(M_{t}(\boldsymbol{\Delta}_{t}))_{t\in\mathbb{N}}$ is a test supermartingale
under any $\mathbb{Q}\in\mathfrak{B}(\Omega)$. ∎

Next, we define

$\widehat{\mathcal{M}}_{t}=\left\\{i\in\mathcal{M}_{0}\mid
C_{t,1-\alpha}\cap\left(\mathbb{R}_{0}^{m\times
m}\right)_{ij,-}\neq\emptyset\textrm{ for all }j\neq i\right\\},\quad
t\in\mathbb{N},$ (14)

where $\left(\mathbb{R}_{0}^{m\times
m}\right)_{ij,-}=\\{\boldsymbol{X}\in\mathbb{R}_{0}^{m\times m}\mid X_{ij}\leq
0\\}$. That is, for each model $i\in\mathcal{M}_{0}$, we check whether
$\Delta_{ij,t}$ may potentially be negative for all $j\neq i$ according to the
confidence region $C_{t,1-\alpha}$. If there exists $j\neq i$ such that we may
reject that $\Delta_{ij,t}\leq 0$ with high probability, then we exclude $i$
from the SMCS.

######  Lemma 3.4.

For any $\alpha\in(0,1)$, the sequence
$(\widehat{\mathcal{M}}_{t})_{t\in\mathbb{N}}$ defined at (14) is a a sequence
of model confidence sets at level $\alpha$ for the sequence of weakly superior
objects $(\mathcal{M}^{\textrm{w},\star}_{t})_{t\in\mathbb{N}}$.

###### Proof.

By Lemma 3.3, for any $\mathbb{Q}\in\mathfrak{B}(\Omega)$, we have

$\displaystyle\mathbb{Q}(\exists
t\in\mathbb{N}:\mathcal{M}^{\textrm{w},\star}_{t}\not\subseteq\widehat{\mathcal{M}}_{t})$
$\displaystyle=$ $\displaystyle\mathbb{Q}(\exists
t\in\mathbb{N}:i\in\mathcal{M}^{\textrm{w},\star}_{t}\textrm{ and
}i\not\in\widehat{\mathcal{M}}_{t}\textrm{ for some }i\in\mathcal{M}_{0})$
$\displaystyle=$ $\displaystyle\mathbb{Q}(\exists t\in\mathbb{N},\exists
i\in\mathcal{M}_{0}:\Delta_{ij,t}\leq 0\textrm{ for all }j\neq i\textrm{ and
}C_{t}\cap\left(\mathbb{R}^{m\times m}\right)_{i\ell,-}=\emptyset\textrm{ for
some }\ell\neq i)$ $\displaystyle\leq$ $\displaystyle\mathbb{Q}(\exists
t\in\mathbb{N}:\boldsymbol{\Delta}_{t}\notin C_{t})\leq\alpha.$

∎

The following Lemma, which is proved in Appendix E, characterizes the shape of
the confidence regions $(C_{t,1-\alpha})_{t\in\mathbb{N}}$ and may be used to
compute the sequence $(\widehat{\mathcal{M}}_{t})_{t\in\mathbb{N}}$ more
efficiently.

######  Lemma 3.5.

For any $\alpha\in(0,1)$ and $t\in\mathbb{N}$, $C_{t,1-\alpha}$ is a convex
upper set. That is, for $\boldsymbol{A},\boldsymbol{B}\in C_{t,1-\alpha}$, it
follows that $\lambda\boldsymbol{A}+(1-\lambda)\boldsymbol{B}\in
C_{t,1-\alpha}$ for all $\lambda\in[0,1]$, and $\boldsymbol{A}\in
C_{t,1-\alpha}$ and $\boldsymbol{A}\leq\boldsymbol{C}$ implies that
$\boldsymbol{C}\in C_{t,1-\alpha}$ for any
$\boldsymbol{C}\leq\boldsymbol{c}/2$.

By Lemma 3.5, we may express $\widehat{\mathcal{M}}_{t}$ equivalently as

$\widehat{\mathcal{M}}_{t}=\left\\{i\in\mathcal{M}_{0}\mid
C_{t,1-\alpha}\cap[\boldsymbol{0},\boldsymbol{c}/2]_{ij,0}\neq\emptyset\textrm{
for all }j\neq i\right\\},\quad t\in\mathbb{N},$

where $\boldsymbol{0}=(0)_{i,j=1}^{m}\in\mathbb{R}_{0}^{m\times m}$ and
$[\boldsymbol{0},\boldsymbol{c}/2]_{ij,0}=\\{\boldsymbol{X}\in[\boldsymbol{0},\boldsymbol{c}/2]\mid
X_{ij}=0\\}$. That is, we exclude model $i\in\mathcal{M}_{0}$ from the SMCS at
$t\in\mathbb{N}$, if there exists $j\neq i$ such that the confidence region
$C_{t,1-\alpha}$ has an empty intersection with the hyperplane
$\\{\boldsymbol{X}\in\mathbb{R}_{0}^{m\times m}\mid X_{ij}=0\\}$.

###  3.5 Predictable bounds and betting schemes

The methods provided in Sections 3.3 and 3.4 require that the score
differences are uniformly bounded over time, which is clearly more restrictive
than the assumption of conditionally bounded score differences imposed in
Section 3.2. However, one can always transform conditionally bounded score
differences into uniformly bounded differences thereby modifying the
underlying loss and the corresponding superior objects.

Assume that $\lvert d_{ij,t}\rvert\leq c_{ij,t}/2$ for all $t\in\mathbb{N}$
for some predictable $(c_{ij,t})_{t\in\mathbb{N}}$. Then, the transformed loss
differences $\tilde{d}_{ij,t}=d_{ij,t}/c_{ij,t},t\in\mathbb{N}$ are uniformly
bounded with constant $c=1$ that is $\lvert\tilde{d}_{ij,t}\rvert\leq 1/2$ for
all $t\in\mathbb{N}$. However, if we use the transformed loss differences for
the methods of the previous sections, then we target (possibly) different
sequences of superior objects

$\tilde{\mathcal{M}}^{\textrm{uw},\star}_{t}=\left\\{i\in\mathcal{M}_{0}\mid\tilde{\Delta}_{ij,r}\leq
0,\forall r\leq t,\forall
j\in\mathcal{M}_{0}\right\\},\quad\tilde{\mathcal{M}}^{\textrm{w},\star}_{t}=\left\\{i\in\mathcal{M}_{0}\mid\tilde{\Delta}_{ij,t}\leq
0,\forall j\in\mathcal{M}_{0}\right\\},$

for

$\tilde{\Delta}_{ij,t}=\frac{1}{t}\sum_{r=1}^{t}\tilde{\mu}_{ij,r}=\frac{1}{t}\sum_{r=1}^{t}\mathbb{E}(\tilde{d}_{ij,r}\mid\mathcal{F}_{r-1})=\frac{1}{t}\sum_{r=1}^{t}\frac{\mu_{ij,r}}{c_{ij,r}}.$

In Section 5, we use the given transformation to convert conditionally bounded
CRPS differences of wind gust forecasts into uniformly bounded loss
differences. Importantly, for any proper scoring rule (consistent scoring
function) ${\rm S}$, the scaled function $\tilde{{\rm S}}={\rm S}/c$ is proper
(consistent) as well for any $c>0$, an observation which justifies the
transformation from a theoretical perspective. Additionally, with the given
scaling, all observations have the same maximal impact on the ranking, an
effect which seems appealing for many applications.

######  Remark 3.

We assume conditionally bounded loss differences in order to work with the
empirical Bernstein e-processes at (13) provided by Choe and Ramdas (2023).
Clearly, for many important loss functions this assumption is not fulfilled;
see, e.g., Example 2.2 where the log-likelihood differences are typically
unbounded. However, if the loss differences emerge from (or may be bounded in
the tails by) a parametric family of distributions, as it is for example the
case in Example 2.2, there are other well-studied e-processes available; see,
e.g., Howard et al. (2020), Pérez-Ortiz et al. (2023) or Hao et al. (2024)
amongst others. In such situations, our methods directly extend with a
suitable choice of e-processes at (13).

We conclude with a comment on the choice of the parameter $\lambda$. For the
uniformly weak hypothesis and the weak hypothesis, we assume a universal bound
$c>0$ on the loss differences, and $\lambda$ may be any fixed value in the
interval $(0,c^{-1})$. If there is no good reason to do otherwise, we suggest
to choose $\lambda=(2c)^{-1}$ by default. In betting language, this
corresponds to bet half of the accumulated evidence at each time step (Shafer,
2021). For the strong hypothesis, we should choose $\lambda$ predictably and
in a smart way in order to increase power. Clearly, there are many possible
betting schemes. In Section 5.1, we propose a particular betting scheme which
makes use of all previous observations and the assumption that the relative
performance of two models does not change too quickly over time. Another
option would be the _method of mixtures_ (Robbins, 1970), where we integrate
over all possible values of $\lambda$ over a particular probability
distribution instead of choosing one specific value. Mixtures of e-processes
are again e-processes and may be expressed in closed form for some
distributions (Howard et al., 2021, Appendix A.3). The method of mixtures is
one of the most widely-studied techniques for constructing uniform bounds
which are shown to be unimprovable under certain conditions (Robbins and
Siegmund, 1970, Howard et al., 2021). However, often it is analytically opaque
and computationally tedious, and does not lead to closed-form confidence
sequences. In accordance with Waudby-Smith and Ramdas (2024), for the strong
hypothesis, our experience shows that choosing $\lambda$ predictably in a
smart ways is generally more promising than the methods of mixtures. Also for
the weaker hypotheses, we could not improve the power of our methods
significantly by using mixtures instead of choosing a fixed $\lambda$ value.
For this reason, and to make the computations more tractable, we suggest to
work with the given e-processes directly rather than with the more complicated
mixture e-processes also for both weak hypotheses.

##  4 Simulations

###  4.1 Simulation 1

We consider the setup given in Example 2.3. That is, we consider $n=1000$
realizations from a sequence $Y_{t}\sim\mathcal{N}(Y_{t-1},1)$ with $Y_{0}=0$,
and assess $m=49$ forecasters which issue sequentially probabilistic
predictions with respect to the CRPS. We write $Y_{t}=Y_{t-1}+Z_{t}$, for an
i.i.d. standard normally distributed sequence $(Z_{t})_{t\in\mathbb{N}}$, and
use the formulas provided by Grimit et al. (2006) to derive

$\displaystyle L_{i,t}$ $\displaystyle=\textrm{CRPS}(f_{i,t},Y_{t})$
$\displaystyle=\sqrt{1+\delta_{i}}\left(\frac{Z_{t}-\varepsilon_{i}}{\sqrt{1+\delta_{i}}}\left(2\Phi\left(\frac{Z_{t}-\varepsilon_{i}}{\sqrt{1+\delta_{i}}}\right)-1\right)+2\varphi\left(\frac{Z_{t}-\varepsilon_{i}}{\sqrt{1+\delta_{i}}}\right)-\frac{1}{\sqrt{\pi}}\right),$
(15)

where $\Phi$ and $\varphi$ denote the cdf and the density function of the
standard normal distribution, respectively. Let
$(\mathcal{F}_{t})_{t\in\mathbb{N}}$ be the canonical filtration, that is
$\mathcal{F}_{t}$ is generated by all observations $Y_{r}$ and forecasts
$f_{i,r}$, $i=1,\dots,m$, $r\leq t$. For any $t\in\mathbb{N}$, we have
$\mathbb{E}(d_{ij,t}\mid\mathcal{F}_{t-1})=\mathbb{E}(d_{ij,t})=\mathbb{E}(d_{ij,1})$
by independence. Since the CRPS is a proper scoring rule, we conclude

$\mathcal{M}^{\textrm{s},\star}_{t}=\mathcal{M}^{\star}=\left\\{i\in\mathcal{M}_{0}\mid\mathbb{E}(d_{ij,1})\leq
0,\forall j\in\mathcal{M}_{0}\right\\}=\\{i_{0}\\}$

for the forecaster $i_{0}$ with $\varepsilon_{i_{0}}=\delta_{i_{0}}=0$.

We construct SMCSs for $(\mathcal{M}^{\textrm{s},\star}_{t})_{t\leq n}$ as
proposed in Section 3.2 with confidence level $\alpha=0.1$. We let
$\lambda_{ij,t}=(2c_{ij})^{-1}$ for the universal bound $c_{ij}>0$ on the
score differences $(d_{ij,t})$ given in Appendix B. Figure 1 shows the average
size of the SMCS over time for $N=100$ realizations of the simulation.
Importantly, we have a 100% coverage rate of the best model $i_{0}$, that is,
in all runs of the simulation, the SMCS never wrongly excludes $i_{0}$.

Figure 1: The average number of models in the SMCS in simulations 1 and 2 over
$N=100$ repetitions. In simulation 1, the size of the SMCS tends to the actual
number of superior objects as the sample size increases, whereas there are
still approximately $10$ models remaining in the SMCS for the uniformly weakly
superior objects after $1000$ observations. For both simulations, the SMCS
never wrongly excludes the best model $i_{0}$.

###  4.2 Simulation 2

We consider the simulation setting as given in Example 2.4. By (4.1),

$\displaystyle\mathbb{E}(L_{i,t})$
$\displaystyle=\mathbb{E}\big{(}Z_{t}\Phi_{\varepsilon_{i},1+\delta_{i}}(Z_{t})\big{)}-\varepsilon_{i}\mathbb{E}\big{(}\Phi_{\varepsilon_{i},1+\delta_{i}}(Z_{t})\big{)}+\varepsilon_{i}+2\sqrt{1+\delta_{i}}\mathbb{E}\big{(}\varphi_{\varepsilon_{i},1+\delta_{i}}(Z_{t})\big{)}-\sqrt{\frac{1+\delta_{i}}{\pi}}$
$\displaystyle=\frac{1+2\sqrt{1+\delta_{i}}}{\sqrt{2\pi(2+\delta_{i})}}\exp\left\\{-\frac{\alpha_{i}^{2}}{2}\right\\}-\varepsilon_{i}\Phi\left(-\alpha_{i}\right)+\varepsilon_{i}-\sqrt{\frac{1+\delta_{i}}{\pi}},$
(16)

where $\alpha_{i}=\varepsilon_{i}/\sqrt{2+\delta_{i}}$ and
$\Phi_{\mu,\sigma^{2}}$ ($\varphi_{\mu,\sigma^{2}}$) denotes the cdf (density)
of $\mathcal{N}(\mu,\sigma^{2})$ for $\mu\in\mathbb{R}$ and $\sigma>0$. For
$i,j\in\mathcal{M}_{0}$ and $t\in\mathbb{N}$, let
$\mathbb{E}(d_{ij,t}\mid\mathcal{F}_{t-1})=\mathbb{E}(d_{ij,t})=\mu_{ij,t}$.
Since the CRPS is a proper scoring rule, we have $\mu_{i_{0}j,t}\leq 0$ for
all $t\notin 7\mathbb{N}$, for all $j\in\mathcal{M}_{0}$. By (16), for
$j\in\mathcal{M}_{0}$ with
$(\varepsilon_{j},\delta_{j})\in\\{-0.2,0.2\\}^{2}$, we have
$\mu_{i_{0}j,t}>0$ for all $t\in 7\mathbb{N}$. Thus,
$\mathcal{M}^{\textrm{s},\star}_{t}=\emptyset$ for $t\geq 7$. However,
equation (16) allows to show numerically that $\Delta_{i_{0}j,t}\leq 0$ for
all $j\in\mathcal{M}_{0}$ and $t\in\mathbb{N}$, and thus
$\mathcal{M}^{\textrm{uw},\star}_{t}=\\{i_{0}\\}$ for all $t\in\mathbb{N}$.

We compute SMCSs for $(\mathcal{M}^{\textrm{uw},\star}_{t})_{t\leq n}$ as
suggested in Section 3.3 with confidence level $\alpha=0.1$. We transform the
conditionally bounded CRPS differences (see Appendix B) as proposed in Section
3.5 into uniformly bounded score differences $\lvert\tilde{d}_{ij,t}\rvert\leq
C/2$ for $C=2$, and let $\lambda=(2C)^{-1}=1/4$. Figure 1 displays the average
size of the SMCS for $N=100$ repetitions of the simulation. Again, we have a
100% coverage rate of the superior model $i_{0}$. Notably, we have less power
compared to the SMCS for the strongly superior objects and still remain with
approximately $10$ surviving models after $1000$ observations.

###  4.3 Simulation 3

We consider the simulation setting as given in Example 2.5. Let $i=1,2,3$ be
forecasters issuing median predictions $m_{i,t}=Y_{t}+\varepsilon_{i,t}$,
$t\in\mathbb{N}$, for i.i.d. standard normally distributed
$(Y_{t})_{t\in\mathbb{N}}$, where the biases $\varepsilon_{i,t}$ are given in
(6). For $t\in\mathbb{N}$, we have by independence

$2\mu_{ij,t}=\mathbb{E}(\Phi(m_{i,t}))-\mathbb{E}(\Phi(m_{j,t}))=\Phi\left(\frac{\varepsilon_{i,t}}{\sqrt{2}}\right)-\Phi\left(\frac{\varepsilon_{j,t}}{\sqrt{2}}\right),\quad
i,j=1,2,3,$

and thus

$\mathcal{M}^{\textrm{w},\star}_{t}=\left\\{i=1,2,3\mid\sum_{r=1}^{t}\Phi\left(\frac{\varepsilon_{i,t}}{\sqrt{2}}\right)\leq\sum_{r=1}^{t}\Phi\left(\frac{\varepsilon_{j,t}}{\sqrt{2}}\right),j=1,2,3\right\\}.$
(17)

We let $\beta=0.6,\gamma=0.998$ and $\delta=0.008$. Then, one can numerically
show by (17) that $\mathcal{M}^{\textrm{w},\star}_{t}=\\{3\\}$ for $t\leq
153$, $\mathcal{M}^{\textrm{w},\star}_{t}=\\{1\\}$ for $153<t<550$, and
$\mathcal{M}^{\textrm{w},\star}_{t}=\\{2\\}$ for $t\geq 550$.

We sample $n=800$ observations and construct SMCSs for
$(\mathcal{M}^{\textrm{w},\star}_{t})_{t\leq n}$ at level $\alpha=0.1$ as
suggested in Section 3.4. We use the fact that $(d_{ij,t})_{t\in\mathbb{N}}$
is uniformly bounded with $c=1$ and let $\lambda_{ij}=1.1^{-1}$, $i,j=1,2,3$.
Figure 2 shows the accumulated observed losses and the resulting SMCS for one
realization of the simulation (left), and the average size of the SMCS over
time for $N=100$ runs (right). Again, the SMCS includes the theoretical
superior objects given at (17) at all time points.

Figure 2: For $i=1,2,3$, the left panel shows the realized accumulated losses
$\sum_{r=1}^{t}L_{i,r}$, $t\leq n$, for one simulation. The black vertical
lines indicate $t=154$ and $t=550$. The resulting SMCS is given in the upper
part with the respective colors: The SMCS correctly excludes the improving
forecaster (blue) at the beginning and the disimproving forecaster (green)
relatively quickly after the first change point until it includes the
improving forecaster again and finally excludes the constantly biased
forecaster (red). The right panel shows the average size of the SMCS over
$N=100$ realizations, where we again have a $100\%$ coverage rate of the
superior models.

##  5 Case Studies

###  5.1 Covid-19 case study

After the outbreak of the Covid-19 pandemic, the Reich lab at the University
of Massachusetts collaborated with the United States centers for disease
control and prevention (CDC) to create the Covid-19 forecast hub, a repository
containing point and probabilistic forecasts for incident cases,
hospitalizations, and Covid-19 related deaths. The repository was founded in
March 2020, and attracted the submission of forecasts by more than $60$
different teams. The collection of forecasts is still ongoing, although most
of the models had stopped submitting forecasts by January 2024, due to the
subsidence of the pandemic.

####  5.1.1 The data

The dataset, which we use in this application, is discussed in detail in
Cramer et al. (2022a) and Cramer et al. (2022b), and publicly available at the
Covid-19 forecast hub GitHub repository and on the Zoltar forecast archive. It
contains forecasts with corresponding observations for key epidemiology
indices for 55 different locations in the US as well as on the aggregated US
national level. The forecasts are issued by a total of 69 different forecast
models and reported at $23$ different quantile levels with forecast horizons
ranging from $1$ day to $20$ weeks. Importantly, many models are not directly
comparable since not all of them have issued their forecasts on the same
weekdays or throughout the whole period.

In this case study, we consider 1-week-ahead forecasts for Covid related
deaths on the national level, issued weekly in the period from 06/06/2020 to
04/03/2024, and focus on the comparison of the following $m=6$ different
forecasting models: Firstly, we consider the benchmark model _baseline_ ,
which naively issues the most recent outcome as the median prediction for the
following week and forms a predictive distribution around this median
prediction based on the past weekly incidences. Secondly, we consider two
epidemic models _PSI-draft_ and _MOBS-gleam_ issuing their forecasts based on
epidemiological assumptions, and the neural network based model _GT-deep_ ,
which was the first purely data-driven model to be included in the Covid-19
forecast hub. Finally, we consider the summary models _ensemble_ and _CDC-
ensemble_ , where the former combines the most recent predictions of all other
models into one predictive distribution, and the latter is supposed to improve
the former by only considering the 10 current best models, measured by the
average weighted interval score over the 12 most recent weeks. For a more
detailed documentation of the models, we refer to Appendix F.

Regarding the nature of the models, we perform our comparison under the strong
hypothesis, since we expect the ensemble models to have consistent superior
performance over the individual models, which in turn are expected to be
better than the naive baseline model issuing its forecasts only based on the
previous outcomes and on no additional information. Furthermore, we can argue
that in such a highly sensible prediction task, we should seek to detect the
models that issue accurate predictions consistently over time and not just on
average.

Notably, the outcomes as well as the forecasts in the given dataset have a
very unstable behaviour and exhibit large variation in scale, see e.g. the
upper left panel of Figure 3, which suggests to apply a log transformation on
all quantities of interest. Moreover, the forecast differences are typically
not very distinctive for the median, whereas they become more apparent when we
look at the tails, compare to the lower panels in Figure 3.

Figure 3: Time progression of predicted and actual Covid related mortality in
linear (upper left) and log scale (upper right). For the forecasts at the
tails, the differences between the different models are more pronounced than
for the median, see the lower panels for $\tau=0.15,0.975$.

####  5.1.2 Methods and implementation details

For $\tau\in(0,1)$, we assess the different quantile forecasts at level $\tau$
with respect to the _generalized piecewise-linear (GPL) quantile scoring
function_ (Gneiting, 2011),

$S_{\tau}(x,y)=\left(\mathbf{1}\\{x\geq
y\\}-\tau\right)\left(\log(x)-\log(y)\right),\quad x,y\in\mathbb{R}.$

The log transformation may thus easily be justified from a statistical point
of view by simply setting $g=\log$ in the formula of the GPL score, see
Appendix A for details. In Appendix B, it is shown that the resulting quantile
score differences are conditionally bounded. More precisely, for
$\tau$-quantile predictions
$(x_{i,t})_{t\in\mathbb{N}},(x_{j,t})_{t\in\mathbb{N}}\subseteq\mathbb{R}$
issued by the models $i,j=1,\dots,m$, we have $\lvert d_{ij,t}\rvert=\lvert
S_{\tau}(x_{i,t},y_{t})-S_{\tau}(x_{j,t},y_{t})\rvert\leq c_{ij,t}/2$ for

$c_{ij,t}=2\max\left\\{\tau,1-\tau\right\\}\lvert\log(x_{1,t})-\log(x_{2,t})\rvert,\quad
t\in\mathbb{N}.$ (18)

We apply our methods of Section 3.2 to construct SMCSs for the strongly
superior objects, where we use the e-processes
$E_{ij,t}=\prod_{r=1}^{t}(1+\lambda_{ij,r}d_{ij,r})$ for
$(\lambda_{ij,t})_{t\in\mathbb{N}}$ predictable with
$0<\lambda_{ij,t}<c_{ij,t}^{-1}$, $t\in\mathbb{N}$. For $i,j=1,\dots,m$ and
$t\in\mathbb{N}$, we suggest to set

$\lambda_{ij,t}=\frac{1}{K_{ij,t}c_{ij,t}+\varepsilon}$

for some small $\varepsilon>0$, included as a safeguard for the scenario where
the two forecaster issue the same prediction and thus $c_{ij,t}=0$, and for

$K_{ij,t}=\frac{2-\left|\tau-\frac{1}{2}\right|}{1+\left|\tau-\frac{1}{2}\right|}\cdot\frac{3\pi/2+\arctan(-d_{ij,t-1})}{\pi}\geq
1.$

The coefficient $K_{ij,t}$ depends on the sign and magnitude of the most
recent score difference $d_{ij,t-1}$ as well as of the centrality of $\tau$:
If $d_{ij,t-1}>0$, then we have evidence against the null hypothesis that
model $i$ is strictly superior than model $j$, and the quantity $K_{ij,t}$
becomes small. That is, the parameter $\lambda_{ij,t}$ tends to be larger, and
we bet more against the null hypothesis. On the other hand, $K_{ij,t}$ becomes
larger for $d_{ij,t-1}<0$, which leads to smaller values of the tuning
parameters $\lambda_{ij,t}$. Finally, the primary multiplicative factor, which
takes values in $[1,2]$, is used to reduce the influence of the term
$\max\left\\{\tau,1-\tau\right\\}$ on the bound in (18). It converges to $1$
as $\tau\downarrow 0$ or $\tau\uparrow 1$, and to $2$ as $\tau\to 0.5$. This
means that we bet more at the tails where the forecasts are typically more
distinctive and less if $\tau$ is close to $0.5$ where the forecasts are
typically less distinguishable.

Choosing a smart betting strategy has an effect on the power of the methods as
discussed in Section 3.5. In Appendix F, we compare the suggested choice of
the predictable $\lambda_{ij,t}$’s with the naive approach where $K_{ij,t}=2$
for all $t\in\mathbb{N}$.

####  5.1.3 Results

Figure 4: The resulting SMCSs at four selected quantile levels with confidence
level $\alpha=0.1$. Even though the SMCS behave differently at the different
quantile levels, we may conclude that the ensemble models performed best as
they are never excluded from the SMCSs. Second best is GT-deep, which only
shows some deficiencies for $\tau=0.3$, and MOBS-gleam and the baseline, which
seem to be approximately equally accurate. PSI-draft performs poorest and is
excluded rather quickly at each quantile level.

Figure 4 shows the resulting SMCSs for four selected quantile levels with a
confidence level of $\alpha=0.1$. As a first observation, we see that the
power of the SMCS varies with the choice of the quantile level $\tau$: Whereas
the SMCS detects four deficient models for $\tau=0.3$, it only has excluded
two models by the end of the observation period for $\tau=0.7$. For
$\tau=0.1$, the SMCS needs considerably more observations to reject the three
models, which are excluded much faster at the median. As argued before, the
lower power at the tails of the distribution is probably due to the factor
$\max\left\\{\tau,1-\tau\right\\}$ in the bound (18). In Appendix F, we show
that the effect of lower power at the tails even increases if we do not
incorporate the quantile level in our betting scheme and that our proposed
betting scheme is capable to mitigate this difficulty. Overall, the results
align with our expectation that the two ensemble models perform best as they
are included in the SMCS for the entire time period at all quantile levels.
The superior performance of the two ensemble models was expected since they
combine the forecasts of all other models and thus have access to the largest
information set. Interestingly, the SMCS does not detect any significant
difference between ensemble and CDC-ensemble, where the latter was supposed to
improve the former by combining just the currently best, rather than all,
forecasting models. The naive baseline performs surprisingly well and even
outperforms the epidemic model PSI-draft, which is excluded earlier than
baseline for all given values of $\tau$. We may conclude that PSI-draft is not
competitive with the other models and should therefore not be considered as an
accurate forecasting model for the pandemic. Finally, the second epidemic
model MOBS-gleam performs better than PSI-draft but not as good as the neural
network based GT-deep according to the SMCSs.

To conclude, we again highlight the key feature that SMCSs allow us to collect
statistical evidence sequentially over time. Hence, if we were back in the
pandemic, by using SMCSs, we could monitor the models performances from the
very beginning on and would not have to wait until we have collected enough
observations to perform a statistical test for predictive ability. Quickly, we
would reject PSI-draft as an accurate model and not consider it in our
analysis and policy-making anymore. Also MOBS-gleam would probably have been
rejected as a competitive model by approximately the end of 2021, and so we
would have been left with the choice between less and less models in the
course of the pandemic. Clearly, in this extreme situation, it has already
been common practice to assess the forecasting models on a regular basis.
However, this practice may lack the theoretical justification, and our
methodology provides statistically safe methods to do so.

###  5.2 Postprocessing of wind gust predictions

Currently, weather services use numerical weather prediction (NWP) models for
predicting the future state of the atmosphere. These NWP models try to
quantify the forecast uncertainty using ensemble predictions, where each
ensemble member represents a different future scenario. However, the resulting
ensemble forecasts fail to adequately quantify the forecast uncertainty, as
they are subject to systematic biases and dispersion errors. Hence, the
forecasts ought to be statistically postprocessed to generate accurate and
reliable predictions (Gneiting and Raftery, 2005; Vannitsem et al., 2018).

Schulz and Lerch (2022) present a comprehensive comparison of multiple
approaches for statistical postprocessing of ensemble forecasts for wind
gusts. Here, we will build on their case study and compare the different
forecasting models for probabilistic wind gust prediction on an extended time
period by using SMCSs under the weak hypothesis.

####  5.2.1 Data

In this application, we use the data set of Primo et al. (2024), which is an
extension of the original data in Schulz and Lerch (2022). The forecast data
is based on an operational NWP model from the German weather service
(Deutscher Wetterdienst, DWD) in the period from 08/12/2010 to 30/06/2023 with
corresponding observations at 174 geographically diverse weather stations in
Germany. The forecasts are initialized at 00 UTC and have a forecast horizon
of up to 21 hours. Here, we consider only those forecasts with a lead time of
18h referring to 18 UTC and the $K=166$ stations with at least 2,600
observations.

Next to the ensemble forecasts (EPS), we include the eight postprocessing
methods applied by Schulz and Lerch (2022) in our comparison. Statistical
postprocessing methods are typically distributional regression models that use
ensemble predictions from the NWP model as input data. The methods in our
comparison can be divided into three groups: Basic techniques that only make
use of the underlying wind gust ensemble predictions (EMOS, MBM, IDR),
established machine learning methods for postprocessing that incorporate
additional covariates (EMOS-GB, QRF), and neural network-based approaches
(DRN, BQN, HEN). While the first two groups fit a separate model for each
station, the neural networks estimate one locally adaptive model for all
stations. Further, the models differ in the resulting forecast distributions,
which range from parametric (EMOS, EMOS-GB, DRN) to semi- and non-parametric
types (all others). The postprocessing models are trained on the period from
2010 to 2015 using the same configuration as in Schulz and Lerch (2022) and
evaluated on the remaining period from 2016 to 2023.

To improve the forecast quality, NWP models are continuously developed further
and frequently updated. Within the period from 2010 to 2023, the underlying
NWP model has undergone several updates, of which three change the systematic
errors of the ensemble forecasts drastically (on 22/03/2017, 16/05/2018 and
10/02/2021). These model updates present a challenge for statistical
postprocessing methods, as the corrections learned on previous NWP model
versions may not lead to the same improvements when applied to the current
model version. In essence, the training data becomes less representative of
the test data after each update. The methods considered here have been trained
on data until the end of 2015, which is before the first major model update in
2017. Still, they are applied to ensemble forecasts that have been generated
by another NWP model version, which has undergone up to three major updates.
Thus, we expect the behavior both of the ensemble predictions and the
postprocessed forecasts to change systematically over time. For further
details on the data and postprocessing methods, we refer to Schulz and Lerch
(2022) and Primo et al. (2024).

####  5.2.2 Methods and implementation details

We assess the $m=9$ different forecasting models by the CRPS and construct
SMCSs for the weak hypothesis as proposed in Section 3.4. By predictability of
the forecasts, the CRPS differences are conditionally bounded. More precisely,
for any predictive cdfs $F_{1}$ and $F_{2}$, we have

$\lvert\textrm{CRPS}(F_{1},y)-\textrm{CRPS}(F_{2},y)\rvert\leq\max\\{\lvert\textrm{CRPS}(F_{1},z)-\textrm{CRPS}(F_{2},z)\rvert\mid
z\in\mathcal{X}\\},\quad y\in\mathbb{R},$ (19)

where $\mathcal{X}$ denotes the crossing points of $F_{1}$ and $F_{2}$
including $\\{\pm\infty\\}$, see Appendix B.

For each station, we thus compute the crossing points of all issued forecasts
at all days in order to obtain predictable bounds of the CRPS differences by
(19). In a second step, we use the transformation of Section 3.5 to obtain
uniformly bounded score differences. That is, $\lvert d_{ij,t}^{k}\rvert\leq
C/2$ for the resulting CRPS differences at station $k\leq K$, for $C=2$,
$i,j\leq m$ and $t\in\mathbb{N}$. We let $\lambda=(2C)^{-1}=1/4$, compare
Section 3.5, and choose $\alpha=0.1$ as our confidence level. In order to draw
inference over all stations simultaneously, we consider the average scaled
CRPS differences ${d}_{ij,t}=\sum_{k=1}^{K}{d}^{k}_{ij,t}/K$, for $i,j\leq 9$
and $t\in\mathbb{N}$.

####  5.2.3 Results

First, we take a look at the SMCS with respect to the loss differences
averaged over all stations, from which we draw conclusions regarding the
performance of postprocessing methods that align with the results in Schulz
and Lerch (2022). The left panel of Figure 5 displays the resulting SMCS and
shows that the two neural network-based approaches DRN and BQN, which
performed best in the former study, are included in the SMCS for the entire
time period. The other models fall out of the SMCS over time, where the order
coincides with that found in the former study, e.g., the underlying ensemble
is omitted first. In this case, where we average the performances over all
stations, the SMCS behaves as for the stronger hypotheses, that is, once we
eliminate a model, we do not include it anymore.

Second, we consider the SMCS station-wise. The right panel of Figure 5 gives a
summary, as it shows the number of stations for which a forecasting model is
included in the corresponding SMCS. We find that the EPS is excluded from most
SMCSs already within the first year, that the basic techniques EMOS, MBM and
IDR are also omitted relatively fast, and that the machine learning approaches
EMOS-GB and QRF fall out of SMCSs successively over time, until early 2021
when the last NWP model update occurred. Interestingly, the number of stations
that include the HEN network approach increases after the last NWP model
update, a behavior only detectable using the weak hypothesis. Finally, the BQN
and DRN approach are included at almost all stations over the entire period,
only with DRN being excluded over the last year at some stations.

Figure 5: SMCS averaged over the loss differences over all stations (left) and
the number of stations where the method is included in the SMCS for each
method dependent on the time (right). The vertical black lines indicate the
three major NWP model updates.

Taking a closer look at the individual stations, we observe different kinds of
local behavior that is not visible when averaging over the loss differences.
According to the findings, we broadly divide the stations into three
categories: Regular (56 stations), dynamic (55) and fluctuating stations (55).
Figure 6 shows one exemplary station for each of the three categories as well
as the mean size of the station-wise SMCSs over time. The SMCS of a _regular_
station behaves (almost) equally as for the stronger hypotheses, which is also
seen in the plot of the mean SMCS size over the time period. The _dynamic_
stations show a different pattern. As for the regular stations, some of the
inferior models are excluded after a certain time period. However, in contrast
to the regular stations, some methods that have been eliminated before are now
included again coming closer to the end of the time period. This behavior
becomes apparent when looking at the mean size of the SMCSs, which increases
after the last update of the underlying NWP model that systematically changes
the predictive performance of the methods. These cases are interesting, since
any test which would be performed just at the end of the evaluation period (or
at some given time point in the middle) would probably not detect the
deficiencies of these models. Consider for example the SMCS at station 10044,
Leuchtturm Kiel, displayed at the topright of Figure 6, where EMOS, IDR and
HEN have been excluded from the SMCS for long periods and have become
competitive again later on. Even though these three models perform not as
strongly as the benchmarks DRN and BQN, the SMCS indicates that they are still
better than the remaining models for this station. At last, we have the
_fluctuating_ stations, where at least one model is repeatedly excluded and
included from the SMCS. It shows a kind of ”borderline”-behaviour of the
method and indicates that it is exactly at the border of being worse than the
other models. As for the other two groups, this behavior is reflected in the
mean size of the SMCS, which is more oscillating for the fluctuating stations
than the other two categories. In all cases, the sequential nature of the SMCS
yields a deeper insight into the performances of the methods than tests for
predictive ability which are issued at single time points only.

Figure 6: SMCSs for one exemplary station in each of the three categories
(topleft: normal; topright: dynamic; bottomleft: fluctuating) and the mean
size of the SMCS averaged over the stations of each category dependent on the
time (bottomright). Figure 8 in Appendix G shows the spatial distribution of
all station categories over Germany.

##  6 Discussion

We introduced sequential model confidence sets and provided techniques to
construct them with respect to three important notions of forecast
superiority. Sequential model confidence sets address a highly relevant
problem and allow to continuously monitor the performances of some given
statistical models. They incorporate the uncertainty of the undergoing
selection by being flexible in size, rely on minimal distributional
assumptions, and achieve safe anytime-valid guarantees.

Concerning future work, there are various possible directions to investigate
and extend the methods developed so far: Firstly, our construction is limited
to (conditionally) bounded score differences. Important loss functions such as
the quadratic score or the logarithmic score do not satisfy this requirement,
and it would be useful to study e-processes and SMCSs for these more general
loss functions, compare to Remark 3. Secondly, the simulation study shows that
our methods are conservative. Since we have to adjust our processes for
multiple testing, we may lose a considerable amount of evidence. More work is
needed to reduce the gap in power between the raw and the adjusted
e-processes, possibly by adaptive weighting schemes of the e-processes.
Alternatively, one could study sequential model confidence sets with time-
uniform FDR control as discussed in Appendix D. However, our impression is
that there is not much gain in power by replacing the FWER with the less
restrictive FDR. Finally, it would be interesting to examine how our results
extend to important information criteria, and to apply sequential model
confidence sets for sequential model selection.

## Data and replication material

Code in R (R Core Team, 2024) for the case studies and replication material
for the simulations are available at https://github.com/GGavrilos/SMCS.

## Acknowledgments

We thank Lutz Dümbgen for valuable discussions and inputs. Sebastian Arnold
and Johanna Ziegel gratefully acknowledge financial support from the Swiss
National Science Foundation. Computations for the simulations and parts of the
wind gust case study have been performed on UBELIX (https://ubelix.unibe.ch/),
the HPC cluster of the University of Bern. Benedikt Schulz gratefully
acknowledges funding within the project C5 “Dynamical feature-based ensemble
postprocessing of wind gusts within European winter storms” of the
Transregional Collaborative Research Center SFB/TRR 165 “Waves to Weather”
funded by the German Research Foundation (DFG). We thank Cristina Primo,
Reinhold Hess and Sebastian Trepte for providing the forecast and observation
data, and Robert Redl for assistance in data handling.

## References

  * Balcan et al. (2009) D. Balcan, V. Colizza, B. Gonçalves, H. Hu, J. J. Ramasco, and A. Vespignani.  Multiscale mobility networks and the spatial spreading of infectious diseases.  _Proc. Natl. Acad. Sci. U.S.A._ , 106(51):21484–21489, 2009.
  * Balcan et al. (2010) D. Balcan, B. Gonçalves, H. Hu, J. J. Ramasco, V. Colizza, and A. Vespignani.  Modeling the spatial spread of infectious diseases: The global epidemic and mobility computational model.  _J. Comput. Sci._ , 1(3):132–145, 2010.
  * Benjamini and Hochberg (1995) Y. Benjamini and Y. Hochberg.  Controlling the false discovery rate: a practical and powerful approach to multiple testing.  _J. R. Stat. Soc. Ser. B: Stat. Methodol._ , 57:289–300, 1995.
  * Benjamini and Yekutieli (2001) Y. Benjamini and D. Yekutieli.  The control of the false discovery rate in multiple testing under dependency.  _Ann. Stat._ , 29:1165–1188, 2001.
  * Choe and Ramdas (2023) Y. J. Choe and A. Ramdas.  Comparing sequential forecasters.  _Oper. Res._ , 2023.
  * Corbet et al. (2019) S. Corbet, B. Lucey, A. Urquhart, and L. Yarovaya.  Cryptocurrencies as a financial asset: A systematic analysis.  _Int. Rev. Financ. Anal._ , 62:182–199, 2019.
  * Cramer et al. (2022a) E. Y. Cramer, Y. Huang, Y. Wang, E. L. Ray, M. Cornell, J. Bracher, A. Brennen, A. J. C. Rivadeneira, A. Gerding, K. House, et al.  The united states covid-19 forecast hub dataset.  _Sci. Data_ , 9(1):462, 2022a.
  * Cramer et al. (2022b) E. Y. Cramer, E. L. Ray, V. K. Lopez, J. Bracher, A. Brennen, A. J. Castro Rivadeneira, A. Gerding, T. Gneiting, K. H. House, Y. Huang, et al.  Evaluation of individual and ensemble probabilistic forecasts of covid-19 mortality in the united states.  _Proc. Natl. Acad. Sci. U.S.A._ , 119(15):e2113561119, 2022b.
  * Diebold and Mariano (1995) F. Diebold and R. Mariano.  Comparing predictive accuracy.  _J. Bus. Econ. Stat._ , 13(3):253–63, 1995.
  * Ding et al. (2018) J. Ding, V. Tarokh, and Y. Yang.  Model selection techniques: An overview.  _IEEE Signal Process. Mag._ , 35(6):16–34, 2018\.
  * Ehm et al. (2016) W. Ehm, T. Gneiting, A. Jordan, and F. Krüger.  Of quantiles and expectiles: consistent scoring functions, choquet representations and forecast rankings.  _J. R. Stat. Soc. Ser. B: Stat. Methodol._ , 78(3):505–562, 2016.
  * Giacomini and White (2006) R. Giacomini and H. White.  Tests of conditional predictive ability.  _Econometrica_ , 74(6):1545–1578, 2006.
  * Gneiting (2011) T. Gneiting.  Making and evaluating point forecasts.  _J. Am. Stat. Assoc._ , 106:746–762, 2011.
  * Gneiting and Raftery (2005) T. Gneiting and A. E. Raftery.  Weather forecasting with ensemble methods.  _Science_ , 310:248–249, 2005.
  * Gneiting and Raftery (2007) T. Gneiting and A. E. Raftery.  Strictly proper scoring rules, prediction, and estimation.  _J. Am. Stat. Assoc._ , 102:359–378, 2007.
  * Gneiting et al. (2005) T. Gneiting, A. E. Raftery, A. H. Westveld, and T. Goldman.  Calibrated probabilistic forecasting using ensemble model output statistics and minimum CRPS estimation.  _Mon. Weather Rev._ , 133:1098–1118, 2005.
  * Grimit et al. (2006) E. P. Grimit, T. Gneiting, V. J. Berrocal, and N. A. Johnson.  The continuous ranked probability score for circular variables and its application to mesoscale forecast ensemble verification.  _Q. J. R. Meteorol. Soc._ , 132:2925–2942, 2006.
  * Grünwald et al. (2019) P. Grünwald, R. de Heide, and W. Koolen.  Safe testing.  _Preprint, arXiv:1906.07801_, 2019.
  * Hansen et al. (2011) P. R. Hansen, A. Lunde, and J. M. Nason.  The model confidence set.  _Econometrica_ , 79(2):453–497, 2011.
  * Hao et al. (2024) Y. Hao, P. Grünwald, T. Lardy, L. Long, and R. Adams.  E-values for k-sample tests with exponential families.  _Preprint, arXiv:2303.00471_, 2024.
  * Henzi and Ziegel (2022) A. Henzi and J. F. Ziegel.  Valid sequential inference on probability forecast performance.  _Biometrika_ , 109(3):647–663, 2022.
  * Henzi et al. (2021) A. Henzi, J. F. Ziegel, and T. Gneiting.  Isotonic distributional regression.  _J. R. Stat. Soc. Ser. B: Stat. Methodol._ , 83:963–969, 2021.
  * Hothorn et al. (2014) T. Hothorn, T. Kneib, and P. Bühlmann.  Conditional transformation models.  _J. R. Stat. Soc. Ser. B: Stat. Methodol._ , 76:3–27, 2014\.
  * Howard et al. (2020) S. R. Howard, A. Ramdas, J. McAuliffe, and J. Sekhon.  Time-uniform chernoff bounds via nonnegative supermartingales.  _Probab. Surv._ , 17:257–317, 2020.
  * Howard et al. (2021) S. R. Howard, A. Ramdas, J. McAuliffe, and J. Sekhon.  Time-uniform, nonparametric, nonasymptotic confidence sequences.  _Ann. Stat._ , 49(2), 2021.
  * Lai et al. (2011) T. L. Lai, S. T. Gross, and D. B. Shen.  Evaluating probability forecasts.  _Ann. Stat._ , 39(5):2356–2382, 2011.
  * Masini et al. (2023) R. P. Masini, M. C. Medeiros, and E. F. Mendes.  Machine learning advances for time series forecasting.  _J. Econ. Surv._ , 37(1):76–111, 2023.
  * Matheson and Winkler (1976) J. E. Matheson and R. L. Winkler.  Scoring rules for continuous probability distributions.  _Manag. Sci._ , 22:1087–1096, 1976.
  * Pérez-Ortiz et al. (2023) M. F. Pérez-Ortiz, T. Lardy, R. de Heide, and P. Grünwald.  E-statistics, group invariance and anytime valid testing.  _Preprint, arXiv:2208.07610_, 2023.
  * Primo et al. (2024) C. Primo, B. Schulz, S. Lerch, and R. Hess.  Comparison of model output statistics and neural networks to postprocess wind gusts.  _Preprint, arXiv:2401.11896_, 2024.
  * Ramdas et al. (2020) A. Ramdas, J. Ruf, M. Larsson, and W. Koolen.  Admissible anytime-valid sequential inference must rely on nonnegative martingales.  _Preprint, arXiv:2009.03167_, 2020.
  * Ramdas et al. (2022) A. Ramdas, J. Ruf, M. Larsson, and W. M. Koolen.  Testing exchangeability: Fork-convexity, supermartingales and e-processes.  _Int. J. Approx. Reason._ , 141:83–109, 2022.
  * Ramdas et al. (2023) A. Ramdas, P. Grünwald, V. Vovk, and G. Shafer.  Game-theoretic statistics and safe anytime-valid inference.  _Stat. Sci._ , 38(4):576–601, 2023.
  * R Core Team (2024) R Core Team.  _R : A language and environment for statistical computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2024\.  URL https://www.R-project.org/.
  * Robbins (1970) H. Robbins.  Statistical methods related to the law of the iterated logarithm.  _Ann. Math. Stat._ , 41(5):1397–1409, 1970.
  * Robbins and Siegmund (1970) H. Robbins and D. Siegmund.  Boundary crossing probabilities for the wiener process and sample sums.  _Ann. Math. Stat._ , 41(5):1410–1429, 1970.
  * Rodriguez et al. (2021) A. Rodriguez, A. Tabassum, J. Cui, J. Xie, J. Ho, P. Agarwal, B. Adhikari, and B. A. Prakash.  Deepcovid: An operational deep learning-driven framework for explainable real-time Covid-19 forecasting.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 35(17):15393–15400, 2021.
  * Schulz and Lerch (2022) B. Schulz and S. Lerch.  Machine learning methods for postprocessing ensemble forecasts of wind gusts: A systematic comparison.  _Mon. Weather Rev._ , 150(1):235–257, 2022.
  * Shafer (2021) G. Shafer.  Testing by betting: A strategy for statistical and scientific communication.  _J. R. Stat. Soc. Ser. A: Stat. Soc._ , 184:407–431, 2021\.
  * Vannitsem et al. (2018) S. Vannitsem, D. S. Wilks, and J. Messner, editors.  _Statistical Postprocessing of Ensemble Forecasts_.  Elsevier, Amsterdam, Oxford, Cambridge, 2018.
  * Ville (1939) J. Ville.  _Étude Critique de la Notion de Collectif_.  PhD thesis, L’Université de Paris, 1939.
  * Vovk and Wang (2021) V. Vovk and R. Wang.  E-values: Calibration, combination and applications.  _Ann. Stat._ , 49(3):1736–1754, 2021.
  * Wang and Ramdas (2022) R. Wang and A. Ramdas.  False discovery rate control with e-values.  _J. R. Stat. Soc. Ser. B: Stat. Methodol._ , 84(3):822–852, 2022.
  * Waudby-Smith and Ramdas (2024) I. Waudby-Smith and A. Ramdas.  Estimating means of bounded random variables by betting.  _J. R. Stat. Soc. Ser. B: Stat. Methodol._ , 86(1):1–27, 2024.
  * Weron (2014) R. Weron.  Electricity price forecasting: A review of the state-of-the-art with a look into the future.  _Int. J. Forecast._ , 30(4):1030–1081, 2014.

##  Appendix A Proper scoring rules and consistent scoring functions

In this section, we review some important concepts for forecast evaluation. We
assume that we are interested in a real-valued unknown quantity and let
$\mathcal{P}(\mathbb{R})$ be the family of all probability distributions on
$\mathbb{R}$.

A _statistical functional_ is a map $T:\mathcal{P}(\mathbb{R})\to
2^{\mathbb{R}}$. Point-valued forecasts for functionals should be compared
using consistent scoring functions (Gneiting, 2011). A measurable map
$S:\mathbb{R}\times\mathbb{R}\rightarrow\bar{\mathbb{R}}$ is a _consistent
scoring function for the functional $T$ relative to the class
$\mathcal{P}\subseteq\mathcal{P}(\mathbb{R})$_ if for all $x\in\mathbb{R}$,
$F\in\mathcal{P}$, the integral $\int S(x,y)\,\mathrm{d}F(y)$ exists and

$\int S(t,y)\,\mathrm{d}F(y)\leq\int S(x,y)\,\mathrm{d}F(y),\quad\textrm{for
all }t\in T(F).$

The _Brier score_ (or _quadratic score_) $\textrm{BS}(x,y)=(x-y)^{2}$,
$x,y\in\mathbb{R}$, and the _quantile score_

$\textrm{QS}_{\tau}(x,y)=(\mathbbm{1}{\\{y\leq x\\}}-\tau)(x-y),\quad
x,y\in\mathbb{R},$ (20)

are main examples of consistent scoring functions for the mean functional, and
for the quantile at level $\tau\in(0,1)$, respectively. More generally, Ehm et
al. (2016) show that a scoring function $S$ is consistent for the quantile
functional at level $\tau\in(0,1)$ if and only if it is of the form

$S_{\tau}(x,y)=(\mathbbm{1}{\\{y\leq x\\}}-\tau)\big{(}g(x)-g(y)\big{)},\quad
x,y\in\mathbb{R},$

for some increasing function $g:\mathbb{R}\to\mathbb{R}$, where we obtain the
quantile score at (20) for $g$ equals the identity, and where the result
extends for the mean and the expectile.

Probabilistic forecasts incorporate the intrinsic uncertainty of future events
and comprise of a probability measure on the possible outcomes. Probabilistic
forecasts should be compared and evaluated using proper scoring rules
(Gneiting and Raftery, 2007). A _proper scoring rule_ is a function ${\rm
S}:\mathcal{P}\times\mathbb{R}\mapsto\bar{\mathbb{R}}$ for some class
$\mathcal{P}\subseteq\mathcal{P}(\mathbb{R})$ such that ${\rm S}(F,\cdot)$ is
measurable for any $F\in\mathcal{P}$, the integral $\int{\rm
S}(G,y)F(\,\mathrm{d}y)$ exists, and

$\int{\rm S}(F,y)\,\mathrm{d}F(y)\leq\int{\rm
S}(G,y)\,\mathrm{d}F(y),\quad\textrm{for all }F,G\in\mathcal{P}.$

That is, for a proper scoring rule ${\rm S}$, the expected score
$\mathbb{E}_{Y\sim F}{\rm S}(G,Y)$ is minimized over all distributional
forecasts $G\in\mathcal{P}$, if we correctly forecast the true distribution of
the random variable $Y$. The most commonly used proper scoring rule for the
evaluation of real-valued outcomes is the _continuous ranked probability
score_ (CRPS; Matheson and Winkler (1976)), which is defined for all
$F\in\mathcal{P}(\mathbb{R})$ and $y\in\mathbb{R}$ as

$\textrm{CRPS}(F,y)=\int(F(x)-\mathbbm{1}\\{y\leq x\\})^{2}\,\mathrm{d}x.$
(21)

The CRPS is a continuous extension of the Brier score for binary events and is
popular across application areas and methodological communities; see, e.g.,
Gneiting et al. (2005), Hothorn et al. (2014) and Henzi et al. (2021).

##  Appendix B Conditionally bounded CRPS and quantile score differences

The following results show that CRPS and quantile score differences are
conditionally bounded by the predictability of the forecasts.

######  Lemma B.1.

For $F_{1},F_{2}\in\mathcal{P}(\mathbb{R})$, the CRPS difference

$\delta(y)=\textrm{CRPS}(F_{1},y)-\textrm{CRPS}(F_{2},y)$

is bounded in $y\in\mathbb{R}$, and attains its minimum and maximum at a
crossing point of $F_{1}$ and $F_{2}$ (including $\\{-\infty,+\infty\\})$,
respectively.

###### Proof.

By the integral representation of the CRPS given at (21),

$\displaystyle\delta(y)$
$\displaystyle=\int\big{(}F_{1}(x)-\mathbbm{1}\\{y\leq
x\\}\big{)}^{2}-\big{(}F_{2}(x)-\mathbbm{1}\\{y\leq
x\\}\big{)}^{2}\,\mathrm{d}x$ $\displaystyle=\int
F_{1}(x)^{2}-F_{2}(x)^{2}\,\mathrm{d}x+2\int_{y}^{\infty}F_{2}(x)-F_{1}(x)\,\mathrm{d}x,$

where the first summand is independent of $y\in\mathbb{R}$. The function
$x\mapsto F_{2}(x)-F_{1}(x)$ is bounded in $[-1,1]$, with limit $0$ as
$x\to\pm\infty$ and changes sign at the crossing points of $F_{1}$ and
$F_{2}$. ∎

######  Lemma B.2.

For $x_{1},x_{2}\in\mathbb{R}$ and $\tau\in(0,1)$, the quantile loss
difference

$\delta(y)=\textrm{QS}_{\tau}(x_{1},y)-\textrm{QS}_{\tau}(x_{2},y)$

is bounded in $y\in\mathbb{R}$ with

$(\tau-1)\lvert x_{1}-x_{2}\rvert\leq\delta(y)\leq\tau\lvert
x_{1}-x_{2}\rvert.$

###### Proof.

Assume w.l.o.g. that $x_{1}\leq x_{2}$. Then, we have

$\displaystyle\delta(y)$ $\displaystyle=\big{(}\mathbbm{1}\\{y\leq
x_{1}\\}-\tau\big{)}(x_{1}-y)-\big{(}\mathbbm{1}\\{y\leq
x_{2}\\}-\tau\big{)}(x_{2}-y)$ $\displaystyle=y\mathbbm{1}\\{x_{1}<y\leq
x_{2}\\}-\tau(x_{1}-x_{2})+x_{1}\mathbbm{1}\\{y\leq
x_{1}\\}-x_{2}\mathbbm{1}\\{y\leq x_{2}\\}$

and the bounds follow directly by differentiating cases. ∎

##  Appendix C Sequential family-wise error rate and family-wise validity

Let $(E_{1,t})_{t\in\mathbb{N}},\dots,(E_{m,t})_{t\in\mathbb{N}}$ be
e-processes for some hypotheses
$\mathcal{H}_{1},\dots,\mathcal{H}_{m}\subseteq\mathfrak{B}(\Omega)$.

Following Vovk and Wang (2021), we call $f$ a _symmetric e-merging function_
if it is invariant with respect to permutations of its arguments and if
$f(E_{i},i\in I)$ is an e-value for any family of e-values $\\{E_{i}\mid i\in
I\\}$. The most important e-merging function is the arithmetic mean as it
essentially dominates any symmetric e-merging function (Vovk and Wang, 2021,
Proposition 3.1.). Let $f$ be any symmetric e-merging function. Then, by the
closure principle (Vovk and Wang, 2021, Appendix A.3.), for any
$t\in\mathbb{N}$,

$E^{\star}_{i,t}=\min\\{f(E_{j,t},j\in I)\mid I\subseteq\\{1,\dots,m\\},i\in
I\\},\quad i=1,\dots,m,$ (22)

are family-wise valid e-variables for the hypotheses
$\mathcal{H}_{1},\dots,\mathcal{H}_{m}$. That is, there exists a dominating
family of e-variables
$\\{E^{\mathbb{Q}}\mid\mathbb{Q}\in\mathfrak{B}(\Omega)\\}$, with
$E^{\mathbb{Q}}\geq E^{\star}_{i,t}$ for any $i=1,\dots,m$ and
$\mathbb{Q}\in\mathcal{H}_{i}$.

If we adjust the e-processes
$(E_{1,t})_{t\in\mathbb{N}},\dots,(E_{m,t})_{t\in\mathbb{N}}$ by (22) with
respect to the same e-merging function $f$ for all $t\in\mathbb{N}$, then it
follows directly that
$(E^{\star}_{1,t})_{t\in\mathbb{N}},\dots,(E^{\star}_{m,t})_{t\in\mathbb{N}}$
are family-wise valid e-processes for the hypotheses
$\mathcal{H}_{1},\dots,\mathcal{H}_{m}$. Importantly, the corresponding
sequential test given at (8) control the _family-wise error rate (FWER)_
uniformly in time, that is, for any $I\subseteq\\{1,\dots,m\\}$ and
$\mathbb{Q}\in\bigcap_{i\in I}\mathcal{H}_{i}$,

$\textrm{FWER}_{I}(\mathbb{Q})=\mathbb{Q}\left(\exists
t\in\mathbb{N}:E_{i,t}^{\star}\geq 1/\alpha\textrm{ for some }i\in
I\right)\leq\alpha,\quad\textrm{for all }\alpha\in(0,1).$ (23)

To verify (23), consider a dominating family of e-processes
$\\{(E_{t}^{\mathbb{Q}})_{t\in\mathbb{N}}\mid\mathbb{Q}\in\mathfrak{B}(\Omega)\\}$.
Then, for any $\alpha\in(0,1)$, $I\subseteq\\{1,\dots,m\\}$ and
$\mathbb{Q}\in\bigcap_{i\in I}\mathcal{H}_{i}$,

$\displaystyle\textrm{FWER}_{I}(\mathbb{Q})\leq\mathbb{Q}\left(\exists
t\in\mathbb{N}:E_{t}^{\mathbb{Q}}\geq 1/\alpha\right)\leq\alpha.$

##  Appendix D SMCSs by FDR control

For $i=1,\dots,m$, let $(\Psi_{i,t})_{t\in\mathbb{N}}$ be a sequential testing
procedure for some hypothesis $\mathcal{H}_{i}\subseteq\mathfrak{B}(\Omega)$.
We let $R_{t}=\sum_{i=1}^{m}\Psi_{i,t}$ be the number of all rejected
hypotheses (discoveries) at $t\in\mathbb{N}$. For
$\mathbb{Q}\in\bigcup_{i=1}^{m}\mathcal{H}_{i}$, we let
$F_{t}(\mathbb{Q})=\sum_{i\in I_{\mathbb{Q}}}\Psi_{i,t}$ be the number of all
true null hypotheses that are rejected at time $t$. The _false discovery rate
(FDR)_ at $t\in\mathbb{N}$ of the given testing procedures is defined as

$\textrm{FDR}_{t}=\sup_{\mathbb{Q}\in\bigcup_{i=1}^{m}\mathcal{H}_{i}}\mathbb{E}_{\mathbb{Q}}\left[\frac{F_{t}(\mathbb{Q})}{\max\\{1,R_{t}\\}}\right].$

We say that the sequential tests
$(\Psi_{1,t})_{t\in\mathbb{N}},\dots,(\Psi_{m,t})_{t\in\mathbb{N}}$ _time-
uniformly control the false discovery rate at level_ $\alpha\in(0,1)$ if, for
all (possibly infinite) stopping times $\tau$,

$\textrm{FDR}_{\tau}\leq\alpha.$ (24)

Wang and Ramdas (2022) study FDR control with e-values in a nonsequential
setting. Looking carefully at the proof of their Theorem 2, we see that their
arguments directly extend to e-processes to obtain sequential tests that
satisfy (24) as explained in the following.

Let $(E_{1,t})_{t\in\mathbb{N}},\dots,(E_{m,t})_{t\in\mathbb{N}}$ be
e-processes for the hypotheses
$\mathcal{H}_{1},\dots,\mathcal{H}_{m}\subseteq\mathfrak{B}(\Omega)$.
Following Wang and Ramdas (2022), we let $E_{[1],t}\geq\dots\geq E_{[m],t}$ be
the order statistics of $E_{1,t},\dots,E_{m,t}$ in decreasing order. For a
given $\alpha\in(0,1)$, we define

$i^{\star}_{t}=\max\left\\{i\in\\{1,\dots,m\\}\;\Big{|}\;\frac{iE_{[i],t}}{m}\geq\frac{1}{\alpha}\right\\},\quad
t\in\mathbb{N}.$

Then, we have a time-uniform FDR control at level $\alpha$ for the sequential
tests

$\Psi_{i,t}=(\mathbbm{1}\\{E_{i,t}\geq E_{[i_{t}^{\star}],t}\\}),\quad
i=1,\dots,m,t\in\mathbb{N}.$ (25)

As our construction for the strong and uniformly weak hypothesis builds on
e-processes, we can use the sequential tests at (25) with respect to the
e-processes given at (9). That is, at each $t\in\mathbb{N}$, we reject the
$i^{\star}_{t}$ hypotheses with the largest values of the corresponding
e-processes, and obtain the alternative SMCS

$\check{\mathcal{M}}_{t}=\left\\{i\in\mathcal{M}_{0}\mid
E_{i\cdot,t}<E_{[i_{t}^{\star}]\cdot,t}\right\\}.$

Importantly, this sequence does not necessarily satisfy the originally imposed
coverage requirement at (4) or (5) anymore. In contrast, it bounds the
expected fraction of wrongly rejected models by the total number of rejected
models at any random time by $\alpha$, that is, for any
$\mathbb{Q}\in\mathfrak{B}(\Omega)$ and any random stopping time $\tau$, we
have

$\mathbb{E}_{\mathbb{Q}}\left[\frac{\lvert\check{\mathcal{M}}^{c}_{\tau}\cap\mathcal{M}^{\star}_{\tau}\rvert}{\max\\{1,\lvert\check{\mathcal{M}}_{\tau}^{c}\rvert\\}}\right]\leq\alpha,$

where
$\check{\mathcal{M}}^{c}_{\tau}=\\{1,\dots,m\\}\setminus\check{\mathcal{M}}_{\tau}$
and $(\mathcal{M}^{\star})_{t\in\mathbb{N}}$ is given by (1) or (2).

##  Appendix E Proof of Lemma 3.5

###### Proof.

For $t\in\mathbb{N}$, we have
$C_{t,1-\alpha}=\\{\boldsymbol{X}\in\mathbb{R}_{0}^{m\times m}\mid\lVert
f_{t}(\boldsymbol{X})\rVert_{1}\leq m(m-1)/\alpha\\}$, for

$f_{t}:\mathbb{R}_{0}^{m\times m}\to\mathbb{R}_{0}^{m\times
m},\quad\boldsymbol{X}\mapsto\left(f_{ij,t}(X_{ij})\right)_{i\neq
j}=\left(\exp\left\\{\lambda_{ij}t\hat{\Delta}_{ij}-\lambda_{ij}tx-\psi_{E,c_{ij}}(\lambda_{ij})V_{ij,t}\right\\}\right)_{i\neq
j},$

and where $\lVert\boldsymbol{X}\rVert_{1}=\sum_{i,j=1}^{m}\lvert X_{ij}\rvert$
for $\boldsymbol{X}\in\mathbb{R}^{m\times m}$. For any $i\neq j$, $f_{ij,t}$
is a nonnegative and nonincreasing convex function. Thus, the second claim
follows immediately since $\boldsymbol{A}\leq\boldsymbol{B}$ implies
$f_{t}(\boldsymbol{A})\geq f_{t}(\boldsymbol{B})$. Let
$\boldsymbol{A},\boldsymbol{B}\in C_{t,1-\alpha}$, then, for any
$\lambda\in[0,1]$,

$\displaystyle\lVert
f_{t}\big{(}\lambda\boldsymbol{A}+(1-\lambda)\boldsymbol{B}\big{)}\rVert_{1}\leq\lVert\lambda
f_{t}(\boldsymbol{A})+(1-\lambda)f_{t}(\boldsymbol{B})\rVert_{1}\leq\lambda\lVert
f_{t}(\boldsymbol{A})\rVert_{1}+(1-\lambda)\lVert
f_{t}(\boldsymbol{B})\rVert_{1}\leq\frac{m(m-1)}{\alpha},$

where the first inequality holds by convexity and the nonnegativity of
$f_{t}$, the second inequality by the axioms of a norm, and the third
inequality by the fact that $\boldsymbol{A},\boldsymbol{B}\in C_{t,1-\alpha}$.
∎

##  Appendix F Additional material for the Covid case study

###  F.1 Documentation of the models

The first forecasting model _baseline_ is a naive model created by the
Covid-19 forecast hub team as a comparison benchmark. Baseline simply outputs
the outcome of the most recent week as a weekly median forecast and forms a
predictive distribution around this median prediction based on the past weekly
incidences. For more information, we refer to the supplementary material of
Cramer et al. (2022b). Secondly, we consider the model _MOBS-gleam_ , a global
epidemic and mobility structured meta-population model (GLEaM), developed by
the MOBS lab. It consists of a population-mobility layer, which uses
transportation data to model the population mobility though a stochastic
process, and of a disease model, which classifies each person in a population
as susceptible, latent, infectious, or recovered, and estimates the rate in
which a susceptible individual contracts the disease. For more details, we
refer to the original MOBS lab publications by Balcan et al. (2009) and Balcan
et al. (2010). The third model _PSI-draft_ was developed by the Predictive
Science Inc. It is a SEIRX model and estimates the reproductive number $R_{0}$
by classifying each individual of a population into susceptible (S), exposed
(E), infected (I), removed (R), or quarantined (X). The model stratifies the
population by age and disease severity, and incorporates the social
restrictions imposed over the course of the pandemic, see the GitHub page of
the Covid-19 forecast hub (https://github.com/reichlab/covid19-forecast-
hub/tree/master/data-processed/PSI-DRAFT) or the Zoltar forecast archive
(https://zoltardata.com/model/254) for more information. Further, we consider
the neural network based model _GT-deep_ , which includes some auto-regressive
terms to ensure temporal consistency across consecutive predictions. GT-deep
was the first purely data-driven model to be included in the Covid-19 forecast
hub, in the sense that it makes no epidemiological assumptions about the
spread of the pandemic. More detailed information can be found in the original
publication by Rodriguez et al. (2021). Finally, we consider two models that
made their predictions based on previously issued forecasts from other models:
Firstly, _ensemble_ , formed by the Covid-19 forecast hub team, which combines
the latest predictions of all models that submitted a full set of forecasts at
all $23$ quantile values. Whereas ensemble outputted the average forecast of
all individual models at each quantile level for the first few weeks of the
pandemic, it switched to the corresponding median forecast later on. As
discussed in Cramer et al. (2022b), ensemble was the model used by the CDC in
its public communications about the pandemic. Finally, _CDC-ensemble_ ,
developed jointly by the Covid-19 forecast hub and CDC, is very similar to
ensemble, but instead of the median of all models (at each quantile value), it
considers only the 10 best models, measured by the average weighted interval
score over the 12 most recent weeks, see Cramer et al. (2022a).

###  F.2 Comparison to the naive betting approach

Figure 7 shows the resulting SMCSs with respect to the naive approach, where
$K_{ij,t}=2$ for all $i,j=1,\dots,m$ and $t\in\mathbb{N}$, compare to Section
5.1.2.

Figure 7: The resulting SMCSs with respect to the naive choice $K_{ij,t}=2$
for all $i,j\leq m$ and $t\in\mathbb{N}$ with confidence level $\alpha=0.1$.
Whereas, there are no significant differences to the adaptive betting scheme
proposed in Section 5.1.2 for $\tau=0.3,0.5,0.9$, the adaptive approach is
clearly superior for $\tau=0.1$.

##  Appendix G Station categories for the wind gust case study

Figure 8: Spatial distribution of the station categories over Germany.