Probabilistic Temperature Forecasting with a Heteroscedastic Autoregressive Ensemble Postprocessing model

TL;DR

This paper introduces a heteroscedastic autoregressive ensemble postprocessing model for temperature forecasts, improving reliability and extending forecast horizons by accounting for forecast error variability and ensemble spread.

Contribution

The paper presents a novel heteroscedastic extension of an autoregressive postprocessing model, incorporating high-resolution forecasts and enabling better long-horizon temperature predictions.

Findings

Significant improvements for forecasts beyond 24 hours.

Enhanced model reliability through variance dependence on ensemble spread.

Implementation available in the R package ensAR.

Abstract

Weather prediction today is performed with numerical weather prediction (NWP) models. These are deterministic simulation models describing the dynamics of the atmosphere, and evolving the current conditions forward in time to obtain a prediction for future atmospheric states. To account for uncertainty in NWP models it has become common practice to employ ensembles of NWP forecasts. However, NWP ensembles often exhibit forecast biases and dispersion errors, thus require statistical postprocessing to improve reliability of the ensemble forecasts. This work proposes an extension of a recently developed postprocessing model utilizing autoregressive information present in the forecast error of the raw ensemble members. The original approach is modified to let the variance parameter depend on the ensemble spread, yielding a two-fold heteroscedastic model. Furthermore, an additional high-resolution forecast is included into the postprocessing model, yielding improved predictive performance. Finally, it is outlined how the autoregressive model can be utilized to postprocess ensemble forecasts with higher forecast horizons, without the necessity of making fundamental changes to the original model. We accompany the new methodology by an implementation within the R package ensAR to make our method available for other researchers working in this area. To illustrate the performance of the heteroscedastic extension of the autoregressive model, and its use for higher forecast horizons we present a case study for a data set containing 12 years of temperature forecasts and observations over Germany. The case study indicates that the autoregressive model yields particularly strong improvements for forecast horizons beyond 24 hours.

Figures

Tables

Table 1. Table 1: Verification metrics of EMOS and (heteroscedastic) AR-EMOS for station Magdeburg aggregated over 4341 verification dates

	CRPS	DSS	RMV	Var(PIT)
EMOS	0.8415	2.0918	1.3670	0.0946
AR-EMOS	0.8309	1.9149	1.3825	0.0876

Table 2. Table 2: Verification metrics of EMOS and (heteroscedastic) AR-EMOS aggregated over 4341 × 76 4341 76 4341\times 76 verification cases (dates and stations)

	CRPS	DSS	RMV	Var(PIT)
EMOS	0.9057	2.1523	1.4907	0.0954
AR-EMOS	0.9033	2.0129	1.5322	0.0883

Table 3. Table 3: Verification metrics of EMOS ∗ , (heteroscedastic) AR-EMOS ∗ , and SLP combination SLP ∗ of both, aggregated over 4251 verification dates at station Magdeburg

	CRPS	DSS	RMV	Var(PIT)
EMOS∗	0.8223	2.0304	1.3714	0.0908
AR-EMOS∗	0.8097	1.8404	1.3663	0.0849
SLP∗	0.8000	1.9043	1.3965	0.0854

Table 4. Table 4: Verification metrics of EMOS ∗ , (heteroscedastic) AR-EMOS ∗ , and SLP combination SLP ∗ of both, aggregated over 4251 × 76 4251 76 4251\times 76 verification cases (dates and stations)

	CRPS	DSS	RMV	Var(PIT)
EMOS∗	0.8712	2.1412	1.4250	0.0931
AR-EMOS∗	0.8685	1.9270	1.4950	0.0854
SLP$\ast$	0.8460	1.9350	1.5031	0.0860

Table 5. Table 5: Verification metrics of EMOS(ENS) and (heteroscedastic) AR-EMOS(ENS) aggregated over 4340 × 76 4340 76 4340\times 76 (48-h) and 4339 × 76 4339 76 4339\times 76 (72-h) verification cases (dates and stations)

	CRPS	DSS	RMV	Var(PIT)
EMOS 48-h	1.0101	2.4147	1.6156	0.0979
AR-EMOS 48-h	0.9897	2.1749	1.7263	0.0872
EMOS 72-h	1.1244	2.6353	1.7831	0.099
AR-EMOS 72-h	1.0949	2.3548	1.9591	0.086

Equations

$$Z_m(t):=Y(t)-x_m(t)$$

$$Z_m(t)-{\alpha }_m=\sum _{j=1}^{p_m}{\beta }_{m,j}[Z_m(t-j)-{\alpha }_m]+{\epsilon }_m(t),$$

$$Y(t)={\tilde{x}}_m(t)+{\epsilon }_m(t),$$

$${\tilde{x}}_m(t)=x_m(t)+{\alpha }_m+\sum _{j=1}^{p_m}{\beta }_{m,j}[y(t-j)-x_m(t-j)-{\alpha }_m]$$

$${\gamma }_m^2(t)$$

$$Y(t)|x_1(t),\dots ,x_M(t)\sim \mathcal{N}(\mu (t),{\sigma }^2(t)),$$

$$\mu (t)=\frac{1}{M}\sum _{m=1}^M{\tilde{x}}_m(t).$$

$$\sigma (t)=w\sqrt{\frac{1}{M}\sum _{m=1}^M{\gamma }_m^2(t)}+(1-w)\sqrt{{\tilde{S}}^2(t)},$$

$$\mu (t)$$

$$\sigma (t)$$

$$F(x)=w_1{G}_1(x)+w_2{G}_2(x),G_l(x)=\Phi \left(\frac{x-{\mu }_l}{{\sigma }_{l}c}\right),$$

$$\text{CRPS}(F,y_{\text{obs}})={\int }_{-\infty }^{\infty }{\left\{F(y)-1(y\ge y_{\text{obs}})\right\}}^2\text{d}y,$$

