L-moments for automatic threshold selection in extreme value analysis

Abstract

In extreme value analysis, sensitivity of inference to the definition of extreme event is a paramount issue. Under the peaks-over-threshold approach, this translates directly into the need of fitting a Generalized Pareto distribution to observations above a suitable level that balances bias versus variance of estimates. Selection methodologies established in the literature face recurrent challenges such as an inherent subjectivity or high computational intensity. We suggest a truly automated method for threshold detection, aiming at time efficiency and elimination of subjective judgment. Based on the well-established theory of L-moments, this versatile data-driven technique can handle batch processing of large collections of extremes data, while also presenting good performance on small samples. The technique’s performance is evaluated in a large simulation study and illustrated with significant wave height data sets from the literature. We find that it compares favorably to other state-of-the-art methods regarding the choice of threshold, associated parameter estimation and the ultimate goal of computationally efficient return level estimation.

Appendixes

1.1 Appendix 1: ALRSM distances for the hybrid illustrative sample

See Table 6.

Table 6 Automatic selection procedure for a simulated sample of size $n=500$ from a Hybrid($0.75,\,0.2$) distribution, with sample quantile threshold candidates $\left\{ u_i\right\} _{i=1}^{20}$, as in (6) : minimum Euclidean distance $d_{u_i}$ between each point $( t _{3,u_i}, t _{4,u_i})$ for the excesses over each candidate threshold $u_i$ and the GPd theoretical curve (4). Selected threshold, $i=14$, in bold

1.2 Appendix 2: Complementary simulation study: GEVd data

In the interest of briefly analysing the performance of the proposed methodology against the presented alternatives in a broader non-GPd scenario, we performed a complementary simulation study to that presented in Sect. 3, using samples drawn from distributions without an exact GPd tail. For this succinct presentation, we considered two GEVd with tail weights given by two of the previously used EVIs—again, both heavy and light tail scenarios were treated. The GEV cumulative distribution function is generally given in terms of its location $\mu \in \mathbb {R}$, scale $\sigma >0$ and shape $\xi \ne 0$ parameters as ${\text {GEV}}(x|\mu ,\,\sigma ,\,\xi ):= \exp \left\{ -\left( 1+\xi \,\frac{x-\mu }{\sigma }\right) ^{-\nicefrac {1}{\xi }}\right\}$, while for the case $\xi =0$ it is taken as the limit of this expression. The remaining parameters were procured as to have equivalence, after some point, of the new distributions with the tail of the Hybrid($0.75,\,\xi$) distributions previously used—this serves the purpose of comparability of results under both studies. Thus, location and scale of the GEVd were defined as functions of $\xi$, respectively as

$$\begin{aligned} \mu =0.75+\frac{1-0.75}{\xi }\left[ \left( \frac{1}{1-0.75}\right) ^{-\xi }-1\right] \;, \qquad \sigma =\left( 1-0.75\right) ^{\xi +1}\;. \end{aligned}$$

(11)

Note, however, that the GEVd does not carry a true threshold, as the Hybrid did, and as such assessment of selection results is more problematic: although we show the mean selected thresholds for each case in Fig. 15, given the lack of interpretability for this case, we choose to focus mainly on the selections’ performance regarding estimation of the EVI and $p=0.001$ exceedance probability quantile.

Analogously to the scheme presented before, estimated entities were computed from a set of 1500 samples generated according to each combination of the proposed scenarios:

• GEVd tail shape parameter: $\xi =-0.2$; $\xi =0.2$;

• Sample size: $n=1000$; $n=500$; $n=200$;

• Sets of candidate thresholds: $I=10$ as in (5); $I=20$ as in (6).

Fig. 15

Mean selected threshold for the 1500 simulations: $n=1000,\,500,\,200$ observations from the GEV($\mu ,\,\sigma ,\,\xi$) with $\xi =-0.2,\,0.2$, $\mu$ and $\sigma$ as in Eq. (11) and sample quantile threshold candidates $\left\{ u_i\right\} _{i=1}^{I}$, $I=10,\,20$

Analysis of the results comprised in Figs. 15 and 16 yields very similar conclusions to the study presented in the main text (colour scheme remains as green for the ALRSM, red for the aSTSM and yellow for the SGFSM). We firstly observe, from Fig. 15, that the mean selected threshold by the ALRSM is consistently higher than that of the competing SGFSM and aSTSM for every scenario posed. In terms of proximity of estimated EVI and 99.9% quantile to their true values, the proposed methodology clearly shows superior performance, in what regards the smaller bias values depicted in Fig. 16a and closer to 1 values for the ratios shown in Fig. 16c. Furthermore, the RMSE associated with estimation of these parameters also appears to be competitive against the alternative selection methodologies, as can be assessed in Fig. 16b, d, resp., corresponding to estimation of the EVI and 99.9% quantile.

