Please note that this documentation is currently under construction, particularly Sections 3 and 4 have not been finished yet. Comments or any kind of feedback is highly appreciated. Please send an e-mail to the authors.

1 Introduction

This document briefly reviews and extends the ideas for verifying decadal prediction presented by Goddard et al. [2013] and describes how the plugin MurCSS of the MiKlip Central Evaluation System (CES) can be used to assess hindcast skill accordingly.

Following Goddard et al. [2013], two key questions can be identified which should guide the quality assessment of a decadal ensemble prediction system:

1. Is the initialized hindcast more accurate compared to the un-initialized projection?
2. Is the models’ ensemble spread representative for the forecast uncertainty?

A way to tackle the first issue is to estimate accuracy of the ensemble mean with the mean squared error skill score (MSESS) and its decomposition, covered in Sec. 2. The second question can be addressed using a particular variant of the continuous ranked probability skill score (CRPSS), see Sec. 3. The basic ideas underlying these two approaches stem from traditional forecast verification, for an overview, see, e.g., [Wilks, 2011] and [Jolliffe and Stephenson, 2012]. The decomposition of the mean squared error skill score goes back to Murphy [1988]. Here, however, we use these methods for forecast verification to assess the quality of an ensemble prediction system by verifying a set of hindcasts against a reanalysis or another data set derived from observations.

In section 2, the assessment of accuracy of the ensemble mean with the mean square error skill score is described. Section 3 discusses the representation of forecast uncertainty with the ensemble spread and Sec. 4 explains input and output of the component MurCSS and gives an example analysis.

2 Accuracy with the mean square error skill score

2.1 Mean square error and skill scores

The accuracy of the ensemble mean can be assessed with the mean squared error (MSE). To evaluate the added value of initialising the ensemble prediction system with an estimated state of the ocean (and the atmosphere), performance is compared an uninitialised set of predictions; in the context of climate research uninitialised predictions are called climate projections or climate simulations and are typically driven by scenarios for the radiative forcing (or emissions) [Kirtman et al., 2013]. This setting calls for using skill scores, which evaluate the performance of one prediction system against a reference prediction, typically a more or less trivial benchmark [e.g., Wilks, 2011]. In weather forecast, this trivial benchmark is the climatological forecast, the persistent forecast or a random forecast. The appropriate benchmark for decadal prediction is the uninitialised prediction (climate projection). For hindcasts, the “climate projection” is called a historical run. A skill score evaluates the gain in performance against the trivial benchmark. The general expression for a skill score is

$${SS}_{ref}=\frac{A-A_{ref}}{A_{perf}-A_{ref}}\cdot 100\%,$$

(1)

and is typically given in percent gain with respect to a reference forecast. Here, $A$ is the value for the accuracy measure (here the MSE, see below) for the prediction system under question, $A_{perf}$ is its value for a perfect prediction and $A_{ref}$ is its value for a reference forecast system. I the prediction system in question obtains the same accuracy as the perfect prediction, it will be assigned 100% skill; skill of 0% or less is interpreted as being equal or worse than the (typically trivial) reference system.

In the following, the mean square error MSE is considered as a more or less intuitive (or at least familiar) score to be the basis of the skill score. The mean square error (MSE) between a hindcast $H$ and observations (reanalysis) $O$ is defined as the sum of the squared deviations of the two at instances (times) $j=1,\dots ,n$, and reads following the notation of Goddard et al. [2013]

$${MSE}_H=\frac{1}{n}\sum _{j=1}^n{\left(H_j-O_j\right)}^2.$$

(2)

The minimum value for the MSE is 0, indicating a (hypothetical) perfect prediction system, leading to $A_{perf}=0$ in Eq. (1). The upper bound of the MSE depends on whether $O_j$ (and $H_j$) is bounded, which is not the case for many continuous variables, leading to no upper bound for the MSE.

Using the MSE (Eq. (2)) as a skill score (Eq. (1)) with the climatological forecast $Ō$ as the reference leads to

$$MSESS\left(H,Ō,O\right)=\frac{MS{E}_H-MS{E}_{Ō}}{0-MS{E}_{Ō}}=1-\frac{MS{E}_H}{MS{E}_{Ō}},$$

(3)

with the MSE of the climatological forecast

$${MSE}_{Ō}=\frac{1}{n}\sum _{j=1}^n{\left(Ō-O_j\right)}^2.$$

(4)

