Evaluation of point forecasts for extreme events using consistent scoring functions

Extreme events occur in many systems, from atmospheric to economic, and present significant challenges to society. Hence the accurate prediction of extreme events is of vital importance. In many such situations, competing forecasts are produced by a variety of forecast systems and it is natural to want to compare the performance of such forecasts with emphasis on the extremes.

In this context, it is critical that the methodology for requesting, evaluating and ranking competing forecasting systems is decision-theoretically coherent. Since the future is not precisely known, ideally forecasts ought to be probabilistic in nature, taking the form of a predictive probability distribution and assessed using a proper scoring rule (Gneiting and Katzfuss, 2014). Nonetheless, in many contexts and for a variety of reasons, point forecasts (i.e. single-valued forecasts taking values from the prediction space) are issued and used. For decision-theoretically coherent evaluation of point forecasts, either the scoring function (such as the squared error or absolute error scoring function) that will be used to assess predictive performance should be advertised in advance, or a specific functional (such as the mean or median) of the forecaster’s predictive distribution should be requested and evaluation then conducted using a scoring function that is consistent for that functional (Gneiting, 2011a ). The use of proper scoring rules and consistent scoring functions encourage forecasters to quote honest, carefully considered forecasts.

To compare competing forecasts for the extremes, one may be tempted to use what would otherwise be a proper scoring rule or consistent scoring function, but restrict evaluation to a subset of events for which extremes were either observed, or forecast or perhaps both. However, such methodologies promote hedging strategies (as illustrated in Section 2) and can result in misguided inferences. This gives rise to the forecaster’s dilemma, whereby “if forecast evaluation proceeds conditionally on a catastrophic event having been observed [then] always predicting a calamity becomes a worthwhile strategy” (Lerch et al., 2017).

Nonetheless, for predictive distributions Gneiting and Ranjan, (2011) showed that a suitable alternative exists in the threshold-weighted continuous ranked probability score, which is a weighted version of the commonly used continuous ranked probability score (CRPS). The weight is selected to emphasise the region of interest (such as the tails or some other region of a variable’s range) and induces a proper scoring rule. This technique extends in a natural way to point forecasts targeting the median functional, since the CRPS is a generalisation of the absolute error scoring function. The UK Met Office has recently applied this method to report the performance of temperature forecasts, with emphasis on climatological extremes (Sharpe et al., 2020).

Despite this progress, Lerch et al., (2017) offer this summary of the general situation for evaluating point forecasts at the extremes: “there is no obvious way to abate the forecaster’s dilemma by adapting existing forecast evaluation methods appropriately, such that particular emphasis can be put on extreme outcomes.” In this paper we remedy this situation. We construct consistent scoring functions that can be used evaluate point forecasts that emphasise performance in the region of interest for important classes of functionals, including expectations and quantiles. Moreover, the relevant specific case of these constructions agrees with the extension of the threshold-weighted CRPS to the median functional.

The main idea of this paper can be illustrated using the squared error scoring function $S(x,y)=(x-y)^{2}$, which is consistent for point forecasts of the expectation functional. Suppose that the outcome space $\mathbb{R}$ is partitioned as $\mathbb{R}=I_{1}\cup I_{2}$, where $I_{1}=(-\infty,a)$ and $I_{2}=[a,\infty)$ for some $a$ in $\mathbb{R}$. Corollary 3.2 gives the decomposition $S=S_{1}+S_{2}$, where each scoring function $S_{i}$ is consistent for expectations whilst also emphasising predictive performance on the interval $I_{i}$. In particular, if $x,y\in I_{1}$ then $S_{2}(x,y)=0$. Given a point forecast $x$ and corresponding observation $y$, the explicit formula for $S_{2}$ is

Here $\mathbbm{1}$ denotes the indicator function, so that $\mathbbm{1}\{x\geq a\}$ equals $1$ if $x\geq a$ and $0$ otherwise.