In Tables 7 and 8 we can see that the conclusions in Sect. 3 concerning the computational efficiency of the methodologies remain true when dealing with these samples: the ALRSM is, among the three considered, by far the fastest methodology to perform the complete 1500 simulations, and although it shows some failures (2.6% in the worst case), as previously mentioned, these are due to the fitting methodology used, contrary to the selection failures of the SGFSM—which is, once again, the most direct competitor to the ALRSM in terms of bias and RMSE of estimates.

Thus, general conclusions regarding validity of our automated method remain valid for these additional scenarios, even without the possibility of evaluating the threshold selection in itself.

Fig. 16

a Mean bias and b RMSE of EVI estimation, c mean ratio of estimates over true 99.9% quantile $\left( \nicefrac {\widehat{\chi _{0.001}}}{\chi _{0.001}}\right)$ and d RMSE of 99.9% quantile estimation for the 1500 simulations: $n=1000,\,500,\,200$ observations from the GEV($\mu ,\,\sigma ,\,\xi$) with $\xi =-0.2,\,0.2$, $\mu$ and $\sigma$ as in Eq. (11) and sample quantile threshold candidates $\left\{ u_i\right\} _{i=1}^{I}$, $I=10,\,20$

Table 7 Total run time, in seconds, of threshold selection and parameter estimation for 1500 simulated samples with size $n=1000,\,500,\,200$ from the GEV($\mu ,\,\sigma ,\,\xi$) with $\xi =-0.2,\,0.2$, $\mu$ and $\sigma$ as in Eq. (11) and sample quantile threshold candidates $\left\{ u_i\right\} _{i=1}^{I}$, $I=10,\,20$; best run times for each scenario in bold

Table 8 Total failure count for 1500 simulated samples with size $n=1000,\,500,\,200$ from the GEV($\mu ,\,\sigma ,\,\xi$) with $\xi =-0.2,\,0.2$, $\mu$ and $\sigma$ as in Eq. (11) and sample quantile threshold candidates $\left\{ u_i\right\} _{i=1}^{I}$, $I=10,\,20$

1.3 Appendix 3: Complementary simulation study: quantized data

A known issue in estimation under the GPd is, as mentioned before, the effect of data quantization. Several authors have addressed this problem more or less extensively, from different perspectives which we can find amongst the already referred works. For example, Bader et al. (2018) attempt to “fix” the quantized data by applying a jittering procedure previous to the application of their threshold selection procedure—the objective is to approximate the data to that obtained from a continuous process, in order to use the known results for continuous distributions without further issues. On the other hand, other works such as that of Deidda (2010) choose to develop methods that can deal with quantized samples in themselves, that is, that perform their task well whether the data is discretized or not. Although it is our view that the first approach would be the most indicated to deal with the problem, in this brief simulation study we will assess the robustness of the proposed ALRSM in the face of considerably discretized data. However, we limit the investigation to rounded samples from the Hybrid distribution, for the reasons presented in Sect. 3.

For illustration purposes, consider again the Hybrid($0.75,\,0.2$) sample of dimension $n=500$ used in the examples of Sects. 2 and 3, Figs. 3 and 4 and Table 1. The data was used with the machine precision of the R software. At this point, we round the data points to 2 decimal places, and consider this to be our quantized sample, with a quantization level $\Delta =0.01$, which we feel is an appropriate representation for a variety of real world discretized data sets (note that the real data sets used in Sect. 4 were recorded to 3 decimal places). Figure 17 is parallel to Fig. 3 of the main text: it shows the LMRD and the selection process for the same quantized sample, and as we can see the results are analogous, hinting at robustness of the ALRSM against quantization. As before, the heuristic’s distances are comprised in Table 9.

Fig. 17

Automatic selection procedure for a simulated sample of size $n=500$ from a Hybrid($0.75,\,0.2$) distribution, quantized with $\Delta =0.01$, sample quantile threshold candidates $\left\{ u_i\right\} _{i=1}^{20}$, as in (6). Theoretical L-skewness and L-kurtosis for continuous GPd with same $\xi =0.2$ plotted as hollow triangle; sample estimates $( t _{3,u_i}, t _{4,u_i})\,,\,i=1,\ldots ,20$, plotted as hollow circles, with corresponding minimum-distance points on the curve as solid circles; $( t _{3,u^*}, t _{4,u^*})$ for selected threshold as solid red

Table 9 Automatic selection procedure for a simulated sample of size $n=500$ from a Hybrid($0.75,\,0.2$) distribution, quantized with $\Delta =0.01$, with sample quantile threshold candidates $\left\{ u_i\right\} _{i=1}^{20}$, as in (6) : minimum Euclidean distance $d_{u_i}$ between each point $( t _{3,u_i}, t _{4,u_i})$ for the excesses over each candidate threshold $u_i$ and the GPd theoretical curve (4). Selected threshold, $i=14$, in bold