$$\text{DSS}(F,y_{\text{obs}})=\frac{(y_{\text{obs}}-{\mu }_F{)}^2}{{\sigma }_F^2}+2\log {\sigma }_F,$$

$$S=\sqrt{T}\frac{\overline{d}}{\sqrt{\sum _{\tau =-(h-1)}^{h-1}{\hat{\gamma }}_d(\tau )}},$$

$$\overline{d}=\frac{1}{T}\sum _{t=1}^{T}d(t),d(t)=s_1(t)-s_2(t),$$

$${\hat{\gamma }}_d(\tau )=\frac{1}{T}\sum _{t=|\tau |+1}^T(d(t)-\overline{d})(d(t-|\tau |)-\overline{d})$$

$$H_1:CRP{S}_{\text{AR-EMOS}}<CRP{S}_{\text{EMOS}}.$$

Full text

Keywords: Ensemble postprocessing; predictive probability distribution; autoregressive process; spread-adjusted linear pool; heteroscedastic model, spread-error correlation, high-resolution forecast

1 Introduction

Today, weather prediction is based on numerical weather prediction (NWP) models. Such models are deterministic in nature and represent the dynamical physics of the atmosphere by a set of differential equations. The current state of the atmosphere is evolved forward in time to predict future atmospheric states. The solutions strongly depend on the initial conditions and model formulations. Thus, NWP models suffer from several sources of uncertainties. Common practice in addressing these uncertainties is the use of ensemble prediction systems (EPS). The NWP model is run multiple times, each time with variations in the model parameterizations and/or initial and boundary conditions (Gneiting and Raftery, 2005, Leutbecher and Palmer, 2008). A forecast ensemble can be viewed as probabilistic forecast that allows to assess forecast uncertainty (Palmer, 2002). However, in practise, NWP ensembles exhibit forecast biases and dispersion errors, and require statistical postprocessing to improve calibration and forecast skill by utilizing recently observed forecast errors and observations. An additional benefit is that statistical postprocessing can yield full predictive probability distributions (Gneiting and Raftery, 2005, Gneiting and Katzfuss, 2014, Wilks and Hamill, 2007).

Statistical postprocessing models have enjoyed increasing popularity and success during the last decades and a variety of models tailored to specific problems have been developed. Many of the recently proposed models are extensions and modifications of two generic state-of-the-art models, namely the Ensemble Model Output Statistics approach (EMOS; Gneiting et al., 2005) and the Bayesian Model Averaging (BMA; Raftery et al., 2005). The original EMOS and BMA models were designed for Gaussian distributed weather quantities, and a variety of modifications for other weather quantities have been developed (see, e.g., Hemri et al., 2014, Gneiting and Katzfuss, 2014, Schefzik et al., 2013, for an overview).

Further, postprocessing models allowing to incorporate inter-variable, spatial, or temporal dependence structures have gained increased interest (see, e.g., Möller et al., 2013, Baran and Möller, 2015, Schefzik et al., 2013, Bouallegue et al., 2015, Berrocal et al., 2007, Feldmann et al., 2015, Kleiber et al., 2011, Wilks, 2015, Vrac and Friederichs, 2015, Schuhen et al., 2012, Pinson, 2012, Möller et al., 2016, Hemri et al., 2015).

A general overview on various aspects of ensemble postprocessing can be found in Vannitsem et al., (2018).

In line with the need for models incorporating dependencies explicitly, Möller and Groß (2016) introduced a postprocessing model for Gaussian distributed weather quantities (such as temperature) that accounts for dependencies of the individual ensemble forecasts across time. In this regard, the model utilizes the autoregressive information present in the forecast error of the individual raw ensemble forecasts to set up corrected ensemble forecasts as well as a predictive distribution.

The work presented here extends the AR-EMOS model of Möller and Groß (2016) to be of heteroscedastic (or “non-homogeneous”) nature, meaning that the variance parameter of the model varies with the (empirical) ensemble spread. A postprocessing model using a heteroscedastic variance parameter accounts for the well-known spread-error correlation (Barker, 1991, Whitaker and Loughe, 1998) of forecast ensembles, stating that there is a positive association between the forecast error (or predictive skill) and the spread of the ensemble. Thus, the extended AR-EMOS model incorporates heteroscedasticity in “two directions”, namely across time (longitudinal) for each individual member, and across the ensemble members (cross-sectional), to account for above mentioned spread-error-correlation. Therefore, the approach allows for features not possible with standard postprocessing models, such as fitting a predictive distribution based only on a single ensemble member. This feature is investigated on the basis of an additional high-resolution forecast added to the ensemble, which is known to improve predictive performance to a great extend. While in the original paper of Möller and Groß (2016) the AR-EMOS model was only applied to 24-h ahead ensemble forecasts, in general it is also applicable to (arbitrary) other forecast horizons. In this follow-up work we explain how the AR-EMOS model can be used for other than 24-h ahead forecast horizons and present results on predictive performance.

The development of postprocessing models accounting for specific problems is a quite active area of research, however, not all of the software carrying out model fitting for these recently proposed methods is publicly available. Prominent examples of postprocessing software implemented under the statistical software environment R (R Core Team, 2018) are the state-of-the-art EMOS and BMA models (Yuen et al., 2018, Fraley et al., 2018), a recently developed heteroscedastic logistic model for ensemble postprocessing (Messner et al., 2013, 2014, 2016b), implemented in the package crch (Messner et al., 2016a), and an implementation of verification metrics to assess probabilistic forecasts in the package scoringRules (Jordan et al., 2016).

In line with the need for publicly available postprocessing software, this follow-up work to the methodology presented in Möller and Groß (2016) is accompanied by an implementation within an R package called ensAR (Groß and Möller, 2019), which can currently be installed from the Git repository hosting service GitHub, e.g. by by using the package devtools (Wickham and Chang, 2016).