The MSESS represents the improvement in accuracy of the hindcasts over the climatological forecast derived from the observations. It is unbounded for negative values (as hindcasts can be arbitrarily bad) but bounded by 1 (if not expressed as percentages) for positive values. It is a function of the hindcast $H$ (which is the ensemble mean in our case), the reference forecast, here the climatological forecast $Ō$, and the observations $O$ themselves, as they enter both, $MS{E}_H$ and $MS{E}_{Ō}$. A positive value of the MSESS indicates gain in hindcast performance compared to the reference forecast; zero or negative values indicates no gain.

It follows that the MSESS against a general reference forecast, e.g., the uninitialised historical runs $R$, reads

$$MSESS\left(H,R,O\right)=1-\frac{MS{E}_H}{MS{E}_R}.$$

(5)

In the following, decompositions of the MSE (Eq. (2)) and its skill score (Eqs. (4) and (5)) are discussed.

2.2 Decomposing the mean square error

Rewriting of the MSE (Eq. (2)) or the MSESS (Eqs. (4) and (5)) leads to elegant decompositions of the MSE and its skill scores into a few contributing terms with intuitive interpretations.

Consider first the decomposition of the MSE: inserting a zero by adding and substracting the mean hindcast over all predictions issued (all times) $\bar{H}=\frac{1}{n}{\sum }_{j=1}^n{H}_j$ and the mean observation $Ō=\frac{1}{n}{\sum }_{j=1}^n{O}_j$ to the right hand side (rhs) of Eq. (2) yields

$${MSE}_H=\frac{1}{n}\sum _{j=1}^n{\left[\left(H_j-\bar{H}\right)-\left(O_j-Ō\right)+\left(\bar{H}-Ō\right)\right]}^2$$

(6)

. Expanding Eq. (6) into single terms leads to $$\begin{array}{lllllll}\hfill {MSE}_H& =\underset{=s_H^2}{\underbrace{\frac{1}{n}\sum _{j=1}^n{\left(H_j-H\right)}^2}}-\underset{=2s_{HO}}{\underbrace{\frac{2}{n}\sum _{j=1}^n\left[\left(H_j-\bar{H}\right)\left(O_j-Ō\right)\right]}}+\underset{=0}{\underbrace{\frac{2\left(\bar{H}-Ō\right)}{n}\sum _{j=1}^n\left(H_j-\bar{H}\right)}}& \hfill & & \hfill & & \hfill \\ \hfill & \underset{=0}{\underbrace{-\frac{2\left(\bar{H}-Ō\right)}{n}\sum _{j=1}^n\left(O_j-Ō\right)}}+\underset{=s_O^2}{\underbrace{\frac{1}{n}\sum _{j=1}^n{\left(O_j-Ō\right)}^2}}+{\left(\bar{H}-Ō\right)}^2,& \hfill & & \hfill \text{(7)}\end{array}$$

with the third and fourth term being zero because the sum of all anomalies are zero. The first term represents the sample variance of the hindcasts $s_H^2$, the second term is the sample covariance of hindcasts and observations $s_{HO}$, the next-to-last term is the sample variance of the observations $s_O^2$ and the last term is the square of the bias $Bias=\bar{H}-Ō$. A shorthand notation of the above equation is

$${MSE}_H={\left(\bar{H}-Ō\right)}^2+s_H^2+s_O^2-2s_{HO}.$$

(8)

Using $s_{HO}=s_H{s}_O{r}_{HO}$ and thus the sample correlation between the hindcasts and the observations $r_{HO}$, Eq. (8) yields

$${MSE}_H={\left(\bar{H}-Ō\right)}^2+s_H^2+s_O^2-2s_H{s}_O{r}_{HO}$$

(9)

In the following, the MSE decomposition Eq. (9) is plugged into the skill score Eq. (1) to yield the decomposition of the MSESS.

2.3 Decomposing the mean square error skill score

The decomposition Eq. (9) is inserted into the MSESS with the climatological forecast $Ō$ as reference (Eq. (4)). Using ${MSE}_{Ō}=s_O^2$, i.e. the MSE of the climatological forecast being an estimate of the climatological variance, we obtain

$$MSESS\left(H,Ō,O\right)=-{\left[\frac{\bar{H}-Ō}{s_O}\right]}^2-{\left[\frac{s_H}{s_O}\right]}^2+2\frac{s_H}{s_O}{r}_{HO}$$

(10)

Adding and subtracting $r_{HO}^2$ and some transformations the MSESS reads

$$MSESS\left(H,Ō,O\right)=r_{HO}^2-{\left[r_{HO}-\frac{s_H}{s_O}\right]}^2-{\left[\frac{\bar{H}-Ō}{s_O}\right]}^2.$$

(11)