The performance of competing point forecasts for expectations can then be compared by computing the mean scores $\bar{S}$, $\bar{S}_{1}$ and $\bar{S}_{2}$ for each forecast system, over the same set events. To illustrate, consider two forecast systems A and B, whose error characteristics are depicted by the scatter plot of Figure 1. System B is homoscedastic (i.e., has even scatter about the diagonal line throughout the variable’s range) whereas System A is heteroscedastic (with relatively small errors over lower parts of the variable’s range and relatively large over higher parts). For this set of events, the mean squared error $\bar{S}$ for each system is very similar and there is no statistical significance in their difference. However, with $a=10$, the mean score $\bar{S}_{2}$ of System A is significantly higher than that of B (i.e., B is superior when emphasis is placed on the region $[10,\infty)$) whilst the mean score $\bar{S}_{1}$ of System A is significantly lower than that of B (i.e., A is superior when emphasis is placed on the region $(-\infty,10)$). Full details for this case study are given in Section 3.5.

This example is illustrative. Decompositions of the outcome space can consist of more than two intervals, and the boundary between intervals can also be ‘blurred’ by selecting suitable weight functions. Each decomposition of the outcome space then induces a decomposition

of a given consistent scoring function $S$, whose summands $S_{i}$ are also consistent for the functional of concern and with each $S_{i}$ emphasising forecast performance in the relevant region of the outcome space. Such decompositions are presented for the consistent scoring functions of quantiles, expectations, expectiles, and Huber means. The approach is unified, in that the same decomposition of the outcome space induces decompositions of the CRPS and of a variety of scoring functions for point forecasts. Details are given in Section 3 with illustrative examples.

Mathematically, the main result of this paper is a corollary of the mixture representation theorems for the consistent scoring functions of quantiles, expectiles (Ehm et al., 2016) and Huber functionals (Taggart, 2021). This furnishes each member $S_{i}$ of the decomposition (1.2) with an interpretation rooted in optimal decision theory. Certain classes of optimal decision rules elicit a user action if and only if a point forecast $x$ exceeds a particular decision threshold $\theta$. The score $S_{i}(x,y)$ is a measure of economic regret, relative to decisions based on a perfect forecast, of using the point forecast $x$ when the observation $y$ realises, averaged over all decision thresholds $\theta$ belonging to corresponding interval $I_{i}$ of the partition. Precise details are given in Section 4 and illustrated with the aid of Murphy diagrams.

We illustrate how a seemingly natural approach for comparing predictive performance at the extremes creates opportunities for hedging strategies.

A meteorological agency maintains two forecast systems A and B, each of which produces predictive distributions for hourly rainfall $Y$ at a particular location. System A is fully automated and hence cheaper to support than B. System B has knowledge of the System A forecast prior to issuing its own forecast. The agency requests the mean value of their predictive distributions with a lead time of 1 day. These point forecasts are assessed using the squared error scoring function, which is consistent for forecasts of the mean (i.e., the forecast strategy that minimises one’s expected score is to issue the mean value of one’s predictive distribution). For a two year period, the bulk of observed and forecast values fall in the interval $[0\,\mathrm{mm},10\,\mathrm{mm}]$ and there is no statistically significant difference between the performance of the two systems when scored using the squared error function.

However, the maintainers of System B claim that B performs better for the extremes and that bulk verification statistics fail to highlight this. The agency decides to use forecasts from A for the majority of cases, but will test whether B is significantly better than A at forecasting the heavy rainfall events. The agency considers four options for selecting hourly events to assess, after which the squared error scoring function will be used to compare predictive performance on those events. If System B does not perform significantly better than A over the next 12 months then it will be decommissioned.