A case study for temperature forecasts of the European Center for Medium Range Weather Forecasts (ECMWF; Buizza et al., 2007) over Germany is carried out to illustrate the performance and properties of the proposed heteroscedastic autoregressive postprocessing model.

2 Methods

2.1 Individual Ensemble Member Postprocessing

Suppose that (ensemble) forecasts are initialized at a fixed time point and predict a weather quantity a fixed time step ahead (forecast horizon), which is (for now) not greater than 24 hours. If $t$ denotes the time point (day and hour) for which the forecast is valid, the data consists of forecasts $x_{1}(t),\ldots,x_{M}(t)$, and a matching weather quantity observation $y(t)$ for $t=1,\ldots,T$. E.g. if a forecast is initialized at 12 UTC and predicts $18$ hours ahead, the forecast is valid at 06 UTC of the following day. Given the initialization time of the forecasts is fixed (which would usually be the case) the data (observation and ensemble members) is a collection of evenly (24 hour) spaced time series referring to the respective validation time point. Let

[TABLE]

be the time series of forecast errors of the individual ensemble members $x_{m}(t)$. In case $x_{m}(t)$ would be obtained from an appropriate model for $Y(t)$, the error series $Z_{m}(t)$ can be assumed to be white noise. However, Möller and Groß (2016) found that the observed individual error series $z_{m}(t)$ exhibit substantial autoregressive behaviour. The authors propose to utilize this residual autoregressive information to obtain a corrected (AR-adjusted) forecast ensemble and to define a predictive distribution based on this AR-adjusted ensemble forecasts.

In this regard, it is assumed that each $\{Z_{m}(t)\}$ follows an autoregressive process of order $p_{m}$, denoted by $\text{AR}(p_{m})$, i.e.

[TABLE]

where $\{\varepsilon_{m}(t)\}$ is white noise with expectation $\text{E}(\varepsilon_{m}(t))=0$ and variance

$\text{Var}(\varepsilon_{m}(t))=\sigma_{m}^{2}$. Then the random variable $Y(t)$ representing the weather quantity can be written as

[TABLE]

where

[TABLE]

can be viewed as a “corrected” forecast member for $y(t)$ based on the original ensemble member $x_{m}(t)$, $x_{m}(t-1)$, …, $x_{m}(t-p_{m})$ at past time points up to and including $t$, and the observation $y(t-1)$, …, $y(t-p_{m})$.

Performing the described procedure for each ensemble member $x_{m}(t)$ individually yields an “AR-adjusted” or “corrected” forecast ensemble $\widetilde{x}_{1},\ldots,\widetilde{x}_{M}$.

This approach of obtaining a corrected forecast ensemble rather than a predictive probability distribution/density has a connection to so-called “member-by-member-postprocssing” (MBMP), Van Schaeybroeck and Vannitsem, 2015, Schefzik, 2017, see, e.g.). MBMP approaches have gained increased interest, as they retain the dependence structure inherent in the original raw forecast ensemble, while this implicit dependence information is often lost when performing (univariate) postprocessing.

2.2 Forecast Error Variance

The variance of the autoregressive process $\{Z_{m}(t)\}$ in Equation (1),

$\text{Var}(Z_{m}(t))=:\gamma_{m}^{2}(t)$, is given as

[TABLE]

where $\rho_{m}(k)$ is the autocorrelation function of the process $\{Z_{m}(t)\}$ at lag $k$, see e.g. Cryer and Chan (2008, Eq. (4.3.31)). In R, the autocorrelation function of an ARMA process can be computed with the function ARMAacf.

2.3 Different Forecast Horizons

If the (ensemble) forecasts are less or equal to 24 hours ahead, it is obvious how to obtain the corrected ensemble forecasts $\widetilde{x}_{m}$ based on the AR-fit to the error series $Z_{m}(t)$, as the values $x_{m}(t)$, $x_{m}(t-1)$, …, $x_{m}(t-p_{m})$, as well as $y(t-1)$, …, $y(t-p_{m})$ are readily available at time point $t$. So the required parameters of the AR process can be directly estimated from the (observed) training series $z(t-s),\ldots,z(t-2),z(t-1)$ of length $s$.

However, if the considered forecast horizon is in the interval $(24\text{h},48\text{h}]$ and $t$ is any time point (day and hour) for which the forecast(s) are valid, then the observed error $z(t-1)$ is not available, as $y(t-1)$ has not yet been observed. Therefore, a preprocessing step is introduced before the AR-adjusted ensemble can be obtained. To compensate for the unavailable observed error $z(t-1)$, an AR process is fitted in advance to $z(t-s),\ldots,z(t-3),z(t-2)$, and $z(t-1)$ is predicted from the respective AR model fit.

After predicting $z(t-1)$, the complete error series $z(t-s),\ldots,z(t-2),z(t-1)$ is available again, and the AR-adjusted ensemble can be obtained as described in the previous section. For forecast horizons greater than 48-h, the observed errors at more time steps than $t-1$ are missing (e.g. for 72-h $z(t-1)$ and $z(t-2)$ are missing, and so on). The same procedure can then be applied to predict the missing errors from a fitted AR-model based on the past errors still available.

2.4 Heteroscedastic Autoregressive Predictive Distribution

Möller and Groß (2016) assume the predictive distribution for $Y(t)$ to be Gaussian, that is

[TABLE]

where $\mu(t)$ may be a function of the ensemble members and $\sigma^{2}(t)$ may be a function of the ensemble variance.

To account for the well-known spread-error correlation in ensemble forecasts, this follow-up work proposes an improved model for the predictive variance $\sigma^{2}(t)$ in a similar fashion as the variance term is defined in the EMOS model (Gneiting et al., 2005).