The first term of Eq. (11) is the square of the sample correlation coefficient, a dimensionless measure of the linear relationship between the hindcasts and the observations. The second term is the conditional bias squared, a measure of reliability. The last term represents the bias or unconditional bias; the latter can be removed by analysing anomalies. The MSESS with the bias removed reads thus

$$MSESS\left(H,Ō,O\right)={r_{HO}}^2-{\left[r_{HO}-\frac{s_H}{s_O}\right]}^2,$$

(12)

hence being basically decomposed into correlation and conditional bias. For the hindcast $H$ against the climatological forecast $Ō$, we denote the left-hand-side (lhs) of Eq. (12) as ${MSESS}_H=MSESS\left(H,Ō,O\right)$ and, analogously for another forecast $R$ against the climatological forecast $Ō$ as ${MSESS}_R=MSESS\left(R,Ō,O\right)$.

A similar decomposition can be obtained for a general reference from Eq. (5) by multiplying and subtracting $s_o^2$. After a few transformations, it reads

$$MSESS\left(H,R,O\right)=\frac{{r_{HO}}^2-{\left[r_{HO}-\frac{s_H}{s_O}\right]}^2-{r_{RO}}^2+{\left[r_{RO}-\frac{s_R}{s_O}\right]}^2}{1-{r_{RO}}^2+{\left[r_{RO}-\frac{s_R}{s_O}\right]}^2}.$$

(13)

Thus, the MSESS of the hindcast against a reference forecast can be described with the MSESS of the hindcast against the climatological forecast and the MSESS of the reference prediction against the climatological forecast

$$MSESS\left(H,R,O\right)=\frac{MSES{S}_H-MSES{S}_R}{1-MSES{S}_R}.$$

(14)

In the following, the interpretation of the various terms of this decomposition are discussed.

2.4 Interpretation of the MSESS and its decomposition

2.4.1 MSESS

The MSESS is an integral measure of the improvement in accuracy of the hindcasts $H$ over the climatological forecast $Ō$ or another reference forecast $R$ with respect to the observations (reanalysis) $O$. The MSESS can take arbitrary large negative values (as hindcast performance can be arbitrarily bad) but it is bounded from by 1 for positive values. Positive values indicate a gain in accuarcy for the hindcast compared to the reference; zero or a negative values indicate no skill.

2.4.2 Correlation

The correlation coefficient represents a dimensionless measure of the linear relationship between hindcast (or reference) and the observations, ranging from -1 to 1. In the case of the MSESS comparing the hindcast against the climatological forecast the component MurCSS of the central evaluation system (CES) yields the correlation coefficient $r_{HO}$ between $H$ and $O$. In the case of an analysis of the hindcast against another reference prediction the gain in correlation $r_{HO}-r_{RO}$ is shown.

2.4.3 Conditional Bias

The conditional bias is a measure of the reliability of the hindcast/reference prediction, i.e. it gives information about the distribution of the observation given the hindcast, $p\left(o\mid h\right)$ [see, e.g., Wilks, 2011]. Its meaning can be illustrated with a linear regression of the hindcast against the observations, see Fig. 1. The slope of this regression line can be expressed by the correlation coefficient and the standard deviations $m=s_O∕s_H\cdot r_{HO}$. Thus, the conditional bias gives an indication of the ratio of the variances of hindcast and observation. Or in other words: given the hindcasts, it gives information about the distribution of the observations. Figure 1 gives an example for a positive correlation coefficient.

---

---

For a positive correlation coefficient a positive/negative conditional bias indicates a smaller/greater variance of the hindcast compared to the observations. For a negative correlation this relationship is reversed.

The component MurCSS of the central evaluation system (CES) gives for the analysis of the hindcast $H$ against a reference prediction $R$ the gain in the conditional bias

$$\left|r_{HO}-\frac{s_H}{s_O}\right|-\left|r_{RO}-\frac{s_R}{s_O}\right|$$

(15)

–

3 Evaluating ensemble spread

3.1 Continuous ranked probabiliy score

Additionally to the question whether the ensemble mean is an accurate prediction, it rests to clarify whether the ensemble spread is an appropriate representation of the forecast uncertainty. To this end, Goddard et al. [2013] suggested to employ the continuous ranked probability skill score (CRPSS). The CRPSS is a skill score based on the continuous ranked probability score (CRPS), a quadratic measure in probability space, [e.g., Wilks, 2011].

$$CRPS={\int \; <!--nolimits-->}_{-\infty }^{\infty }{\left[F\left(y\right)-F_0\left(y\right)\right]}^2,$$

(16)

with $F\left(y\right)$ being the cumulative distribution function (CDF) derived from the hindcast and, as the observations are not assigned a probability, $F_0\left(y\right)$ is the Heaviside (or switch) function. The Heavyside function is defined as