For a given forecast case, let $F_{\mathrm{A}}$ and $F_{\mathrm{B}}$ denote the predictive distributions produced by each system, and $x_{\mathrm{A}}$ and $x_{\mathrm{B}}$ the respective point forecasts issued. Suppose that the random variable $Y$ has a distribution specified by $F_{\mathrm{B}}$. For each option, the maintainers of System B have a strategy to optimise their expected score; that is, to choose $x_{\mathrm{B}}$ such that $\mathbb{E}[(x_{\mathrm{B}}-Y)^{2}]$ is minimised.

1. Option 1:

   Only assess events where the observation is at least 20 mm.
   Strategy: If $\mathbb{P}(Y\geq 20)>0$ then $F_{\mathrm{B}}|\{Y\geq 20\}$, which denotes the predictive distribution of $\mathrm{B}$ conditioned on the event $\{Y\geq 20\}$, exists. In this case, forecast $x_{\mathrm{B}}=\mathrm{Mean}(F_{\mathrm{B}}|\{Y\geq 20\}))$ and otherwise forecast $x_{\mathrm{B}}=20$.

2. Option 2:

   Only assess events where either $x_{\mathrm{A}}$ or $x_{\mathrm{B}}$ is at least 20 mm.
   Strategy: If $\max(x_{\mathrm{A}},\mathrm{Mean}(F_{\mathrm{B}}))\geq 20$ then $x_{\mathrm{B}}=\mathrm{Mean}(F_{\mathrm{B}})$. Otherwise if $\mathbb{E}[(20-Y)^{2}]<\mathbb{E}[(x_{\mathrm{A}}-Y)^{2}]$ then $x_{\mathrm{B}}=20$ else $x_{\mathrm{B}}=\mathrm{Mean}(F_{\mathrm{B}})$. This will ensure that the event is assessed whenever System B expects that a forecast of 20 will receive a more favourable score than than a forecast of $x_{\mathrm{A}}$.

3. Option 3:

   Only assess events where either $x_{\mathrm{A}}$, $x_{\mathrm{B}}$ or the observation is at least 20 mm.
   Strategy: If $\max(x_{\mathrm{A}},\mathrm{Mean}(F_{\mathrm{B}}))\geq 20$ then $x_{\mathrm{B}}=\mathrm{Mean}(F_{\mathrm{B}})$. Else if $\mathbb{E}[(20-Y)^{2}]<\mathbb{E}[(x_{\mathrm{A}}-Y)^{2}]$ then $x_{\mathrm{B}}=20$. Otherwise, the only other way the event will be assessed is if the observation is at least 20mm, so forecast $x_{\mathrm{B}}=19.9$.

4. Option 4:

   Only assess events where $x_{\mathrm{A}}\geq 20$.
   Strategy: In this case there is nothing that System B can to do influence which events will be assessed. Therefore $x_{\mathrm{B}}=\mathrm{Mean}(F_{\mathrm{B}})$.

Of these, only Option 4 does not expose the agency to a hedging strategy from System B. However, under Option 4 the developers of A may employ the strategy of forecasting $x_{\mathrm{A}}=\max(19.9,\mathrm{Mean}(F_{\mathrm{A}}))$, so that no further comparison of systems is made.

There is an Option 5: use a scoring function that is consistent for the mean functional and that emphasises performance at the extremes. We turn attention to this now.

We work in a setting where point forecasts $x$ and observations $y$ belong to some interval $I$ of the real line $\mathbb{R}$ (possibly with $I=\mathbb{R}$). In Section 3.1 we introduce partitions of unity, which are used to ‘subdivide’ the outcome space $I$ into subregions of interest. We then illustrate in Section 3.2 how such subdivisions induce decompositions of the CRPS, where each member of the decomposition emphasises performance on the corresponding subregion of $I$ whilst retaining propriety. To obtain analogous decompositions for consistent scoring functions, we recall in Section 3.3 that the consistent scoring functions for quantiles and expectations (among others) have general mathematical forms. The aim is to find which specific instances of the general form emphasise performance on the subregions specified by the partition of unity. This is answered in Section 3.4. Section 3.5 contains examples of such decompositions, and opens with a worked example showing how to find the formula for the squared error decomposition of Equation (1.1).

