Evaluation of time series models under non-stationarity with application to the comparison of regional climate models

TL;DR

This paper introduces a methodology for evaluating time series models under non-stationarity by using moving scores based on changepoint analysis, enabling more accurate model comparison in changing conditions.

Contribution

It proposes a novel approach to model evaluation that accounts for non-stationarity through moving windows and changepoint analysis, improving fairness and accuracy.

Findings

Methodology effectively captures non-stationary behavior.

Application to climate models demonstrates practical utility.

Simulation and case study validate the approach.

Abstract

Different disciplines pursue the aim to develop models which characterize certain phenomena as accurately as possible. Climatology is a prime example, where the temporal evolution of the climate is modeled. In order to compare and improve different models, methodology for a fair model evaluation is indispensable. As models and forecasts of a phenomenon are usually associated with uncertainty, proper scoring rules, which are tools that account for this kind of uncertainty, are an adequate choice for model evaluation. However, under the presence of non-stationarity, such a model evaluation becomes challenging, as the characteristics of the phenomenon of interest change. We provide methodology for model evaluation in the context of non-stationary time series. Our methodology assumes stationarity of the time series in shorter moving time windows. These moving windows, which are selected based on a changepoint analysis, are used to characterize the uncertainty of the phenomenon/model for the corresponding time instances. This leads to the concept of moving scores allowing for a temporal assessment of the model performance. The merits of the proposed methodology are illustrated in a simulation and a case study.

Equations

$$s_{\mathrm{SE}}(F,y):=(E_F[X]-y{)}^2$$

$$s_{\mathrm{CRPS}}(F,y):={\int }_{-\infty }^{\infty }{\left(F(z)-\text{mathds}1\left\{z\ge y\right\}\right)}^2\text{d}z=E_F\left|X-y\right|-\frac{1}{2}{E}_F\left|X-\tilde{X}\right|,$$

$$s_{\mathrm{SE}}(\mathbf{x},y)={\left(\frac{1}{n}\sum _{j=1}^n{x}_j-y\right)}^2,$$

$$s_{\mathrm{CRPS}}(\mathbf{x},y)=\frac{1}{n}\sum _{j=1}^n\left|x_j-y\right|-\frac{1}{2n^2}\sum _{j=1}^n\sum _{k=1}^n\left|x_j-x_k\right|.$$

$$\sum _{j=0}^m\mathcal{C}({\mathbf{y}}_{({\tau }_j+1):{\tau }_{j+1}})+\kappa m.$$

$${\hat{\mu }}_j=\frac{1}{{\tau }_{j+1}-{\tau }_j}\sum _{i={\tau }_j+1}^{{\tau }_{j+1}}{y}_i\text{and}{\hat{\sigma }}_j^2=\frac{1}{{\tau }_{j+1}-{\tau }_j}\sum _{i={\tau }_j+1}^{{\tau }_{j+1}}{\left(y_i-{\hat{\mu }}_j\right)}^2,$$

$$\mathcal{C}({\mathbf{y}}_{({\tau }_j+1):{\tau }_{j+1}}):=({\tau }_{j+1}-{\tau }_j)\left\{\ln \left[\frac{2\pi }{{\tau }_{j+1}-{\tau }_j}\sum _{i={\tau }_j+1}^{{\tau }_{j+1}}{\left(y_i-\frac{{\sum }_{k={\tau }_j+1}^{{\tau }_{j+1}}{y}_k}{{\tau }_{j+1}-{\tau }_j}\right)}^2\right]+1\right\}.$$

$$\mathcal{W}(t)=\{t-{\delta }_t^{-},\dots ,t,\dots ,t+{\delta }_t^{+}\}\subset \{1,\dots ,N\},$$

$$s({\mathbf{x}}_{\mathcal{W}(t)},y_t)$$