$$F_0\left(y\right)=\Theta \left(x-y\right)\left\{\begin{array}{cc}0\hfill & \hfill y<x\\ ,1\hfill & \hfill y\ge x,\end{array}\right.,$$

(17)

and basically jumps from 0 to 1 at the observation $y$. CDFs of hypothetical hindcasts and the Heaviside function are schematically depicted in Fig. 2(left). The other two panels of Fig. 2 illustrates how the CRPS works, the middel panel gives the difference of the two forecasts $F\left(y\right)$ and an observation $F_{𝜃}\left(y\right)$, and the right panel shows the squared difference.

---

---

For this example, we deliberately chose two forecasts distributions with the same mean but different variance to illustrate that the CRPS is sensitive to the variance of the forecasted probability distribution. For an ensemble prediction system, the forecast distribution is derived from a discrete (and often small) ensemble. This can be done by, e.g., deriving an empirical cumulative distribution function (ECDF) based on the ranks of the individual members, or assuming an underlying Gaussian distribution and estimate its mean and variance, or by various other methods, such as ensemble dressing [Broecker and Smith, 2008]. Goddard et al. [2013] suggested to assume the Gaussian hypothesis and use the parametric form of the CRPS derived by Gneiting and Raftery [2007]. For a hindcast ensemble with mean ${Ĥ}_j$ and variance ${\sigma }_{{Ĥ}_j}^2$ at hindcast dates $j$, it reads

$$CRPS\left({Ĥ}_j,{\sigma }_{{Ĥ}_j}^2,O_j\right)={\sigma }_{{Ĥ}_j}\left[\frac{1}{\sqrt{\pi }}-2\varphi \left(\frac{O_j-{Ĥ}_j}{{\sigma }_{{Ĥ}_j}}\right)-\frac{O_j-{Ĥ}_j}{{\sigma }_{{Ĥ}_j}}\left(2\Phi \left(\frac{O_j-{Ĥ}_j}{{\sigma }_{{Ĥ}_j}}\right)-1\right)\right],$$

(18)

with $\varphi \left(.\right)$ and $\Phi \left(.\right)$ denoting the Gaussian probability density and distribution function, respectively. This setting allows to assess the ensemble spread for the hypothetical case of an unbiased forecast: the ensemble mean ${Ĥ}_j$ is set to the climatological mean derived from the observations, ${Ĥ}_j=Ō$.

3.2 Continuous ranked probability skill score

As the goal is to assess whether the ensemble spread is a useful measure for the forecast uncertainty, we compare the ensemble forecast to a reference forecast with the desired spread. Goddard et al. [2013] suggested an unbiased forecast with spread equal to the MSE between the ensemble mean and the observations. The latter is a valid representation of the forecast uncertainty. Hence, the CRPSS for this particular case reads

$${CRPSS}_{ES}=1-\frac{{CRPS}_H\left(Ō,{\sigma }^2,O\right)}{{CRPS}_{MSE}\left(Ō,{MSE}_H,O\right)},$$

(19)

and it basically compares the ensemble spread ${\sigma }^2$ to the ${MSE}_H$. Please note that for the desired case ${\sigma }^2={MSE}_H$, the ratio on the rhs is 1 and the ${CRPSS}_{ES}=0$. Thus, for the case that the ensemble spread compares well to the ${MSE}_H$ (our measure for forecast uncertainty), i.e. for the optimal case, the skill score used here attains 0 not 1, as one would typically expect. Figure 3 illustrates this behaviour based on a simulation study with simple stochastic processes emulating the observations and hindcasts.

---

---

[ENSEMBLE SPREAD SCORE HAS BEEN RECENTLY IMPLEMENTED BUT IS NOT DESCRIBED HERE YET.]

4 Verification with MurCSS

4.1 Input

The plugin MurCSS for the CED calculates the MSESS and its decomposition. This sections describes the various options for the plugin. It is useful to search the data for the analysis beforehand with sol_search to reduce the possibilities of the options. Table 1 list and describes all possible options.

4.2 Output

Here, the output of the plugin is described, for comparing hindcasts against the climatological forecast (Sec. 4.2.1) and hindcasts versus another reference forecast (Sec. 4.2.2).

4.2.1 Hindcast vs. climatological forecast

Table 2 lists the output obtained for the analysis of the hincasts versus the climatological forecast.

4.2.2 Hindcast vs. reference prediction

Table 3 lists the output obtained for the analysis of the hincasts versus another reference forecast, e.g., the uninitialised historical runs.

4.3 Example analysis

This section gives an example analysis with the plugin MurCSS. Suppose, we want to estimate the skill of baseline1 against baseline0 of the near-surface air temperature tas with the HadCrut data set as observations. First, solr_search identifies the respective data series. For baseline1 the command looks like

The output is a list of the paths of all available data files, for example:

/miklip/integration/data4miklip/model/global/miklip/ baseline1/output/
MPI-M/MPI-ESM-LR/decs4e1960/mon/atmos/tas/r1i1p1/
tas_Amon_MPI-ESM-LR_decs4e1960_r1i1p1_196101-197012.nc

With the knowledge of the directory structure we have made the following adjustments for the analysis Baseline1 vs. Baseline0 :

After the analysis is completed the netcdf-files and the produced plots (section 4.2) are available. Some of these plots are depicted in Fig. 4 for the analysis Baseline0/Baseline1 vs. Climatology.

---

---

In the analysis with the standard deviation of the observations as reference the MSESS describes how well the observations can be reproduced. A positive value indicates a good reproduction of the observations with a perfect hindcast at 1 and negative values the opposite. The near-surface air temperature is better represented by Baseline1 especially in the indian ocean and southeast asia with respect to Baseline0. But both models can not reproduce the observations in the pacific ocean sufficiently.

The MSESS can be decomposed into the correlation and the conditional bias (see Sec. 2.3). High postive values of the correlation indicate a good accordance of the model and the observations. High negative values suggest a strong linear relationship as well, but the different signs show a poor accordance of model and observations. The improved MSESS of Baseline1 with respect to Baseline0 is mainly produced by a higher correlation coefficient at most of the grid points.

The observations can be best represented with a conditional bias of 0. In Fig. 4 it can be seen that Baseline1 shows the lowest conditional bias, but there are still some areas, for example the pacific ocean and the nothern atlantic, with a pronounced conditional bias.

Another approach to understand the conditional bias and its relationsship to the standard deviations of the hindcast and the observations is the slope of the linear regression. In Fig. 4 the biasslope -1 is depicted as well. It is noticable that the pattern are quite similar to the pattern of the conditional bias. For a positive correlation coefficient a positive slope and a positive conditional bias suggest a higher standard deviation of the observations with respect to the standard deviation of the hindcast. If the correlation coefficient is negative a positive slope and conditional bias indicate a higher standard deviation of the hindcast. For the near-surface air temperatur the correlation and the conditional bias is positive in Central Europe for Baseline1 and Baseline0, which suggest a higher standard deviation of the observations. The negative correlation and conditional bias over the indian ocean in Baseline0 indicate a higher standard deviation of the observations as well, but in Baseline1 the variations of the hindcast is higher due to the positive correlation and the negative conditional bias.

The ratio of the standard deviations $S_O∕S_H$ is also shown in Fig. 4.

In Fig. 5 the output of MurCSS is depicted for the analysis Baseline1 vs. Baseline0. A positive value of the MSESS suggest a improved accuracy of the hindcast compared to the reference and a negative value indicates the opposite.

---

---

In Fig. 5 a) it can be seen, that Baseline1 is more accurate in most of the regions. Only a few small areas, for example in South America or the Labbrador Sea the accordance of Baseline0 and the observations is higher. Fig. 5 b) shows the gain of the correlation coefficient of Baseline1 compared to Baseline0. Positive values indicate a improved correlation of Baseline1. The analysis of the near-surface air temperature show an improved accuracy of Baseline1 with respect to Baseline0 for almost every gridpoint as well. In Fig. 5 c) the gain of the conditional bias is illustrated. Here it is noticable, that negative values represent a reduction of the bias and thus a improved accuarcy. For example in the Labrador Sea the conditional bias is higher in Baseline0, the correlation is, however, slightly increased. Allover the contribution of the bias overbalance the effect of the correlation, thus the MSESS shows a decreased accuracy in the Labrador Sea.

An additional approach to illustrate the conditional bias is the gain of the biasslope -1, which is shown in Fig. 5 d). At this figure the negative values suggest a reduction of the conditional bias as well. Thus the pattern is quite similar to Fig. 5 c). Due to the calculation of the gain with absolute values ($|m_{HO}-1|-|m_{RO}-1|$) a statement about the gain in the standard deviations is no longer possible.

A third method the illustrate the condiational bias is the bias Skill Score $1-\frac{|m_{HO}-1|}{|m_{RO}-1|}$, Fig. 5 e). This time, positive values suggest a higher conditional bias of Baseline1 compared to Baseline1 and negative values the opposite. Nevertheless, the Figs. 5 c)-e) illustrating the conditional bias have similar patterns.

References