Mixture representations and elementary scoring functions facilitate a decision-theoretic interpretation of the scoring function decompositions of Corollary 3.2.

To avoid unimportant technical details, in this section assume that $I=\mathbb{R}$, $g$ is a nondecreasing differentiable function on $\mathbb{R}$ and $\phi$ is convex twice-differentiable function on $\mathbb{R}$. Suppose that $S$ is any one of the scoring functions $\mathrm{QSF}(g,\alpha)$, $\mathrm{ESF}(\phi,\alpha)$ or $\mathrm{HSF}(\phi,\nu)$. Thus $S$ is consistent for some functional $\mathrm{T}$ (either a specific quantile, expectile or Huber mean). Then $S$ admits a mixture representation

where each $S_{\theta}^{\mathrm{T}}$ is an elementary scoring function for $\mathrm{T}$ and the mixing measure $\mathrm{d}M(\theta)$ is nonnegative (Ehm et al., 2016; Taggart, 2021). That is, $S$ is a weighted average of elementary scoring functions. Explicit formulae for these elementary scoring functions and mixing measures are given in Tables 2 and 3.

The elementary scoring functions $S_{\theta}^{\mathrm{T}}$ and their corresponding functionals $\mathrm{T}$ arise naturally in the context of optimal decision rules. In a certain class of such rules, a predefined action is taken if and only if the point forecast $x$ for some unknown quantity $Y$ exceeds a certain decision threshold $\theta$. The usefulness of the forecast for the problem at hand can be assessed via a loss function $S_{\theta}$, where $S_{\theta}(x,y)$ encodes the economic regret, relative to a perfect forecast, of applying the decision rule with forecast $x$ when the observation $y$ is realised. Typically $S_{\theta}(x,y)>0$ whenever $\theta$ lies between $x$ and $y$ (e.g., the forecast exceeded the decision threshold but the realisation did not, resulting in regret), and $S_{\theta}(x,y)=0$ otherwise (i.e., a perfect forecast would have resulted in the same decision, resulting in no regret).

For the decision rule to be optimal over many forecast cases, the point forecast $x$ should be one that minimises the expected score $\mathbb{E}S_{\theta}(x,Y)$, where $Y\sim F$ for the predictive distribution $F$. Binary betting decisions (Ehm et al., 2016) and simple cost–loss decision models (where the user must weigh up the cost of taking protective action as insurance against an event which may or may not occur (Richardson, 2000)) give rise to some $\alpha$-quantile of $F$ being the optimal point forecast. The $\alpha$-expectile functional gives the optimal point forecast for investment problems with cost basis $\theta$, revenue $y$ and simple taxation structures (Ehm et al., 2016). The Huber mean generates optimal point forecasts when profits and losses are capped in such investment problems (Taggart, 2021). Though the language here is financial, such investment decisions can be made using weather forecasts (Taggart, 2021, Example 5.4). In each case, the particular score $S_{\theta}(x,y)$ is, up to a multiplicative constant, the elementary score $S_{\theta}^{\mathrm{T}}(x,y)$ in Equation (4.1) for the appropriate functional $\mathrm{T}$.

Thus the representation (4.1) expresses $S(x,y)$ as a weighted average of elementary scores, each of which encodes the economic regret associated with a decision made using the forecast $x$ with decision threshold $\theta$. The weight on each $\theta$ is determined by the mixing measure $\mathrm{d}M(\theta)$, as detailed in Table 3. Relative to $\mathrm{QSF}(g,\alpha)$, the scoring function $\mathrm{QSF}(g_{j},\alpha)$ weighs the economic regret for decision threshold $\theta$ by a factor of $\chi_{j}(\theta)$, thus privileging certain decision thresholds over others.