Analogously to the study in Sect. 3, we will more thoroughly evaluate performance of the methodologies in selecting the threshold and estimating the EVI and 99.9% quantile, regarding bias and RMSE of estimation, computed from the same set of 1500 Hybrid($0.75,\,\xi$) samples, now subjected to rounding (to the second decimal place, as mentioned—$\Delta =0.01$), generated according to each combination of the proposed scenarios:

• GP tail shape parameter: $\xi =-0.2$; $\xi =0.2$;

• Sample size: $n=1000$; $n=500$; $n=200$;

• Sets of candidate thresholds: $I=10$ as in (5); $I=20$ as in (6).

Looking at Fig. 18 we can see that the ALRSM again outperforms the competitors in terms of threshold selection bias and RMSE across the considered sample sizes and number of candidates, as well as for both heavy and bounded tails. From comparison with Fig. 5 of the main text, the suggested methodology appears to be robust to quantization of data in this point, while the SGFSM seems to show the most impacted selection performance, being outdone by even the aSTSM in several cases.

Fig. 18

a Mean threshold selection bias and b RMSE of threshold selection for the 1500 simulations: $n=1000,\,500,\,200$ observations from a Hybrid($0.75,\,\xi$), with $\xi =-0.2,\,0.2$, quantized with $\Delta =0.01$, and sample quantile threshold candidates $\left\{ u_i\right\} _{i=1}^{I}$, $I=10,\,20$

Concerning EVI estimation, Fig. 19a, b shows once more close similarity to the bottom and middle panels of Fig. 6, with common underestimation of the shape parameter to the three selection methodologies. Again, both in terms of bias and RMSE, the ALRSM performs either better or comparably well to the alternatives for every considered scenario—EVI estimation through these selections does not appear here to be very sensitive to the quantization level applied.

For the heavy tail samples ($\xi =0.2$), estimation of the 99.9% quantile seems to have suffered some impact from the discretization, when we compare Fig. 19c with Fig. 7b. Here, the ALRSM seems to be slightly more sensitive to the rounding of data than the SGFSM. However, it still performs satisfactorily well, specially for the bounded tail case and regarding robustness against sample size and number of candidate thresholds. Estimation RMSE for this return level, seen in Fig. 19d, is also comparable between the 3 processes, across the proposed scenarios and similar to the depicted in Fig. 8b for the continuous case.

Fig. 19

a Mean bias and b RMSE of EVI estimation, c mean ratio of estimates over true 99.9% quantile $\left( \nicefrac {\widehat{\chi _{0.001}}}{\chi _{0.001}}\right)$ and d RMSE of 99.9% quantile estimation for the 1500 simulations: $n=1000,\,500,\,200$ observations from a Hybrid($0.75,\,\xi$), with $\xi =-0.2,\,0.2$, quantized with $\Delta =0.01$, and sample quantile threshold candidates $\left\{ u_i\right\} _{i=1}^{I}$, $I=10,\,20$

Conclusions drawn in Sect. 3 (and also Appendix 2) regarding superiority of the ALRSM, in what concerns computational efficiency, over the competing selection methodologies, remains true when considering this quantized data scenario, both in terms of total run time and number of failed simulations. Corresponding results are summarized in Tables 10 and 11.

Table 10 Total run time, in seconds, of threshold selection and parameter estimation for 1500 simulated samples with size $n=1000,\,500,\,200$ from a Hybrid($0.75,\,\xi$), with $\xi =-0.2,\,0.2$, quantized with $\Delta =0.01$, and sample quantile threshold candidates $\left\{ u_i\right\} _{i=1}^{I}$, $I=10,\,20$; best run times for each scenario in bold

Table 11 Total failure count for 1500 simulated samples with size $n=1000,\,500,\,200$ from a Hybrid($0.75,\,\xi$), with $\xi =-0.2,\,0.2$, quantized with $\Delta =0.01$, and sample quantile threshold candidates $\left\{ u_i\right\} _{i=1}^{I}$, $I=10,\,20$

As we can infer from this small study, our selection methodology reveals itself as a good tool for choosing a single threshold even when in the presence of rounded data. However, it is our view, as mentioned in the Concluding Remarks (Sect. 5), that it can be preferable as a general principle to firstly submit the raw, discretized data to an appropriate jittering process, before further extreme value analysis—this was amply studied, for example, in the Supplementary Material of Bader et al. (2018).

1.4 Appendix 4: Mean residual life plots: significant wave height data sets

See Fig. 20.

Fig. 20

MRL plots for the hindcasts of storm peaks significant wave heights data sets, with WLS fitted line corresponding to the first minimum WMSE of the fit above each candidate (Langousis et al. 2016); dotted axis ticks at the 0.1-step used candidates, with corresponding mean excess estimates as ‘$\varvec{\times }$’; sample quantiles indicated on the upper axis scale. Left: GoM data—selected $u=3.304$; Right: NS data—selected $u=1.952$