The predictive mean $\mu(t)$ is defined as the average over the AR-adjusted ensemble members

[TABLE]

However, an extended variance model is suggested, combining the estimated dispersion of the error process (longitudinal variation) retrieved solely from the past of each ensemble member $x_{m}$ with the spread of the corrected members $\widetilde{x}_{m}$ (cross-sectional variation). To be more specific, the model is given by

[TABLE]

where $\gamma_{m}^{2}(t)$ is the variance of the AR process as given in Equation (2), $\widetilde{S}^{2}(t)$ is the empirical variance of the AR-corrected ensemble, and $w\in[0,1]$ is a weight obtained by minimizing the CRPS over an (additional) training period. A simple (longitudinal) variance model can be obtained as a special case by setting $w=1$.

The classical EMOS model is also known by the term non-homogeneous regression, due to a variance model that is non-constant (non-homogeneous/heteroscedastic) with respect to the spread in the ensemble. Therefore, the modification of the AR-EMOS model with the variance model (5) incorporating the ensemble spread, is called heteroscedastic AR-EMOS - following the EMOS nomenclature. As the AR-EMOS method is based on a time-series model, and does not only consider the variation in the ensemble spread (cross-sectional part), but also the variation across time (longitudinal part), the more general (and in statistics more common) term heteroscedastic is used rather than the term non-homogeneous.

2.5 Postprocessing a Single Forecast

When there exists a single distinguished forecast $x_{\ast}(t)$, e.g. the high-resolution forecast $x_{\text{hres}}(t)$ described below, it is still possible to obtain a corresponding predictive distribution by the described AR-EMOS method.

The parameters $\mu_{\ast}(t)$ and $\sigma_{\ast}(t)$ are in principle estimated in the same way as those corresponding to the regular members. However, the mean (4) and the longitudinal part of the variance formula (5) reduce to a single summand, and the second part of the variance model in Equation (5) becomes zero. Nonetheless, the variance $\gamma_{\ast}^{2}$ of the error series corresponding to the individual member $x_{\ast}(t)$ can still be computed - on the basis of past values and the AR-fit. Thus, the original AR approach and its refined version presented here both allow to estimate the variance parameter and fit a predictive distribution based only on a single forecast. This is a new and useful feature not available in other state-of-the-art postprocessing models.

2.6 High-Resolution Forecast

In this section we describe how the high-resolution forecast known to improve predictive performance (see, e.g. Kann et al., 2009, Gneiting, 2014, Persson, 2015) can be included into the AR-EMOS postprocessing model. We call the AR-EMOS model including the additional forecast an extended model.

Let again $x_{1}(t),\ldots,x_{M}(t)$ denote the forecast ensemble. However, this time the $M$ members comprise of the regular (exchangeable) forecasts described in Section 2.1 and the additional high-resolution forecast $x_{\text{hres}}(t)$. That means the (total) number $M$ of forecasts utilized is actually increased by one. As described in Section 2.1 it is again assumed that the forecast errors $Z_{m}(t)$, $m=1,\ldots,M$ follow an $\text{AR}(p_{m})$ process, yielding the AR-corrected ensemble $\widetilde{x}_{1}(t),\ldots,\widetilde{x}_{M}(t)$, $\widetilde{x}_{\text{hres}}(t)$, which is the basis for estimating mean and variance of the predictive distribution.

As the high-resolution forecast has somewhat different properties than the regular ensemble members, an apparent approach may be to treat them as two different (exchangeable) groups with respect to the parameters of the predictive distribution. This course of action is quite common in ensemble post-processing models, ensemble members belonging to a certain (exchangeable) group are assumed to share the same coefficients in the model (Gneiting, 2014).

To account for the above mentioned groups, each parameter in the AR-EMOS model is defined as (equally weighted) sum of the respective group-wise parameters, i.e.,

[TABLE]

Here, $\mu_{\text{ens}}(t)$ is estimated as already stated in Equation (4). The parameters $\mu_{\text{hres}}(t)$ and $\sigma_{\text{hres}}(t)$ corresponding to the high-resolution forecast are estimated as described in the precious subsection.

Assigning each of the group-specific parameters fixed and equal weights is an initial and relatively straightforward approach for demonstrating the general idea. In the case study it will be shown that this simple version already yields good results for predictive performance.

Of course this rather simple method can be modified to be more data driven, that is using weights for the group-specific parameters directly estimated from data, e.g. by minimum CRPS estimation. Furthermore, the approach for two groups (regular ensemble members and high-resolution forecast) can be generalized to include multiple (exchangeable) groups, which also need not necessarily contribute equally to predictive performance, which makes it reasonable to estimate the weights for the group-specific parameters from data.

One possibility for a more general definition of an AR-EMOS group model can be accomplished for example by combining (group-wise) predictive distributions with the spread-adjusted linear pool already employed by Möller and Groß (2016) and shortly described in the following subsection.

2.7 Combination of Predictive Distributions

Möller and Groß (2016) proposed to combine the predictive distribution of classical EMOS and AR-EMOS in a spread-adjusted linear pool (SLP, Gneiting and Ranjan, 2013). For the special case of combining $n=2$ predictive distributions, the SLP combination has cumulative distribution function (CDF)

[TABLE]

$l=1,2$, where $w_{1}$ is a non-negative weight parameter, $w_{2}=1-w_{1}$, and $c$ is a strictly positive spread adjustment parameter. Here, $\phi$ and $\Phi$ denote the probability density function (PDF) and the cumulative distribution function (CDF) of the standard normal distribution, respectively. The two distributions $G_{1}$ and $G_{2}$ can be fitted separately by postprocessing models of choice, and the weights are obtained by minimizing a verification score (specifically the CRPS) over a training period, for fixed and given $G_{1},G_{2}$.