We illustrate these ideas using the point forecasts generated by Systems A and B from the Synthetic data example (Section 3.5). These point forecasts target the mean functional, and so induce optimal decision rules for investment problems with cost basis $\theta$. The decision rule is to invest if and only if the forecast future value $x$ of the investment exceeds the fixed up-front cost $\theta$ of the investment. As previously noted, neither forecast system performs significantly better than the other when scored by squared error. However, for certain decision thresholds $\theta$, the mean elementary score $\bar{S}_{\theta}$ (which measures the average economic regret) of one system is clearly superior to the other. This is illustrated by the Murphy diagram of Figure 6, which is a plot of $\bar{S}_{\theta}$ against $\theta$ (Ehm et al., 2016). Since a lower score is better, users with a decision threshold exceeding about $8$ should use forecasts from B, while those with a decision threshold less than $8$ should use forecasts from A.

Let $S=S_{1}+S_{2}$ denote the same decomposition of the squared error scoring function used for the Synthetic data example (Section 3.5), induced from a rectangular partition of unity. To aid simplicity of interpretation, assume for the user group under consideration that there is one constant $k$ such that the each elementary score equals $k$ multiplied by the economic regret. We observe the following.

For $S$, $\phi(t)=2t^{2}$ and so $\mathrm{d}M(\theta)=4\,\mathrm{d}\theta$ by Table 3. So by Equation (4.1), the mean squared error $\bar{S}$ for each forecast system is proportional to the total area under the respective Murphy graph. Hence mean squared error $\bar{S}$ can be interpreted as being proportional to the mean economic regret averaged across all decision thresholds.

For $S_{2}$, $\mathrm{d}M(\theta)=4\chi_{2}(\theta)\,\mathrm{d}\theta$, where $\chi_{2}$ is the weight function $\chi_{R}$ on the top panel of Figure 6. From Equation (4.1), the mean score $\bar{S}_{2}$ for each forecast system is proportional to the area under the respective Murphy curve for decision thresholds $\theta$ in $[10,\infty)$. This is illustrated on the top panel of Figure 6 with shaded regions. System B has the smaller area and mean score and is hence superior. So $\bar{S}_{2}$ can be interpreted as being proportional to the mean economic regret averaged across all decision thresholds $\theta$ in $[10,\infty)$. Similarly, $\bar{S}_{1}$ can be interpreted as being proportional to the mean economic regret averaged across all decision thresholds $\theta$ less than $10$.

The bottom panel of Figure 6 illustrates the effect of using a trapezoidal weight function $\chi_{T}$. The corresponding scoring function applies zero weight to economic regret experienced by users with decision thresholds less than 0, and increasing weight to thresholds greater than zero until a full weight of 1 is applied for decision thresholds greater than or equal to 20. The areas of the shaded regions are proportional to the mean score for Systems A and B.

We have shown that the scoring functions that are consistent for quantiles, expectations, expectiles or Huber means can be decomposed as a sum of consistent scoring functions, each of which emphasises predictive performance in a region of the variable’s range according to the selected partition of unity. In particular, this allows the comparison of predictive performance for the extreme region of a variable’s range in a way that cannot be hedged and is not susceptible to misguided inferences.

The decomposition is consonant with the analogous decomposition of the CRPS using the threshold-weighting method, and hence by extension with known decompositions for the absolute error scoring function.

Such decompositions are a consequence of the mixture representation for each of the aforementioned classes of consistent scoring functions and their associated functionals. This has two implications. First, each score in the decomposition can be interpreted as the weighted average of economic regret associated with optimal decision rules, where the average is taken over all user decision thresholds and the weight applied is given by the corresponding weight function in the partition of unity. Second, such decompositions will also exist for the consistent scoring functions of other functionals that have analogous mixture representations.

The author thanks Jonas Brehmer, Deryn Griffiths, Michael Foley and Tom Pagano for their constructive comments, which helped improve the quality of this paper.