$${\delta }_t^{\mathrm{OF}}:=\{\begin{array}{ll}t-1,& \text{for }t=1,\dots ,{\delta }^{\mathrm{OF}},\\ {\delta }^{\mathrm{OF}},& \text{for }t=1+{\delta }^{\mathrm{OF}},\dots ,N-{\delta }^{\mathrm{OF}},\\ N-t,& \text{for }t=N-{\delta }^{\mathrm{OF}}+1,\dots ,N,\end{array}$$

$${\mathcal{W}}^{\mathrm{OF}}(t):=\left\{t-{\delta }_t^{\mathrm{OF}},\dots ,t,\dots ,t+{\delta }_t^{\mathrm{OF}}\right\},$$

$${\mathcal{W}}^{\mathrm{OV}}(t):=\left\{t-{\delta }_t^{\mathrm{OV}},\dots ,t,\dots ,t+{\delta }_t^{\mathrm{OV}}\right\},$$

$${\delta }_t^{\mathrm{OV}}:=\{\begin{array}{ll}t-1,& \text{for }t=1,\dots ,\lfloor {\gamma }_0\rfloor ,\\ \lfloor \left(\varsigma (t)-1\right)/2\rfloor ,& \text{for }t=\lceil {\gamma }_0\rceil ,\dots ,\lfloor {\gamma }_{m+1}\rfloor ,\\ N-t,& \text{for }t=\lceil {\gamma }_{m+1}\rceil ,\dots ,N.\end{array}$$

$${\mathcal{W}}^{\mathrm{DV}}(t):=\left\{{\tau }_j+1,\dots ,{\tau }_{j+1}\right\},\text{for }t={\tau }_j+1,\dots ,{\tau }_{j+1},j=0,\dots ,m.$$

$$Y_t$$

$$X_t^k$$

$$Y_t$$

$$X_t^1$$

$$X_t^2$$

$$X_t^3$$

$$X_t^4$$

$$X_t^5$$

$${\mathbf{\mu }}^0=(0,1,0),$$

$${\mathbf{\mu }}^1={\mathbf{\mu }}^0,$$

$${\mathbf{\mu }}^2=(0.25,0.25,0.25),$$

$${\mathbf{\mu }}^3={\mathbf{\mu }}^0,$$

$${\mathbf{\mu }}^4={\mathbf{\mu }}^2,$$

$${\mathbf{\mu }}^5=(0.1,0.9,0.1),$$

$$Y_t$$

$$X_t^1$$

$$X_t^2$$

$$X_t^3$$

$$X_t^4$$

$$X_t^5$$

$${\mathbf{\mu }}^0=(0,1/3,2)/N,$$

$${\mathbf{\mu }}^1={\mathbf{\mu }}^0,$$

$${\mathbf{\mu }}^2=(0,1/3,1.9)/N,$$

$${\mathbf{\mu }}^3={\mathbf{\mu }}^0,$$

$${\mathbf{\mu }}^4={\mathbf{\mu }}^2,$$

$${\mathbf{\mu }}^5={\mathbf{\mu }}^2,$$

$$Y_t$$

Full text

Keywords: Changepoints; Climate models; ENSEMBLES; Moving scores; Performance measures, Proper scoring rules

1 Introduction

Ever since $1990$, the Intergovernmental Panel on Climate Change (IPCC) has published regular reports on the current state of knowledge on climate change. Its latest Assessment Report (AR5, see IPCC,, 2014, for the synthesis report) consolidates our current understanding of climate change, its causes and impacts, and discusses possible adaptation and mitigation strategies. Climate models which provide a mathematical description of certain processes in the Earth’s climate system are the main tool to learn about the possible future changes of the climate.

A multitude of different climate models is used to simulate future states of the climate as described by variables such as surface temperature and precipitation, amongst others. In order to evaluate the accuracy of different climate models, simulations from these models are compared to historical observations. As Knutti et al., (2010) note, no general all-purpose metric has been found that unambiguously identifies a best model. Furthermore, different metrics may produce different model rankings (e.g. Gleckler et al.,, 2008).

One way to compare model simulations to observations is to perform the comparison separately for each time unit (and each spatial unit). Such a point-wise comparison (e.g. based on one of the metrics discussed in Hyndman and Koehler,, 2006) is obviously the most straightforward approach to model evaluation. However, this approach is not adequate for evaluating climate models based on a daily or sub-daily temporal resolution. Climate models aim to model the long-term evolution of the climate, and can not provide accurate simulations/predictions on a daily basis. Rather, they intend to model/simulate the characteristics of the climate for longer time periods, up to several decades. In other words, the climate model output for a specific time instance is associated with uncertainty which is neglected in a point-wise evaluation.

In this setting, three questions arise:

How can we assess the model/prediction distribution, if all we have is one time series of realizations from the model/deterministic forecasts? 2. 2.

How can we assess the distribution of the observed quantity, if all we have is one time series of realizations of that quantity? 3. 3.

Can we assume one and the same distribution for different time instances, that is, is the distribution stationary?

To answer these questions we cling to our application of climate model evaluation. Under the presence of climate change and seasonality, it is obvious that we have to negate the third question. This makes it more difficult to answer the first two questions. Our answer assumes that the characteristics of the variables of interest change only gradually and can be considered as (approximately) stationary for short time windows. Then, for each of these time windows, we can construct empirical distributions based on the corresponding realizations. These empirical distributions are subsequently evaluated using the framework of proper scoring rules (see e.g. Gneiting and Raftery,, 2007) where an observed quantity is compared against a predicted distribution.

The (moving) time windows, for which we assume stationarity, are selected using a changepoint detection algorithm (pruned exact linear time (PELT) algorithm; Killick et al.,, 2012). We propose and compare three different window selection strategies. Based on the samples corresponding to these (moving) windows we compute time series of (moving) scores. This allows us to assess the model performance over time. The introduced evaluation technique is not restricted to the evaluation of climate models. Hence, we introduce it in a general setting. To learn how the approach operates in different settings we conduct a simulation study (Section 4) covering a wide range of cases which are relevant in practice. As a case study (Section 5) we evaluate daily mean surface temperature output of four regional climate models (RCMs) on a fine resolution grid covering Europe (see the ENSEMBLES project: van der Linden and Mitchell,, 2009). This article and all its results are based on Chapter 5 of Erhardt, (2017).

In the ENSEMBLES project, European climate research institutes jointly compiled a large data set of RCM simulations. One goal of this project was it to evaluate and compare the different models (see e.g. Kjellström et al.,, 2010; Lorenz and Jacob,, 2010). While Lorenz and Jacob, (2010) compare linear trends in seasonal temperatures, Kjellström et al., (2010) attempt to compare full probability distributions of temperature and precipitation against corresponding distributions from gridded observations (E-OBS, Haylock et al.,, 2008). They construct empirical estimates of probability density functions by binning the data into a certain number of bins (Perkins et al.,, 2007). As Kjellström et al., (2010) note, this approach requires several subjective choices.

Further ENSEMBLES evaluation studies include Landgren et al., (2014) who compare local temperature and precipitation model output against different data products. They provide a model ranking based on the root mean square deviation between model output and reference data, and a measure for differences in inter-annual variability. In a study of RCM output for Canada, Eum et al., (2012) investigate seasonal model performance based on five attributes: (i) relative absolute mean error on a daily time scale, (ii) differences in the annual variability of monthly means, (iii) differences in the spatial pattern of mean values in a certain region, (iv) discrepancy between $0.1$ and $0.9$ quantiles (“extremes”) of daily observations and model output, and (v) differences in long-term linear trends. This is an attempt to evaluate different features/moments of probabilistic output rather than considering the full distribution which comes along with the difficulty of judging the relative importance of each attribute. Alternatively, the approach proposed here considers the full distribution in time windows which are chosen such that seasonal differences are evaluated in an appropriate manner.

The remainder of the article is organized as follows. In Section 2 we provide the background information on proper scoring rules and the changepoint detection algorithm (PELT algorithm; Killick et al.,, 2012) for the specification of the moving time windows. Section 3 describes the new moving score methodology, whose properties are investigated in the simulation study in Section 4. Section 5 provides the case study on RCMs, showing how moving scores can be applied in practice. Our conclusions are summarized in Section 6.

2 Preliminaries

2.1 Proper scoring rules

To outline the theory behind proper scoring rules (see Gneiting and Raftery,, 2007) we consider the following setup: Let $\mathcal{F}$ be a convex class of probability measures on a sample space $\Omega$. We consider an (observed) phenomenon with the random outcome $Y$ with (unknown) distribution $G\in\mathcal{F}$ and realization $y\in\Omega$. Moreover, we consider models/predictions for $Y$ given through a random variable $X$ with (modeled) distribution $F\in\mathcal{F}$ and realization $x\in\Omega$.

Then a scoring rule $s:\mathcal{F}\times\Omega\rightarrow\mathbb{R}$ is a (negatively oriented) proper scoring rule if $S(G,G)\leq S(F,G)$ for all $F,G\in\mathcal{F}$, where we define the expectation of the score $s(F,Y)$ as $S(F,G)\coloneqq\operatorname{\mathbb{E}}_{G}\left[s(F,Y)\right]$. Hence, a scoring rule is proper if, in expectation, the random score $s(F,Y)$ is optimized if our model/prediction $F$ equals the true distribution $G$ of $Y$. Popular examples of proper scoring rules are the Squared Error (SE)

[TABLE]

and the Continuous Ranked Probability Score (CRPS)

[TABLE]

where $\operatorname{\mathds{1}}\left\{z\geq y\right\}$ equals $1$ if $z\geq y$, otherwise [math], and $\widetilde{X}\sim F$ is an independent copy of $X\sim F$.

In practice, we often do not know the distribution $F$ explicitly. Instead, we have a sample ${\bm{x}}\coloneqq(x_{1},\ldots,x_{n})$, $n\in\mathbb{N}$, from $F$. In that case, we consider sample versions of the scoring rules from above. The sample version of the SE in (1) is

[TABLE]

and the sample version of the CRPS in (2) is

[TABLE]

2.2 Detection of multiple changepoints using the PELT method

To introduce the Pruned Exact Linear Time (PELT) algorithm (see Killick et al.,, 2012), we consider a time series ${\bm{y}}_{1:N}\coloneqq(y_{1},\ldots,y_{N})$. Generally speaking, a time instance $\tau\in\{1,\ldots,N-1\}$ is considered a changepoint, if the statistical properties of the sub-series ${\bm{y}}_{1:\tau}=(y_{1},\ldots,y_{\tau})$ and ${\bm{y}}_{(\tau+1):N}=(y_{\tau+1},\ldots,y_{N})$ differ. To explain how multiple changepoints are detected, we first introduce further notation. Let $\bm{\tau}_{1:m}\coloneqq(\tau_{1},\ldots,\tau_{m})$ denote the ordered sequence of $m\in\{0,\ldots,N-1\}$ changepoints of ${\bm{y}}_{1:N}$, where $\tau_{j}\in\mathbb{N}$, $1\leq\tau_{j}\leq N-1$, $j=1,\ldots,m$. Defining $\tau_{0}\coloneqq 0$ and $\tau_{m+1}\coloneqq N$, the changepoints $\bm{\tau}_{1:m}$ split the sequence ${\bm{y}}_{1:N}$ into the $m+1$ segments ${\bm{y}}_{(\tau_{j}+1):\tau_{j+1}}$ with segment lengths $(\tau_{j+1}-\tau_{j})$, $j=0,\ldots,m$. We call $\bm{\tau}_{0:(m+1)}\coloneqq(\tau_{0},\bm{\tau}_{1:m},\tau_{m+1})$ an $(m+1)$-segmentation of ${\bm{y}}_{1:N}$.

The PELT algorithm detects multiple changepoints by minimizing a target function

[TABLE]

The target function is a sum of cost functions $\mathcal{C}:\mathbb{R}^{n}\to\mathbb{R}$, $n\in\mathbb{N}$, assigning a cost to the segment ${\bm{y}}_{(\tau_{j}+1):\tau_{j+1}}$ for all $j=0,\ldots,m$, and a penalty term $\kappa m$ which is supposed to prevent overfitting by penalizing the number of changepoints $m$.

As we are particularly interested in detecting changes in mean and variance, we assume that the time series observations $y_{1},\ldots,y_{N}$ come from a normal distribution $\mathcal{N}(\mu_{j},\sigma_{j}^{2})$, where either the (unknown) mean $\mu_{j}$ and/or the (unknown) variance $\sigma_{j}^{2}$ change after certain (unknown) time instances (changepoints) $\bm{\tau}_{1:m}$. In this case, we can select the cost function $\mathcal{C}$ as twice the negative log-likelihood corresponding to a normal distribution, where the unknown parameters $\mu_{j}$ and $\sigma_{j}^{2}$ are replaced by their maximum likelihood estimators

[TABLE]

respectively. Hence, the cost $\mathcal{C}({\bm{y}}_{(\tau_{j}+1):\tau_{j+1}})$ of a segment ${\bm{y}}_{(\tau_{j}+1):\tau_{j+1}}$ of an arbitrary $(m+1)$-segmentation $\bm{\tau}_{0:(m+1)}$ of the ordered sequence ${\bm{y}}_{1:N}$ is given by

[TABLE]

Note, that the estimation of the variance demands a minimum segment length $(\tau_{j+1}-\tau_{j})$ of $2$ for all segments ${\bm{y}}_{(\tau_{j}+1):\tau_{j+1}}$, $j=0,\ldots,m$. Common choices for $\kappa$ are $\kappa=2p$ (cp. Akaike information criterion) or $\kappa=p\ln(N)$ (cp. Bayesian information criterion), which depend on the number $p$ of additional parameters needed per additional segment.

The PELT algorithm itself builds on the Optimal Partitioning (OP) algorithm introduced by Jackson et al., (2005). It minimizes (5) based on a recursive detection of the changepoints of the sub-series ${\bm{y}}_{1:s}$, $s\leq N$. Killick et al., (2012) found that the performance of the OP algorithm can be improved by reducing the set of candidate changepoints in the recursion (so called pruning), in order to avoid irrelevant computations.

3 Methodology

In this section we propose a methodology for the (proper) evaluation (see Section 2.1) of time series models in a non-stationary context. In comparison to a purely point-wise evaluation where a model output $x_{t}$ at time point $t$ is compared against the corresponding observation $y_{t}$, our approach also considers higher order structures of the time series. Subjective choices regarding the importance of different model characteristics are avoided. As the approach is designed to deal with seasonality, it does not require a separate consideration of different selected seasons.

In the following, we consider the following setup: Let $\mathcal{F}$ be a convex class of probability measures on a sample space $\Omega$. We consider an (observed) phenomenon with the random outcomes $Y_{t}$ with (unknown) distributions $G_{t}\in\mathcal{F}$ and realizations $y_{t}\in\Omega$ ($t=1,\ldots,N$). Further, we consider a model/prediction for $Y_{t}$, $t=1,\ldots,N$, given by random variables $X_{t}$ with (modeled) distributions $F_{t}\in\mathcal{F}$ and realizations $x_{t}\in\Omega$ ($t=1,\ldots,N$).

3.1 Discussion of naive evaluation approaches

Evaluation under a stationarity (ST) assumption

Let us first assume that $Y_{t}$ and $X_{t}$, $t=1,\ldots,N$, are stationary. Hence, $G_{t}=G$ and $F_{t}=F$ for all $t=1,\ldots,N$. Under this assumption, we can evaluate the model $X_{t}$ for $Y_{t}$, using scores $s^{\operatorname{ST}}(F,y_{t})$, $t=1,\ldots,N$. To reflect that they are calculated based on a stationarity assumption, we call them ST-scores. Their corresponding sample version is given by $s^{\operatorname{ST}}({\bm{x}},y_{t})$, $t=1,\ldots,N$, where ${\bm{x}}=(x_{1},\ldots,x_{N})$. In practice the stationarity assumption is violated in most cases. Hence, ST-score based evaluation is inappropriate. A possible exception are very short time series data, such as annual series, where non-stationarity cannot be estimated reliably.

Point-wise (PW) evaluation

Let us now consider the non-stationary case. Then, it generally holds that $G_{s}\neq G_{t}$ and $F_{s}\neq F_{t}$ for $s\neq t$. A naive approach for model evaluation is to use sample scores $s^{\operatorname{PW}}(x_{t},y_{t})$, $t=1,\ldots,N$. As they compare the realizations of $Y_{t}$ and $X_{t}$ point-wise (separately for each $t=1,\ldots,N$), we call them PW-scores. From Equations (3) and (4) we obtain $s_{\operatorname{SE}}^{\operatorname{PW}}(x_{t},y_{t})=\left(x_{t}-y_{t}\right)^{2}$ and $s_{\operatorname{CRPS}}^{\operatorname{PW}}(x_{t},y_{t})=\left|x_{t}-y_{t}\right|$ for all $t=1,\ldots,N$. For a comprehensive model evaluation, PW-scores are also not satisfying, since they do not account for higher order structures of the observed/modeled phenomenon. They evaluate only how similar the phenomenon and the model/prediction behave in terms of their (temporarily varying) mean. Such point-wise model evaluation completely neglects differences in higher order moments (e.g. variance). A PW-score treats a discrepancy in $x_{t}$ and $y_{t}$ that occurs due to a falsely specified model/prediction mean in the same way as a discrepancy that occurs due to a high uncertainty of the phenomenon.

3.2 Moving scores

To deal with non-stationarity and to consider higher order structures of the observed phenomenon in the evaluation we suggest the following. We assume that the phenomenon $Y_{t}$ and the corresponding model/prediction $X_{t}$ are (approximately) stationary for short time intervals. Then, for each time instance $t=1,\ldots,N$ (window location), we select integers $\delta_{t}^{-},\delta_{t}^{+}\in\{0,\ldots,N-1\}$ which determine the width of a moving (time) window

[TABLE]

such that it completely lies within the observation period and such that $Y_{s}$ and $X_{s}$ are (approximately) stationary for all $s\in\mathcal{W}(t)$. Thus, we assume that for $s\in\mathcal{W}(t)$, $Y_{s}$ and $X_{s}$ are distributed according to distributions $G_{t}$ and $F_{t}$ (depending on the window location $t$), respectively. Hence, in order to evaluate the model/prediction $X_{t}$ for the phenomenon $Y_{t}$, for a specific time instance $t=1,\ldots,N$, we substitute the theoretical scores $s(F_{t},y_{t})$ by sample scores

[TABLE]

based on the sub-samples ${\bm{x}}_{\mathcal{W}(t)}\coloneqq\{x_{s}:s\in\mathcal{W}(t)\}$. Since the scores (7) are based on moving time windows, we call them moving scores. In contrast to the theoretical scores, the moving scores (7) are empirical scores, since they are sample-based. Same holds for the PW- and ST-scores introduced in Section 3.1. In order to rank different models, we can consider the score averages over $t=1,\ldots,N$.

Bearing in mind the discussion of Section 3.1 it becomes clear that the moving score methodology is a compromise between a point-wise evaluation (which does not account for higher order structures) and an evaluation which ignores non-stationarity. There is a trade-off between small and large moving windows $\mathcal{W}(t)$: In order to not violate the stationarity assumption, we have to keep the moving windows small enough. Considering the formula (4) for the computation of sample CRPS, we see that its usage in a moving window based evaluation results in a computation time which grows quadratically with increasing window width $|\mathcal{W}(t)|$. Hence, if the full data set consists of multiple long time series for several models (and possibly many different spatial locations), a moving window based evaluation may become infeasible if the moving windows are too large111For such cases, see Hersbach, (2000) who proposes an equivalent yet more computationally efficient computation of the sample version of the CRPS.. However, if the moving windows are too small (small sample size), the samples ${\bm{x}}_{\mathcal{W}(t)}$ coming from these small windows $\mathcal{W}(t)$ can only achieve an inaccurate approximation of the true distributions $F_{t}$. We see that there are reasons supporting both small and large sizes of the moving windows $\mathcal{W}(t)$. In the subsequent section we introduce three different window selection strategies, all of which make a compromise between small and large windows.

3.3 Selection of moving windows

Since the aim of the model evaluation is to compare several models against the same time series of realizations $y_{1},\ldots,y_{N}$, it is meaningful to determine the moving windows $\mathcal{W}(t)\subset\{1,\ldots,N\}$, $t=1,\ldots,N$ once based on the realized time series. Below we introduce three alternative window selection approaches based on the changepoint analysis described in Section 2.2. In a first step, $m<N-1$ changepoints $\bm{\tau}_{1:m}=(\tau_{1},\ldots,\tau_{m})$ of $y_{1},\ldots,y_{N}$ are detected using the PELT algorithm with the cost function in Equation (6) which assumes varying means and variances for different segments of the time series. Hence, it allows to detect changes in mean and/or variance. For the penalty constant $\kappa$ we choose $\kappa=p\ln(N)$ with $p=3$, since we count three additional parameters (mean, variance, changepoint) per segment. Moreover, we demand that the minimum segment length is $11$ (see our discussion in Section 3.2), that is it must hold $(\tau_{j+1}-\tau_{j})>10$, $j=0,\ldots,m$. In a second step, (different types of) moving windows are specified, based on the selected $(m+1)$-segmentation $\bm{\tau}_{0:(m+1)}=(0,\bm{\tau}_{1:m},N)$:

Overlapping windows with fixed width (OF)

For the first approach, we consider moving time windows $\mathcal{W}^{\operatorname{OF}}(t)$ of a fixed width $\omega^{\operatorname{OF}}\coloneqq\left|\mathcal{W}^{\operatorname{OF}}(t)\right|$ except for edge effects at both ends of the time interval $\{1,\ldots,N\}$. We center the overlapping (symmetric) windows $\mathcal{W}^{\operatorname{OF}}(t)$ around their window locations $t$. To obtain $\mathcal{W}^{\operatorname{OF}}(t)$, we first calculate the median segment length $\lambda\coloneqq\operatornamewithlimits{median}_{j=0,\ldots,m}\left(\tau_{j+1}-\tau_{j}\right)$, where $\operatornamewithlimits{median}_{\mathcal{I}}\left(\cdot\right)$ denotes the sample median of quantities indexed by the index set $\mathcal{I}$. Then, we compute the fixed window width parameter defined as $\delta^{\operatorname{OF}}\coloneqq\left\lfloor(\lambda-1)/2\right\rfloor$, where $\left\lfloor\cdot\right\rfloor$ rounds a number to its next smaller integer. Defining the window width parameter $\delta_{t}^{\operatorname{OF}}$ as

[TABLE]

we ensure that the (moving) overlapping windows with fixed width (OF-windows) defined as

[TABLE]

are symmetric and $\mathcal{W}^{\operatorname{OF}}(t)\subset\{1,\ldots,N\}$ for all $t=1,\ldots,N$. The corresponding window widths equal $\omega_{t}^{\operatorname{OF}}=2\delta_{t}^{\operatorname{OF}}+1$. The windows $\mathcal{W}^{\operatorname{OF}}(t)$ for $t=1,\ldots,\delta^{\operatorname{OF}}$ and $t=N-\delta^{\operatorname{OF}}+1,\ldots,N$ are the border cases. Moving scores obtained based on OF-windows will be referred to as OF-scores.

The top panel of Figure 1 illustrates two time series ($y_{t}$ and $x_{t}$, $t=1,\ldots,N$) of realizations from a phenomenon $Y_{t}$ and a corresponding model/prediction $X_{t}$. Moreover, it shows which changepoints $\bm{\tau}_{1:m}$ were detected by the PELT algorithm. The middle panel shows which window width was selected by the OF-method for each time instance $t=1,\ldots,N$. The corresponding moving windows $\mathcal{W}^{\operatorname{OF}}(t_{1})$ and $\mathcal{W}^{\operatorname{OF}}(t_{2})$ for two selected time instances $t_{1}$ and $t_{2}$ are illustrated as gray boxes in the upper panel. The lower panel shows for both OF-windows (located at $t_{1}$ and $t_{2}$), how the empirical CDFs for both the phenomenon ($\widehat{G}_{t}$) and the model/prediction ($\widehat{F}_{t}$) differ from each other and how much they deviate from the observation ($y_{t}$) in the window location ($t$).

Overlapping windows with varying width (OV)

For the second approach, we consider moving time windows $\mathcal{W}^{\operatorname{OV}}(t)$ with varying width $\omega_{t}^{\operatorname{OV}}\coloneqq\left|\mathcal{W}^{\operatorname{OV}}(t)\right|$. Again, we center the overlapping windows around their window location $t$. To obtain the varying windows $\mathcal{W}^{\operatorname{OV}}(t)$, we first determine the centers of the segments $\{\tau_{j}+1,\ldots,\tau_{j+1}\}$, $j=0,\ldots,m$, of the $(m+1)$-segmentation $\bm{\tau}_{0:(m+1)}$ defined as $\gamma_{j}\coloneqq(\tau_{j}+1+\tau_{j+1})/2$, $j=0,\ldots,m$. Then, we linearly interpolate the corresponding segment lengths $\lambda_{j}\coloneqq\left(\tau_{j+1}-\tau_{j}\right)$ between the segment centers $\gamma_{j}$, $j=0,\ldots,m$. Hence, the interpolated segment lengths are given by $\varsigma(t)\coloneqq[(\gamma_{j+1}-t)\lambda_{j}+(t-\gamma_{j})\lambda_{j+1}]/(\gamma_{j+1}-\gamma_{j})$, for $t\in[\gamma_{j},\gamma_{j+1}]$, $j=0,\ldots,m$. Then, the (symmetric) overlapping windows with varying width (OV-windows) are defined by

[TABLE]

where

[TABLE]

Here, $\left\lfloor\cdot\right\rfloor$ and $\left\lceil\cdot\right\rceil$ round a number to its closest smaller and larger integer, respectively. Note that for the border cases ($t=1,\ldots,\left\lfloor\gamma_{0}\right\rfloor$ and $t=\left\lceil\gamma_{m+1}\right\rceil,\ldots,N$) the moving windows are again defined such that they completely lie within $\{1,\ldots,N\}$. The varying window widths equal $\omega_{t}^{\operatorname{OV}}=2\delta_{t}^{\operatorname{OV}}+1$. Moving scores obtained based on OV-windows will be referred to as OV-scores. An illustration of the OV-method in analogy to Figure 1 is provided by Figure S.2 in the supplementary material.

Disjoint windows with varying width (DV)

For the third approach, we consider disjoint moving time windows $\mathcal{W}^{\operatorname{DV}}(t)$ with varying width $\omega_{t}^{\operatorname{DV}}\coloneqq\left|\mathcal{W}^{\operatorname{DV}}(t)\right|$ given by the segments of the $(m+1)$-segmentation $\bm{\tau}_{0:(m+1)}$. Hence, the disjoint windows with varying width (DV-windows) are defined by

[TABLE]

The varying window widths equal the corresponding segment lengths, $\omega_{t}^{\operatorname{DV}}\coloneqq\left(\tau_{j+1}-\tau_{j}\right)$, for $t=\tau_{j}+1,\ldots,\tau_{j+1}$, $j=0,\ldots,m$. Moving scores based on DV-windows will be referred to as DV-scores. An illustration of the DV-method in analogy to Figure 1 is provided by Figure S.3 in the supplementary material.

4 Simulation study

4.1 Organization of simulation study

To judge the applicability of moving scores based on different moving window selection approaches, we conduct a simulation study under three different scenarios. We first consider a changepoint scenario (C, Section 4.2), where the characteristics (mean and variance) of the phenomenon of interest change at certain (unknown) time instances. Moreover, with climate model evaluation in mind, we consider a trend scenario (T, Section 4.3) and a periodicity scenario (P, Section 4.4).

General setup

For each scenario, we consider a phenomenon $Y_{t}$ and five models $X_{t}^{k}$ for that phenomenon ($t=1,\ldots,N$), numbered by $k=1,\ldots,5$. We assume that $Y_{t}$, $X_{t}^{k}$, $k=1,\ldots,5$, are normally distributed for all $t=1,\ldots,N$, that is

[TABLE]

Since we simulate the phenomenon, we also refer to $Y_{t}$, $t=1,\ldots,N$, as the data generating process (DGP). Each result of the simulation study is based on $10{,}000$ replications of the data generating process (11) and the corresponding model (12).

Methods

For each replication, we compute time series of moving scores (7) for the different window selection approaches dicussed in Section 3.3. Moreover, we also compute point-wise (PW) scores. In case of the changepoint (C) and the trend scenario (T), we also consider ST-scores which assume stationarity of the time series. To get an idea of the temporal evolution of these (empirical) scores, we average the obtained score time series over all $10{,}000$ replications. Furthermore, we compare them to the corresponding (averaged) time series of theoretical (Theo.) scores, which can be obtained, since we know the distributions of the data generating process (11) and each of the five models (12).

4.2 Changepoint scenario

In order to describe the changepoint scenario (C) we consider time varying means and standard deviations which change after certain changepoints $\bm{\tau}_{1:m}\in\mathbb{R}^{m}$ (see Section 2.2). Hence, the time varying means and standard deviations are defined as $M_{\text{C}}(t;\bm{\tau}_{0:(m+1)},\bm{\mu})\coloneqq\mu_{j}$ and $S_{\text{C}}(t;\bm{\tau}_{0:(m+1)},\bm{\sigma})\coloneqq\sigma_{j}$, for $t\in\{\tau_{j}+1,\ldots,\tau_{j+1}\}$ and parameter vectors $\bm{\mu}=(\mu_{0},\ldots,\mu_{m})\in\mathbb{R}^{m+1}$ and $\bm{\sigma}=(\sigma_{0},\ldots,\sigma_{m})\in\mathbb{R}^{m+1}$, respectively.

Changepoint scenario (C)

The subsequent equations define the data generating process (C0) and the five models (C1)–(C5), where we consider time series of length $N=200$ with $m=2$ changepoints and corresponding 3-segmentation $\bm{\tau}_{0:3}=(0,80,130,200)$.

[TABLE]

[TABLE]

Whereas the true model (C1) is equivalent to the DGP (C0) in terms of its parametrization, models (C2)–(C5) differ from (C0) in at least one parameter. Model (C2) assumes a constant mean, but captures the change in the variance. While model (C3) models the change in the mean correctly, it assumes a constant variance. Models (C4) and (C5) are misspecified both in terms of mean and variance with the misspecification of the mean being less severe for (C5). Figure S.4 in the supplementary material illustrates one replication of the DGP (C0) and the corresponding selection of moving windows. The PELT algorithm finds two changepoints in nearly all the simulated series with an estimated fixed window width of $69$ or $71$ time points for most cases.

Time series of (moving) scores

Figure 2 compares models (C1)–(C5) based on time series of theoretical/empirical CRPS. We observe that the temporal evolution of the PW- and the ST-scores differs considerably from their theoretical counterparts. The observed relative model rankings are wrong for all time instances. We further observe that the moving OF- and OV-scores provide faulty model rankings for time instances close to the changepoints. The DV-scores mirror the temporal evolution of the theoretical scores best and mostly provide the correct model rankings.

Model rankings

Table 1 shows which evaluation approaches rank the models (C1)–(C5) correctly and which do not. It provides average SE and CRPS over time instances and replications for the different evaluation approaches (OF, OV, DV, PW and ST) and compares them to the corresponding theoretical (Th) values. The corresponding model rankings under each method are given in the right half of the table. If the score averages of two or more models are equal, all of these models get the same (minimum) rank. We expect the true model (C1) to be ranked lowest (best), and the ranking should be identical to that based on the theoretical (Th) scores. These two criteria are fulfilled for the CRPS calculated based on the OF, OV and DV approach. The SE score is unable to differentiate between all five models as its evaluation is based on the mean value only. In particular, the empirical evaluation approaches rank model (C3) higher than the true model (C1). In the theoretical SE in (1), the realization is evaluated against the expected value of the model/prediction. As model (C3) has the correct mean but a smaller variance than the true model (C1), the realized values under (C3) are likely to be closer to the true mean of the data generating process. The PW- and the ST-approach yield erroneous model rankings under both scores.

4.3 Trend scenario

For the trend scenario (T) we consider trends in the mean and the standard deviation and model these using the function $h_{\text{T}}(t;\bm{\theta})\coloneqq\theta_{0}+\theta_{1}t\exp\left(\theta_{2}t\right)$, where $\bm{\theta}\coloneqq(\theta_{0},\theta_{1},\theta_{2})\in\mathbb{R}^{3}$. Hence, for $t=1,\ldots,N$, we model the mean as $M_{\text{T}}(t;\bm{\mu})\coloneqq h_{\text{T}}(t;\bm{\mu})$, and the standard deviation as $S_{\text{T}}(t;\bm{\sigma})\coloneqq h_{\text{T}}(t;\bm{\sigma})$, where $\bm{\mu}=(\mu_{0},\mu_{1},\mu_{2})\in\mathbb{R}^{3}$ and $\bm{\sigma}=(\sigma_{0},\sigma_{1},\sigma_{2})\in\mathbb{R}^{3}$, respectively.

Trend scenario (T)

The subsequent equations define the data generating process (T0) and the five models (T1)–(T5), where we consider time series of length $N=200$.

[TABLE]

[TABLE]

Whereas the true model (T1) is equivalent to the data generating process (T0) in terms of its parametrization, models (T2)–(T5) differ from (T0) in at least one parameter. The mean of model (T2) is always too small while the variance is modeled correctly. Model (T3) models the exponentially increasing mean correctly with too small a variance. Models (T4) and (T5) are misspecified in both mean and variance, where the misspecification of the variance is less severe for model (T4). Figure S.5 in the supplementary material illustrates one replication of the DGP (T0) and the corresponding selection of moving windows. Figure S.7 compares models (T1)–(T5) based on time series of theoretical/empirical CRPS. Here, the PELT algorithm estimates between four and eight changepoints in the simulated series resulting in a fixed window width ranging from $15$ to $51$.

Model rankings

The resulting scores and the associated model rankings are shown in Table 2. As the theoretical SE evaluates the mean of the model/prediction only, the competing models are here ranked in only two groups. Similarly as for the changepoint scenario above, the misspecified model with correct mean and too small a variance is ranked higher than the true model using the emprical moving scores and the same effect is observed for the CRPS. Again we observe erroneous model rankings under the PW and the ST approaches.

4.4 Periodicity scenario

For the periodicity scenario (P) we consider time series with periodically varying mean and standard deviation. To model the periodicity, we consider the function $h_{\text{P}}(t;\bm{\theta})\coloneqq\theta_{0}+\theta_{1}\sin\left(2\pi t\theta_{2}\right)$, where $\bm{\theta}\coloneqq(\theta_{0},\theta_{1},\theta_{2})\in\mathbb{R}^{3}$. Then, for $t=1,\ldots,N$, we model the mean as $M_{\text{P}}(t;\bm{\mu})\coloneqq h_{\text{P}}(t;\bm{\mu})$, and the standard deviation as $S_{\text{P}}(t;\bm{\sigma})\coloneqq\exp\left(h_{\text{P}}(t;\bm{\sigma})\right)$, with parameter vectors $\bm{\mu}=(\mu_{0},\mu_{1},\mu_{2})\in\mathbb{R}^{3}$ and $\bm{\sigma}=(\sigma_{0},\sigma_{1},\sigma_{2})\in\mathbb{R}^{3}$, respectively.

Periodicity scenario (P)

The subsequent equations define the data generating process (P0) and the five models (P1)–(P5), where we consider time series of length $N=730$.

[TABLE]

[TABLE]

Whereas the true model (P1) is equivalent to the data generating process (P0) in terms of its parametrization, models (P2)–(P5) differ from (P0) in at least one parameter. Model (P2) underestimates the magnitude of the mean oscillations, whereas it captures the variance oscillations correctly. Model (P3) captures the mean oscillations correctly, however, the oscillations in the variance are underestimated. Models (P4) and (P5) are misspecified both in terms of mean and variance. Model (P4) models some of the variation in the variance while Model (P5) assumes a constant variance. Figure S.6 in the supplementary material illustrates one replication of the DGP (P0) and the corresponding selection of moving windows. Figure S.8 compares models (P1)–(P5) based on time series of theoretical/empirical CRPS. For this scenario, the PELT algorithm estimates between $24$ and $34$ changepoints resulting in the fixed window width ranging from $13$ to $23$.

Model rankings

Table 3 shows scores and model rankings for models (P1)–(P5) under the different evaluation approaches. Here, the empirical rankings for the CRPS calculated based on the OF, OV and DV approach are identical to the corresponding theoretical rankings. As before, the SE is unable to differentiate between models with identical mean value structure and the empirical approaches fail to recognize the true model compared to a model with correct mean structure and too small spread. Again, the PW approach yields wrong model rankings.

5 Case study: Evaluation of Regional Climate Models

In an application of the presented methodology, we evaluate four selected Regional Climate Models (RCMs) from the ENSEMBLES project (van der Linden and Mitchell,, 2009). The purpose of RCMs is to dynamically downscale gridded climate data from global models with a coarse spatial resolution (say $\sim 200$–$300$ km) to obtain climate simulations of a finer resolution ($\sim 25$–$50$ km) for a specific region. As a part of the ENSEMBLES project, ensembles of RCM simulations from $16$ different (European) climate research institutes were compiled. In an initial $40$-year experiment ($1961$–$2000$) covering Europe, the RCMs were driven by the ERA-40 reanalysis data set. This experiment was used to evaluate the RCMs (see also Section 1). Furthermore, RCM experiments driven by Global Circulation Model (GCM) output were conducted to create an ensemble of regional climate change projections for Europe. Most of these RCM simulations cover the period $1951$–$2050$, and are available on a daily temporal resolution. For our case study we consider only the output variable 2 meter temperature (daily mean temperature). For the control period $1961$–$1990$ and a spatial resolution of $25$ kilometers, we compare four RCMs:

• •

DMI-HIRHAM of the Danish Meteorological Institute (DMI),

• •

KNMI-RACMO2 of the Royal Netherlands Meteorological Institute (KNMI),

• •

MPI-M-REMO of the Max-Planck-Institute for Meteorology (MPI), and

• •

SMHIRCA of the Swedish Meteorological and Hydrological Institute (SMHI).

In the following, we refer to the different models by the acronym of the corresponding institute. We evaluate the models for both ERA-40 and GCM boundary conditions. All four RCM simulations we consider were driven by the ECHAM5 GCM. To differentiate between the two boundary conditions, we use the acronyms ERA-40 and ECHAM5.

The RCM model output are evaluated against the E-OBS observational data product (version 13.1, Haylock et al.,, 2008). This data product is a gridded data set interpolated from station observations. It covers the European continent (land only) and is available for different grids and spatial resolutions. For our evaluation we use the gridded daily mean temperature observations provided in version 13.1 of the data set. Our evaluation is based on all days in the period $1961$-$1990$ at a total of $10937$ grid points.

5.1 Evaluation approaches

To evaluate and compare the different RCMs, we calculate moving CRPS based on OF-windows (see Equation (8)), OV-windows (see Equation (9)), and DV-windows (see Equation (10)). For comparison, we also compute the point-wise (PW) scores for comparison. Recall that the point-wise CRPS is equivalent to the absolute error.

5.2 Overall model assessment

Table 4 provides overall (spatial and temporal) averages of moving CRPS and point-wise CRPS. As expected, the ERA-40 driven models achieve better average scores than the GCM driven models. While the CRPS calculations yield identical model rankings for all three types of moving windows (OF, OV, DV), the point-wise scores yield a different ranking. As we learned in our simulation study in Section 4, point-wise scores might lead to a spurious model ranking, since they do not account for higher order structures of the observed/modeled phenomenon. The results of Table 4 corroborate these findings. Moving scores provide a more holistic model assessment and hence should be preferred over a point-wise evaluation. It seems, that for the application at hand the type of the moving window is not crucial.

5.3 Temporal evaluation

To assess the model accuracy over time we consider spatial average scores over the entire study area, see Figure 3. For illustration purposes, we further smooth the averaged daily series by computing monthly averages. Figure 3 compares the resulting time series for the moving CRPS with OV-windows, and for the point-wise CRPS. The corresponding figures for OF- and DV-windows can be found in the supplementary material (Figure S.9).

All four models show pronounced seasonal oscillations in the scores, indicating a better performance in the summer and a worse performance in the winter. The osciallations are less pronounced for the ERA-40 driven models than for the ECHAM5 driven models. From that we conclude that the driving GCM (here ECHAM5) is not able to accurately capture the (regional) seasonal temperature dynamics. While for the moving CRPS the differences between the different models and the different boundary conditions (ECHAM5 and ERA-40) appear to be rather small, the point-wise evaluation indicates larger differences. This is due to the fact that the ECHAM5 boundary conditions are simulations from a GCM and, as such, are not designed to capture the weather in a particular month or a season, while the ERA-40 boundary conditions are a reanalysis based on observations which should capture the changes in weather.

Focusing on the moving CRPS, we observe time periods of weaker (e.g. $1985$–$1987$) and of better (e.g. $1974$–$1977$) model performance. Moreover, we find that for the DMI model the discrepancy between the two different boundary conditions is more obvious compared to the other models. While the ERA-40 driven DMI model apparently performs best in comparison to the other models, the ECHAM5 driven DMI model performs worst (see also Table 4).

5.4 Spatial evaluation

To judge the model performance for different locations across Europe we look at the temporal averages of the scores. The temporal averages corresponding to the moving CRPS based on OV-windows and the point-wise CRPS are depicted in Figure 4. The corresponding figures for OF- and DV-windows are provided in the supplementary material (Figure S.10).

Again the distinction between the ERA-40 and the ECHAM5 boundary conditions is more pronounced for the point-wise score which does not account for higher order structures in the temperature series. In particular, a point-wise evaluation does not account for the associated uncertainty which may vary across space. Among the ECHAM5 driven models, the DMI model performs worst, especially regarding the continental climate in eastern Europe. The MPI model under ECHAM5 performs badly in north-eastern Europe where lakes dominate the landscape. Figure 4 indicates that all the models have difficulties in modeling the temperature over large water bodies (e.g. lakes Ladoga and Onega in Russia). This is identified to a lesser degree by the point-wise score. Furthermore, we observe that most models have problems for mountainous areas (e.g. Alps, Pyrenees, Carpathians). To show for each pixel of the study area which model has the lowest (best) average score, we provide Figure S.11 in the supplementary material.

5.5 Evaluation based on linear trends

For comparison, we now consider one of the evaluation approaches undertaken in the literature. Following the approach of Lorenz and Jacob, (2010), we compare linear temperature trends in the RCM output and the E-OBS reference data.

To quantify the trends, we first aggregate the daily time series to yearly time series $y_{t}$, $t=1961,\ldots,1990$, by computing annual means. We then assume a linear regression model $y_{t}=\alpha+\beta t+{\varepsilon}_{t}$, with intercept $\alpha$, linear trend parameter $\beta$ and residuals ${\varepsilon}_{t}$ ($t=1961,\ldots,1990$). The parameter estimates are obtained using least-squares estimation separately for all $n=10937$ pixels/locations $i=1,\ldots,n$ of the study area and each data set. Subsequently, we distinguish between estimates $\widehat{\beta}_{i}^{\text{mod}}$ corresponding to one of the RCM outputs and estimates $\widehat{\beta}_{i}^{\text{ref}}$ corresponding to the reference data set. To compare and rank the different RCMs, we finally compute absolute trend errors of the form $\left|\widehat{\beta}_{i}^{\text{mod}}-\widehat{\beta}_{i}^{\text{ref}}\right|$, $i=1,\ldots,n$.

To summarize the results of the linear trend-based evaluation, we consider spatial averages $1/n\sum_{i=1}^{n}\left|\widehat{\beta}_{i}^{\text{mod}}-\widehat{\beta}_{i}^{\text{ref}}\right|$ of the absolute trend errors, see Table 5. For the ECHAM5 boundary conditions, the DMI model is considered best. For ERA-40 boundary conditions, the models are ranked in the order SMHI (best), KNMI, DMI and MPI (worst). These results differ considerably from the score-based evaluation in Table 4. Here, direct comparison of the DMI model for the two different boundary conditions further suggests that the model driven by ECHAM5 performs better. To see for each pixel of the study area which of the four RCMs has the smallest absolute trend error, we provide Figure S.12 in the supplementary material.

6 Conclusions and outlook

We propose methodology for the evaluation of time series models/predictions under the presence of non-stationarity. Our approach utilizes proper scores on moving time windows where stationarity is assumed, in order to provide a fair and holistic comparison of models/predictions. A simulation study explores the moving score technique under the presence of changepoints, trends and periodicity. A case study, comparing model output from regional climate models, illustrates the utility of the technique for practical applications. We summarize our results and conclusions in the following.

Being based on proper scores, our approach allows for a fair comparative model/prediction evaluation. Evaluation based on moving windows instead of a point-wise comparison accounts for higher order structures of the modeled/predicted phenomenon. Compared to most other evaluation approaches used in practice our evaluation method does not require a separate consideration of different seasons when considering sub-seasonal evaluation, or to make a (subjective) decision on the importance/weight of different features of the model/prediction. To base the choice of the moving windows on (detected) changepoints is meaningful, since changepoints divide a time series into stationary segments.

Note that the sequential comparison of individual observations within a time window against an empirical distribution of the model output from the same time window using a proper scoring rule is equivalent to a direct comparison of the model distribution and the corresponding observation distribution using an associated score divergence (Thorarinsdottir et al.,, 2013). In a multi-model comparison, the proper scoring rule and its score divergence will yield the same model rankings. However, the application of the divergence function in our setting requires some care, as overlapping moving windows and windows of varying sizes may result in aggregated values that include unequal weighting of the observations. Such weighting may destroy the propriety of the metric (Gneiting and Ranjan,, 2011).

Depending on the phenomenon of interest, certain types of moving windows might be more adequate than others. For instance, disjoint windows with varying width (DV) for fixed changepoints, overlapping windows with fixed (OF) or varying width (OV) under the presence of trends or periodicity, etc. Similar as in the simulated examples, our case study showed identical overall model rankings based on moving CRPS for different window selection approaches. For the application of the moving score methodology to continuous outcomes utilization of the CRPS is preferred over the SE score as it evaluates the full predictive distribution rather than the predictive mean only. In general, the simulation study showed that moving scores are able to approximate their theoretical counterparts considerably well, are better suited than point-wise (PW) scores or scores under a stationarity assumption (ST) to yield an adequate (true) model/prediction ranking and yield the correct model ranking in most cases (see the trend scenario (Section 4.3) for an exception).

An extension of the work presented here to categorical outcomes is of interest. For example, one might want to judge models which differentiate only if a certain (extreme) event occurs or not. Further extensions of our approach might concern the judgment of model accuracy in a spatial context, where we can think of the moving window as a spatial neighborhood considering either a fixed number of spatial neighbors or all locations within a certain radius.

Acknowledgements

Tobias Erhardt is supported by Deutsche Forschungsgemeinschaft (DFG) through the TUM International Graduate School of Science and Engineering (IGSSE). Claudia Czado is supported by the DFG (grant CZ86/4-1). Thordis L. Thorarinsdottir is supported by the Research Council of Norway through grant number 243953 “Physical and Statistical Analysis of Climate Extremes in Large Datasets” (ClimateXL). The ENSEMBLES and the E-OBS data used in this work were funded by the EU FP6 Integrated Project ENSEMBLES (Contract number 505539, http://ensembles-eu.metoffice.com) which we gratefully acknowledge. Moreover, we acknowledge the data providers of the E-OBS dataset in the ECA&D project (http://www.ecad.eu). The numerical computations were performed on a Linux cluster supported by DFG grant INST 95/919-1 FUGG.