The original approach was proposed with the aim to improve predictive performance by combining two predictive distributions coming from different sources (see also Gneiting and Ranjan, 2013). This approach can directly be extended to a (finite) number $n>2$ of predictive distributions.

3 Tools to Asses Predictive Performance

3.1 Scoring Rules

A common tool to assess the quality of probabilistic forecasts are scoring rules. They assign a scalar to a pair ($y$,$F$), where $y$ is the verifying observation and $F$ the forecasting distribution (Gneiting, 2011, Gneiting and Raftery, 2007, Gneiting et al., 2007). Scoring rules are negatively orientated such that smaller values indicate better performance.

A well-known and popular score is the continuous ranked probability score (CRPS), assessing calibration and sharpness simultaneously. For a predictive distribution $F(y)$ and observation $y_{\text{obs}}$ the CRPS is given as

[TABLE]

where $\bm{1}(y\geq y_{\text{obs}})$ equals $1$ if $y\geq y_{\text{obs}}$ and [math] otherwise, see also Wilks (2011, Sect. 8.5.1).

If $\mu_{F}$ and $\sigma_{F}^{2}$ denote the mean and variance of $F$, the Dawid and Sebastiani (1999) score is given by

[TABLE]

see also Gneiting and Katzfuss (2014).

3.2 Visual Assessment

To visually assess calibration of a probabilistic forecast, verification rank histograms and PIT histograms are employed (Wilks, 2011).

Here, the verification rank histogram (VRH) or Talagrand diagram is used to assess a forecast ensemble $x_{1},\ldots,x_{M}$. It can be obtained by computing the rank of the observation $y$ within the ensemble (for each forecast case). If the ensemble members $x_{1},\ldots,x_{M}$ and the observation $y$ are statistically indistinguishable (exchangeable), the rank of the observation with respect to the ensemble members has a discrete uniform distribution on $\{1,\ldots,M+1\}$. The VRH then plots the empirical frequency of the observation ranks.

To assess calibration of a full predictive probability distribution, the frequencies of the Probability Integral Transform (PIT) values are plotted in equidistant bins. An observation $y$ can be interpreted as a random sample from the “true” distribution $F$ for the respective weather quantity. If the predicted distribution $F_{0}$ is identical to $F$, then $p=F_{0}(y)$ can be considered as realization of a uniformly distributed random variable on $[0,1]$ and the plot of the frequencies of the PIT values $p$ results in a uniform histogram.

3.3 Further Verification Measures

The variance of the PIT values provides further information on the dispersion properties of the predictive distribution, a neutral dispersion being indicated by a variance equal to $\frac{1}{12}=0.0833$, the variance of the uniform distribution on $[0,1]$, see (Gneiting and Ranjan, 2013).

The root mean variance (RMV) is used as a sharpness measure of predictive probability distributions. A main principle of probabilistic forecasting is “maximizing the sharpness of the predictive distribution subject to calibration” Gneiting and Katzfuss (see e.g. 2014), therefore the sharpness should be investigated in conjunction with the calibration.

3.4 Testing for Improvement in Predictive Performance

The statistical relevance of improvement in the verification scores may be investigated by testing for equal predictive performance of the two considered methods with the Diebold-Mariano test for time series, see Gneiting and Katzfuss (2014).

Let $s_{1}(t),s_{2}(t)$ denote the time series of score values obtained from two competing methods as for example EMOS and AR-EMOS, for a verification period of length $T$ (say). Then the large-sample standard normal test statistic adapted from Diebold and Mariano (1995) is given as

[TABLE]

where

[TABLE]

is the average CRPS differential and

[TABLE]

are the empirical autocovariances.

4 Application to ECMWF temperature forecasts

4.1 Data Description and Data Preprocessing

The considered data for our case study contains the forecast ensemble with 50 members of the European Centre for Medium-Range Weather Forecasts (ECMWF, see e.g. Molteni et al., 1996, Buizza et al., 2007). The data consists of $24$-h, $48$-h and $72$-h ahead forecasts initialized at 12 UTC for 2-m surface temperature in Germany along with the verifying observations at 187 different stations in the time period ranging from 2002-01-01 to 2014-03-20, see also (Hemri et al., 2014). In addition, there is one high-resolution forecast and one control forecast.

For the application of time series methods it is of importance to investigate whether the dates appear in chronological order and if some dates are missing. From the full data set with 187 stations only those stations are retained, which do not reveal NA gaps longer than 1. There are 76 stations which do match this rather strict specification. The new data set still contains missing values.

The remaining missing values can in R be replaced by values obtained from linear interpolation using the function na.approx from the package zoo (Zeileis and Grothendieck, 2005).

Figure 1 shows the 76 stations retained for the subsequent analysis, where the station map was produced with the R package ggmap (Kahle and Wickham, 2013). In the subsequent analysis the station Magdeburg in Eastern Germany will be considered for illustration purposes in the case study.

To find possible outliers in the data, the test statistic from Chang et al. (1988) for detecting additive outliers is applied, being implemented as function detectAO in the R package TSA (Chan and Ripley, 2012). It requires the fit of an ARIMA model to the series, which can e.g. be achieved by the function auto.arima from the R package forecast (Hyndman and Khandakar, 2008).

The analysis reveals one indicated suspicious observation. From further inspection, the detected value of $30.6\,^{\circ}\mathrm{C}$ at 2011-08-26 appears a bit too high with respect to preceding and succeeding observations, but may nonetheless be a true value. Thus, it is concluded that the observation should not be altered or removed.

In order to find strong outlying observations, the significance level in detectAO is put to the very small value $\alpha=0.00001$. By applying the above procedure to each station, it is found that 6 stations reveal suspicious values. From these, only two observations are removed from the data set. At date 2003-01-20 and station Hannover the temperature observation is $-90.8\,^{\circ}\mathrm{C}$, which is clearly impossible. At date 2002-11-23 and station Nuernberg, the temperature observation is $-20.1\,^{\circ}\mathrm{C}$, which is very unusual with respect to preceding and succeeding temperature values and, in addition, is by far the smallest value in the complete series.

Removal is done by setting the outlying value to NA and then applying linear interpolation.

4.2 Comparison of EMOS and Heteroscedastic AR-EMOS

First, the state-of-the-art EMOS model is compared to the heteroscedastic AR-EMOS model (called AR-EMOS in the following) presented in Section 2.4, where each model is based on the 50 regular ECMWF 24-h ahead ensemble forecasts. The parameters of the postprocessing models are estimated station-wise, based only on the data available for the respective station (so-called local approach).

A first analysis investigates the performance of the different models at the station Magdeburg in Eastern Germany. Then, in a second step the analysis is carried out for all 76 stations in the data set.

4.2.1 Settings

To fit the AR-EMOS model presented in Section 2.4 (with the basic methodology in Section 2.1), first the model parameters $\alpha_{m},\beta_{m,1},\ldots,\beta_{m,p_{m}}$ for each member $x_{m}$ are estimated by fitting an $\text{AR}(p_{m})$ process to the observed error series $\{z_{m}\}$ from a training period by Yule-Walker estimation as carried out by the function ar, see also Shumway and Stoffer (2006, Section 3.6). The order $p_{m}$ is automatically selected by a modified Akaike information criterion (AIC) proceeding as if the required variance estimate were obtained by maximum likelihood, cf. R Core Team (2018, function ar).

This procedure is invoked with a default training length of $90$ days, which has been found appropriate by Möller and Groß (2016).

To obtain the predictive mean $\mu(t)$ and standard deviation $\sigma(t)$ (as defined in Equation (5)) a second training period is required to estimate the weight parameter $w$ used in the heteroscedastic variance model (5) such that the average CRPS with respect to the predictive distribution $N(\mu(t),\sigma(t)^{2})$ is minimized for the training period. In the case study following later, the default for the additional training period is set to 30 days length.

The EMOS model is fitted with the R package ensembleMOS. Estimating parameters of EMOS usually requires a training period of length between $20$ to $40$ days (Gneiting, 2014, Sect. 4), where in the subsequent study $30$ days are chosen.

4.2.2 Comparison at a Single Station

The plots displayed in Figure 2 show a comparison of the predictive mean and standard deviation of the EMOS and (heteroscedastic) AR-EMOS predictive distribution for Magdeburg. The predictive means of both methods exhibit a strong relationship (squared correlation coefficient equal to 0.9963), while the predictive standard deviations differ to a certain extent (squared correlation coefficient equal to 0.5142). So, although, AR-EMOS proceeds in a quite different way to estimate the predictive mean, the result does apparently not differ from EMOS very much. The different approaches to estimate the variance obviously also yield different results, the AR-EMOS standard deviations have a tendency to be smaller than those of EMOS. The Boxplot for the stand deviations also show that for EMOS there is much more variation in the estimated standard deviations than for AR-EMOS.

Table 1 presents the CRPS, the DSS, the root mean variance (RMV) and the PIT variance for both methods at station Magdeburg. The top row in Figure 3 additionally shows the PIT histograms of both methods at station Magdeburg.

The PIT values of both models have a variance greater than $1/12=0.0833$, indicating underdispersion of the predictive distributions. This underdispersion is visible in the respective PIT histograms as well. However, the PIT variance of AR-EMOS is much closer to $1/12$ than the PIT variance of EMOS. When looking at the PIT histograms, the EMOS histogram indicates a slightly more pronounced bin for small PIT values, indicating a stronger forecast bias. On the contrary, the EMOS predictive distribution is slightly sharper than the AR-EMOS one, however obviously at an expense in dispersion accuracy.

When looking at the verification scores providing an overall judgment on predictive performance (CRPS, DSS), we can further conclude that AR-EMOS performs better than EMOS with respect to the average CRPS as well as to the average DSS at station Magdeburg.

As the difference in CRPS values is relatively small, we investigate whether AR-EMOS provides a statistically significant improvement in CRPS over state-of-the-art EMOS, by a one-sided Diebold-Mariano test for the alternative

[TABLE]

For station Magdeburg, the resulting p-value is given as 0.01722. Thus, the test shows that the CRPS values of AR-EMOS are (on average) indeed significantly smaller than those of EMOS.

4.2.3 Comparison at all Stations

In a second step the analysis performed for a single example station is carried out for all 76 stations in the data set, and the results are aggregated. For this, a local approach is used, that is, the considered models are estimated at each station individually, resulting in location-specific model parameters based only on the data of a specific station.

When aggregating over all 76 stations, for each station 4341 verification days are available, yielding in total 329916 forecast cases to aggregate over. Table 2 shows the resulting verification metrics for EMOS and AR-EMOS.

The verification scores indicate that the predictive distribution of AR-EMOS has comparable but slightly different properties than EMOS, similar to the comparison at a single station. While CRPS and DSS of AR-EMOS are (slightly) smaller than those of EMOS, the RMV is slightly larger, indicating the AR-EMOS predictive distribution to be less sharp than its EMOS counterpart. However, the sharper distribution of EMOS comes at the expense of calibration, the PIT variance of EMOS is much larger than 1/12, indicating underdispersion, while the PIT variance of AR-EMOS is relatively close to $1/12$, indicating a distribution with dispersion properties close to neutral dispersion. These observations are consistent with the PIT histograms in (the bottom row of) Figure 3, where the EMOS PIT histogram clearly exhibits a U-shape, with indicates underdispersion, while the U-shape is much less pronounced in the AR-EMOS PIT histogram.

To find statistical evidence about the significance of the difference in predictive performance between the two methods, the one-sided Diebold-Mariano test is conducted again, this time for the CRPS time series at each of the 76 stations individually.

The upper left panel of Figure 4 shows the station-wise p-values of the Diebold-Mariano test for EMOS vs. AR-EMOS. Small p-values give statistical evidence for the alternative hypothesis that the values of the AR-EMOS CRPS series are on average smaller than the values of the EMOS CRPS series, thus indicating superior performance of AR-EMOS compared to EMOS. At 31 stations the p-value is $\leq 0.1$, thus indicating superior performance of AR-EMOS, see the left vertical red dashed line in (the upper left panel of) Figure 4.

Large p-values give statistical evidence for the converse alternative hypothesis, thus indicating inferior performance of AR-EMOS compared to EMOS. However, there are only 10 stations with p-value greater than 0.9, see right vertical dashed red line.

Thus, for a large portion of the considered stations, AR-EMOS exhibits significantly better predictive performance than EMOS in terms of the CRPS, while EMOS is significantly superior to AR-EMOS only at a few stations.

4.3 Incorporating the High Resolution Forecast with Group Approach

The ECMWF ensemble also comprises a single high-resolution run, whose importance for statistical postprocessing is described e.g. by Gneiting (2014). As indicated in Section 2.6, extended postprocessing models are now considered, which include the high-resolution forecast.

For the ECMWF data considered here, there are 50 exchangeable (that is statistically indistinguishable) forecast members forming one group, while the high-resolution forecast is regarded as second group due to its different properties.

Within the ensembleMOS package, the group membership of each ensemble forecast can be directly specified. As described in Section 2.6, a straightforward ad-hoc way to implement a 2-group AR-EMOS model is to represent the model parameters as a sum of the two group specific parameters and assign the group-specific parameters (fixed) equal weight.

Furthermore, the SLP combination of EMOS and AR-EMOS proposed in Möller and Groß (2016) is revisited. However, in contrast to the original analysis, here the SLP combination of EMOS and AR-EMOS based on the 50 exchangeable members and the additional high-resolution forecast is investigated (which we called extended models, denoted by EMOS*∗* and AR-EMOS*∗*, respectively).

As additional training data is needed to estimate the weights in the SLP combination, the final number of verification days considered differs from the above analyses comparing only EMOS and AR-EMOS. Here, the results at the station Magdeburg are aggregated over 4251 verification days. When aggregating over all 76 stations (each with 4251 verification days) as well, the results are based on 323076 forecast cases in total.

Results for verification scores at the station Magdeburg are presented in Table 3. It is clearly visible that in terms of CRPS and DSS AR-EMOS*∗* improves over EMOS*∗* to a large extent. The improvement is much more pronounced than in the case where the high-resolution forecast was not incorporated into both models. The SLP*∗* combination of the two models improves both scores even more in comparison to EMOS*∗*.

Concerning sharpness as measured by the RMV the AR-EMOS*∗* model yields the sharpest predictive distribution, with a PIT variance extremely close to $1/12$ at the same time. EMOS*∗* and the SLP combination of both models are less sharp (they exhibit nearly the same level of sharpness), however, while EMOS*∗* as a PIT variance indication underdispersion (larger $1/12$), the SLP*∗* combination also has a PIT variance close to neutral dispersion. Therefore, at the station Magdeburg, the sharpness-calibration properties of AR-EMOS*∗* seem to be extremely appropriate and better than those of the other predictive distributions.

Although the improvement in CRPS of AR-EMOS*∗* compared to EMOS*∗* is much more obvious as in the respective analysis at Magdeburg presented in Table 1, a one-sided Diebold-Mariano test for H_{1}:CRPS_{\text{AR-EMOS{}^{\ast}}}<CRPS_{\text{EMOS{}^{\ast}}} at Magdeburg is performed to investigate the significance of the improvement. The resulting p-value is $0.01233$, showing that AR-EMOS*∗* is indeed performing significantly better than EMOS*∗* at Magdeburg in terms of CRPS.

4.3.1 Illustration of Predictive Distributions

For a visual illustration the three predictive probability density functions (PDFs) corresponding to EMOS*∗, AR-EMOS∗, and SLP∗* (two normals and a mixture of two normals as described in Section 2.7) are graphically displayed for the randomly chosen forecast day 2009-06-24, see Figure 5.

Here, it is clearly visible that the PDF of AR-EMOS*∗* and the SLP*∗* combination are more centered around the observation in comparison to the PDF of EMOS,indicating a higher consistency of the respective distributions with the true observation, thus improved calibration. This observation is in line with the verification results presented in Table 3.

4.3.2 Comparison of all Stations

Next, the above described comparison of EMOS*∗, AR-EMOS∗* and SLP*∗* is conducted for all 76 stations, where again the models are estimated station-wise. Due to the need for additional training data for the SLP combination, the number of verification cases considered again differs from the analysis presented in Table 2 as mentioned for the station Magdeburg. Here, the aggregation in Table 4 is performed over 323076 verification cases (4251 verification dates at each of the 76 stations). The aggregated verification metrics show that AR-EMOS*∗* performs better than EMOS*∗* with respect to the CRPS and DSS. With respect to the CRPS the SLP combination performs best, with respect to DSS it also performs clearly better than EMOS*∗. EMOS∗* has the sharpest predictive distribution in terms of the RMV, while AR-EMOS*∗* and SLP${\ast}$ exhibit a similar level of sharpness. However, the PIT variance of AR-EMOS*∗* is much closer to that of neutral dispersion than EMOS*∗*, having a PIT variance indicating underdispersion.

Again, station-wise p-values of the one-sided Diebold-Mariano test for the CRPS series are computed to investigate whether the improvement in CRPS is significant. The upper right panel of Figure 4 shows the resulting p-values at all 76 stations for the one-sided Diebold-Mariano test for the alternative H_{1}:CRPS_{\text{AR-EMOS{}^{\ast}}}<CRPS_{\text{EMOS{}^{\ast}}}.

Again, the dashed red lines denote the significance levels 0.1 and 0.9, respectively. When incorporating the high-resolution forecast into the models, the number of stations where AR-EMOS*∗* performs significantly better (at significance level 0.1) than EMOS*∗* increases to 41. At 16 stations the p-value is greater than 0.9, indicating that EMOS*∗* performs better than AR-EMOS*∗*.

4.4 Postprocessing for Higher Forecast Horizons

Finally we investigate the performance of the heteroscedastic AR-EMOS model for higher forecast horizons. In many applications, typically 24-h ahead forecasts are investigated, but often higher forecast horizons are not considered. To illustrate the effect, we present results for 48-h and 72-h ahead forecasts.

In Section 2.3 the procedure for applying AR-EMOS to forecast horizons greater than 24-h is explained. EMOS is capable of dealing with other than 24-h ahead forecasts as well, and the forecast horizon considered can explicitly specified within the ensembleMOS package.

Table 5 shows the verification metrics for EMOS and AR-EMOS, based on 48-h and 72-h ahead ensemble forecasts. For 48-h ahead forecasts the verification metrics and PIT histograms are based on a total of 329840 verification cases (4340 verification days for each of the 76 stations), for 72-h ahead forecast, they are based on 329764 verification cases (4339 verification days for each station).

For both forecast horizons it is clearly visible, that AR-EMOS improves on EMOS in terms of CRPS and DSS, with the improvement being even more pronounced for 72-h ahead forecasts. Compared to the results on 24-h ahead forecasts, the improvement of AR-EMOS over EMOS becomes clearer the larger the forecast horizon. For both considered horizons in Table 5, the EMOS predictive distribution is a bit sharper than its AR-EMOS counterpart, however, in each case the PIT variance of EMOS indicates underdispersion to a larger extent than the PIT variance of AR-EMOS.

Figure 6 presents the respective PIT histograms of EMOS and AR-EMOS, where the top panel refers to 48-h ahead, and the bottom panel to 72-h ahead forecasts.

To investigate whether the improvement of AR-EMOS over EMOS is indeed a significant one, again the Diebold-Mariano test is performed at each station for the same one-sided alternative as in the previous paragraphs. The lower panel of Figure 4 displays the resulting p-values for 48-h (left panel) and 72-h (right panel). For both forecast horizons the improvement in predictive performance of AR-EMOS over EMOS in terms of CRPS is highly significant for most of the stations: For 48-h ahead forecasts the p-value is $\leq 0.1$ at 64 stations (and $\leq 0.05$ still at 60 stations), for 72-h ahead forecasts even at 70 stations (and $\leq 0.05$ at 67 stations). In both cases, there is only one station with p-value greater than 0.9. Furthermore, the proportion of stations with an insignificant p-value is extremely small for both forecast horizons.

Although for 24-h ahead forecasts there is a larger portion of stations where AR-EMOS performs significantly better than EMOS, there is also a non-negligible number of stations where neither EMOS was significantly better than AR-EMOS nor the other way around. When moving to higher forecast horizons, the number of stations where AR-EMOS is significantly superior to EMOS increases heavily, and at more and more stations the level of significance even gets smaller. This indicates that the performance of AR-EMOS seems increases in comparison to EMOS for higher the forecast horizon.

5 Conclusion

This follow-up work presents some new features and extensions of the AR-EMOS model introduced by Möller and Groß (2016), and is accompanied by an implementation of the method within an R package called ensAR (Groß and Möller, 2019). The original model for the predictive variance is extended to incorporate the ensemble spread, yielding a heteroscedastic model implicitly accounting for the spread-error correlation, in slightly different way than the EMOS model. The (heteroscedastic) AR-EMOS model allows to fit a predictive distribution to a single ensemble member, as the longitudinal part of the (extended) predictive variance can still be computed for a single ensemble forecast, which is an advantage over standard postprocessing approaches such as EMOS.

Additionally, incorporation of a high-resolution forecast is investigated. The conducted case study indicates that this forecast improves predictive performance to a large extent. To incorporate the high-resolution forecast, an AR-EMOS group model is defined, which follows a somewhat different approach than the EMOS group model. In the case study, only a simple heuristic form of the AR-EMOS group model is considered, which already yields excellent results. However, extensions to a more general and data driven form are relatively straight forward and subject to future research.

Finally, a feature of the AR-EMOS model not discussed in the original work is presented. The model allows to fit predictive distributions based on ensemble forecasts with arbitrary forecast horizons. In the original work a case study based only on 24-h ahead forecasts is presented. However, the AR-EMOS model can postprocess ensemble forecasts with arbitrary forecast horizons. For forecast horizons smaller or equal to 24-h ahead the model can be directly employed without any additional modifications. For horizons larger than 48-h ahead, the model can be applied by adding only one small preprocessing step, namely predicting the days between the last validation date of the forecast and the verification date with the AR model, also used to set up the AR-EMOS method itself.

The conducted case study indicates that for forecast horizons beyond 24-h ahead (with 48-h and 72-h ahead considered as examples) the AR-EMOS performs particularly well, and improves significantly over EMOS. Therefore, the autoregressive postprocessing approach shows potential for extremely accurate prediction at higher forecast horizons.

Acknowledgements

The authors acknowledge support of the research by Deutsche Forschungsgemeinschaft (DFG) Grant No. MO 3394/1-1, and by the Hungarian National Research, Development and Innovation Office under Grant No. NN125679. Furthermore, we are grateful to the European Centre for Medium-Range Weather Forecasts (ECMWF) and the German Weather Service (DWD) for providing forecast and observation data, respectively.