Automatic Threshold Selection for Generalized Pareto and Pareto–Poisson Distributions in Rainfall Analysis: A Case Study Using the NOAA NCDC Daily Rainfall Database

Abstract

Both extreme-excess modeling and extreme-value analysis of precipitation events frequently utilize the Generalized Pareto (GP) distribution to model peaks above a selected threshold. However, selecting an appropriate threshold remains a complex and challenging task, which has discouraged many practitioners from employing Pareto or Pareto–Poisson distributions for extreme-value analysis. Recent analyses of threshold selection methods proposed in the technical literature, particularly when applied to rainfall records with high quantization levels, have shown that nonparametric methods are often unreliable. Additionally, methods relying on the asymptotic properties of the GP distribution tend to produce unrealistically high threshold estimates. In contrast, graphical methods and goodness-of-fit (GoF) metrics that account for the pre-asymptotic behavior of the GP distribution have demonstrated better performance. Despite these improvements, there remains no automatic and statistically robust methodology for threshold selection. This study develops an automatic, statistically sound procedure for optimal threshold selection, leveraging weighted mean square errors and internally studentized residuals. The proposed method outperforms existing approaches in terms of accuracy, as demonstrated through numerical experiments and its application to real-world data from the NOAA NCDC Daily Rainfall Database. Results indicate that the method not only improves threshold estimation precision but also enhances the reliability of extreme-value analysis for precipitation records, making it a valuable tool for hydrological applications. The findings emphasize the practical implications of the method for analyzing extreme rainfall events and its potential for broader climatological studies.

1. Introduction

Accurate estimation of extreme rainfall events is critical in hydrology for assessing flood risks and managing water resources. Two common approaches for modeling such extremes are the Annual Maxima Series (AMS) and the Partial Duration Series (PDS). The AMS approach focuses on extracting block-maxima, typically using annual maxima, while the PDS method retains values that exceed a given threshold, known as the Peaks Over Threshold (POT) method [1,2].

The choice of the threshold in PDS analysis significantly influences the accuracy of estimates derived from extreme-value theory (EVT). The Generalized Pareto Distribution (GPD) is commonly used to model exceedances over such thresholds [3,4]. However, determining an appropriate threshold remains a challenging task. Numerous methods, including graphical diagnostics and formal statistical tests, have been proposed to assist in threshold selection [5,6]. This is not surprising, as no universal solution can emerge in statistics; the choice always depends upon the specific metric or criterion being considered.

Langousis et al. [7] provide an excellent comparison of representative threshold selection methods for GPD modeling. They apply these methods to over 1700 daily rainfall records from the NOAA NCDC database, aiming to develop a robust methodology for threshold detection in extreme-value analysis. This research is particularly valuable given the variability and dependence of rainfall extremes on large-scale climatic features [8,9]. Their findings suggest that GP threshold estimates for daily rainfall range between 2 and 12 mm/d, with a mean value of around 6.5 mm/d, highlighting the effectiveness of the PDS approach in capturing extremes compared to AMS.

Despite these contributions, the work by Langousis et al. [7] fails to explicitly address the important assumption of independence between Peaks Over Threshold (POT). Ensuring that the exceedances used to fit the GPD model are independent is a key theoretical requirement in extreme-value theory. Neglecting this assumption can lead to biased parameter estimates and inaccurate results. This paper incorporates the independence hypothesis into the threshold detection process and assesses its influence on the final results, providing a more robust framework for extreme-value analysis in hydrology.

Moreover, the selection of an appropriate threshold is non-trivial. The threshold must be high enough to ensure the validity of the GPD assumptions while retaining sufficient data to avoid excessive variability in parameter estimates. Existing methods, such as the “failure-to-reject” method based on goodness-of-fit (GoF) tests are often sensitive to data quantization and lack statistical power for small sample sizes [1,10]. Additionally, methods relying on the Hill estimator and other asymptotic properties of the GPD often perform poorly under pre-asymptotic conditions [2,3], particularly for rainfall data where the shape parameter is small (0.1–0.2) [7]. In contrast, graphical methods and GoF metrics that account for pre-asymptotic behavior tend to perform better but lack automatic and statistically rigorous procedures for threshold detection [7].

The aim of this paper is twofold. First, the method proposed by Langousis et al. [7] is enhanced to ensure its application to independent peaks and to enable the automatic generation of optimal threshold values. Additionally, a new method is proposed, based on weighted least squares residuals, which incorporates two well-known techniques for bad data detection and employs the largest normalized residual test for outlier identification. This new strategy demonstrates comparable performance to the method by Langousis et al. [7] in terms of threshold estimation while offering a statistically sound decision process that can be adjusted using pre-established significance levels. The proposed method is applied for threshold selection and compared using extensive daily precipitation records from the Global Historical Climatology Network-Daily (GHCN-Daily) dataset [11,12], ensuring a robust evaluation of its effectiveness across different climatic regions.

2. Extreme-Value Analysis

In many fields, such as hydrology, statistical analysis of extremes from empirical records is carried out using either the Annual Maxima Series (AMS) or Partial Duration Series (PDS). With AMS, the block-maxima from different years are extracted, while with PDS, the peak of excesses above a sufficiently high threshold (u) is retained. This is why the approach is also known as Peaks Over Threshold (POT). It is important to note that hydrologists have been using the POT method for a long time to deal with time-dependent processes. In practice, high exceedances occur in clusters, which can represent a single storm in a meteorological context. By separating the peaks within clusters, one can assume that they are approximately independent. This allows models that require independent data, such as Pareto–Poisson, to be applied.

2.1. GEV Distribution

AM of successive years are assumed to be (i) independent random variables and (ii) identically distributed. Under these assumptions, the random variable annual maximum X follows a GEV distribution with location parameter $\mu$, scale parameter $\psi$, and shape parameter $\xi$, with a cumulative distribution function (CDF) given by:

(1)$F_X(x;\mu ,\psi ,\xi )=\left\{\begin{array}{c}exp\left\{-{\left[1+\xi \left({\displaystyle \frac{x-\mu }{\psi }}\right)\right]}_{+}^{{\displaystyle -\frac{1}{\xi }}}\right\};\xi \ne 0,\hfill \\ exp\left\{-exp\left[-\left({\displaystyle \frac{x-\mu }{\psi }}\right)\right]\right\};\xi =0,\hfill \end{array}\right.$

where ${\left[a\right]}_{+}=max(0,a)$, and the support is $x\le \mu -\psi /\xi$, if $\xi <0$, or $x\ge \mu -\psi /\xi$, if $\xi >0$. The GEV family includes three distributions corresponding to the different types of tail behavior: Gumbel ($\xi =0$) with a light tail decaying exponentially; Fréchet distribution ($\xi >0$) with a heavy tail decaying polynomially; and Weibull ($\xi <0$) with a short tail.

Remark 1.

Please note that climate change could challenge the validity of the identically distributed assumption. From a practical perspective, it is important to highlight that the researchers who laid the foundations for extreme-value analysis models likely did not believe in the strict stationary assumption. Rather, these models are often employed to design structures or systems expected to operate under conditions that represent a “snapshot in time” relative to the age of the Earth and ongoing changes. In fact, the author has worked in the past with non-stationary extreme-value analysis methods [13,14], and making robust estimates of shape parameters is already challenging in a stationary context. Extending these estimates to non-stationary models makes the task even more complex. This pseudo-stationary assumption is generally reasonable for making estimates within the expected timeframe of applicability for the analysis.

2.2. Pareto and Pareto–Poisson Distributions

The Pareto–Poisson model combines the Generalized Pareto Distribution (GPD) for the study of exceedances over a threshold u and the Poisson distribution for the occurrence of such exceedances. This model is based on the following assumptions:

1. The number of exceedances over the level u during the year has a Poisson distribution with parameter $\lambda$.

2. Those exceedances follow the GPD distribution:

   (2)$G_u(x;\psi ,\xi )=\left\{\begin{array}{c}1-{\left[1+\xi \left({\displaystyle \frac{x-u}{\psi }}\right)\right]}_{+}^{{\displaystyle -\frac{1}{\xi }}};\xi \ne 0,\hfill \\ 1-exp\left[-\left({\displaystyle \frac{x-u}{\psi }}\right)\right];\xi =0,\hfill \end{array}\right.$

   where $\sigma ,\xi$ are the scale and shape parameters, respectively, and the support is $u<x\le u+\sigma /\left|\xi \right|$ if $\xi <0$, $x\ge u$ if $\xi >0$, or $u<x<\infty$ if $\xi =0$. The Pareto distribution also includes three distribution families corresponding to the different types of tail behavior: Exponential family ($\xi =0$); Pareto tail distribution ($\xi >0$); and one analogous to the Weibull family ($\xi <0$) and a bounded tail.

Under this hypothesis, the cumulative probability distribution (CDF) of the annual maximum can be expressed as:

(3)$F_X(x;\lambda ,\psi ,\xi )=\left\{\begin{array}{c}exp\left\{-\lambda {\left[1+\xi \left({\displaystyle \frac{x-u}{\psi }}\right)\right]}_{+}^{{\displaystyle -\frac{1}{\xi }}}\right\};\xi \ne 0,\hfill \\ exp\left\{-\lambda exp\left[-\left({\displaystyle \frac{x-u}{\psi }}\right)\right]\right\};\xi =0,\hfill \end{array}\right.$

where ${\left[a\right]}_{+}=max(0,a)$, and the support is $x>u$, $x\le \psi /\left|\xi \right|$ if $\xi <0$, or $x\le \infty$ if $\xi >0$. The Pareto–Poisson family includes, as the Pareto family does, three distributions corresponding to the different types of tail behavior: Exponential ($\xi =0$); traditional Pareto tail ($\xi >0$), and an analogous Weibull distribution ($\xi <0$) with a bounded tail.

An advantage of this model is that it is entirely consistent with GEV. In fact, expression (3) is equal to (1) if the following condition holds [3]:

(4)$\begin{array}{cccc}\sigma =\psi +\xi (u-\mu );\hfill & \lambda ={\left[1+\xi \left({\displaystyle \frac{u-\mu }{\sigma }}\right)\right]}^{{\displaystyle -\frac{1}{\xi }}}\hfill & \mathrm{if}& \xi \ne 0,\hfill \\ \sigma =\psi ;\hfill & \lambda =e^{-{\displaystyle \frac{u-\mu }{\sigma }}}\hfill & \mathrm{if}& \xi =0,\hfill \end{array}$

which, in the end, allows us to work with GEV distribution for annual maxima but using data associated with the peaks over a threshold. Analytically, there is a well-known connection between the GEV distribution G and the GP distribution H, given by $H=1+ln\left(G\right)$. This equivalence ensures consistency between the two approaches and highlights the flexibility of using peaks over threshold methods for extreme-value analysis.

Please note that although this manuscript focuses on threshold selection, it is essential to validate the assumptions underlying the fitted models once the threshold is determined. Specifically, the Poisson assumption for exceedances and the overall goodness of fit of the GEV or Pareto–Poisson models should be rigorously tested using the real-world sample.

From a practical perspective, it is recommended to use the Pareto–Poisson (PP) model because (a) it takes into account more than just a single data point (e.g., annual maxima) per year, as it considers all exceedances, and (b) it excludes years where the annual maximum value is not representative of the upper tail of the distribution. For instance, in years with no extreme precipitation events, the annual maxima method may select a relatively low value that does not belong to the upper tail. Such values can distort the analysis, as they do not reflect the true behavior of extreme events. Additionally, contrary to the belief that POT values converge to the Pareto distribution as the threshold approaches infinity [15], they actually converge when the threshold approaches the upper end of the random variable [16]. However, there are also drawbacks, the most important being the need for proper selection of the threshold u to obtain reliable estimates, as pointed out by [7].

It is important to note that the expectation and variance of the excesses above a threshold may not exist for certain values of the shape parameter $\xi$. Specifically, the variance is undefined when $\xi \ge 0.5$, and the expectation does not exist for $\xi \ge 1$. This limitation is also present in the GEV approach and highlights the importance of carefully considering the shape parameter when applying extreme-value models. However, in the context of precipitation, this issue is typically not encountered, as the shape parameter is generally estimated to be within the range of $0.1$ to $0.2$ [7], far from the critical values where the expectation or variance would be undefined. This ensures that the models remain mathematically valid and practically applicable for precipitation data.

2.3. Mean Residual Life Plot

The most popular graphical method for selecting a threshold is based on the mean excess or mean residual life plot (MRLP) [1,3]. Assuming that the conditional excesses $[X-u^{*}|X>u^{*}]$ of a random variable X above threshold $u^{*}$ follow a GP distribution with scale and shape parameters $\sigma$ and $\xi$, respectively, then for any threshold $u>u^{*}$, the conditional excesses $[X-u|X>u]$ also follow a GP distribution with the same shape parameter $\xi$ and a scale parameter that depends linearly on the threshold u:

(5)${\sigma }_u=\sigma +\xi (u-u^{*}).$

The expectation of the excesses corresponds to:

(6)${\displaystyle E[X-u|X>u]=\frac{{\sigma }_u}{1-\xi }=\frac{\sigma +\xi (u-u^{*})}{1-\xi }={\beta }_{1}u+{\beta }_0,}$

where the slope and the intercept of the linear relationship can be obtained as a function of the GP distribution parameters:

(7)$\begin{array}{ccc}\hfill {\beta }_0& =& {\displaystyle \frac{\sigma -\xi u^{*}}{1-\xi }}\hfill \end{array}$

(8)$\begin{array}{ccc}\hfill {\beta }_1& =& {\displaystyle \frac{\xi }{1-\xi },}\hfill \end{array}$

while the variance of excesses is equal to:

(9)${\displaystyle V[X-u|X>u]=\frac{{\sigma }_u^2}{{(1-\xi )}^2(1-2\xi )}=\frac{{\left(\sigma +\xi (u-u^{*})\right)}^2}{{(1-\xi )}^2(1-2\xi )}.}$

As stated by [7] in their research, this method for selecting thresholds provides better benefits when compared to other methods. However, there are two main concerns with this method. The first concern is related to automation, as it requires visual inspection. The second concern is subjective judgment, which could lead to significant bias if the selection is not appropriate. In an attempt to automate threshold selection, ref. [7] suggested an iterative procedure given a sample $x_1,x_2,x_3,\dots ,x_n$ of size n:

1. Sort the sample in ascending order, resulting in $x_{1:n},x_{2:n},x_{3:n},\dots ,x_{4:n}$ where $x_{i:n}$ corresponds to the i-th order statistic.

2. For each order statistic $x_{i:n}$ from $i=1$ to $i=n-10$, set the threshold equal to the corresponding order statistic, i.e., $u_i=x_{i:n}$, and compute the sample mean of excesses ${\overline{e}}_i$ above the threshold $u_i$ and the sample standard deviation $s_i$. Note that a minimum of 10 excesses is ensured to compute the conditional mean and standard deviation in the worst-case scenario.

3. For each threshold $u_i;i=1,2,\dots ,n-20$, use the weighted least squares (WLS) method to fit a linear regression model to points $(u_j,{\overline{e}}_j);\forall j\ge i$. To account for the change in variance with increasing threshold according to (9) and the sample size variability, the weight $w_j$ associated with each point corresponds inverse of the variance, i.e., $w_j=(n-j)/s_j^2$.

4. The optimal threshold selection is determined by finding the lowest threshold that corresponds to a local minimum of the weighted mean square error (WMSE) function in linear regression.

Although the method proposed by [7] allows for some level of automation and produces acceptable results for the stations they analyzed, the authors acknowledge that the rule cannot be considered unique because it lacks statistical support.

3. Weighted Least Squares (WLS)

Consider the standard linear regression model

(10)$\mathit{Y}=\mathit{X}\mathit{\beta }+\mathit{\epsilon },$

where $\mathit{Y}={(Y_1,Y_2,\dots ,Y_n)}^T$ is a $n\times 1$ response variable vector, X is a $n\times k$ matrix of predictor variables often called “design matrix”, $\mathit{\beta }$ is a $k\times 1$ vector of regression coefficients or parameters, and $\mathit{\epsilon }={({\epsilon }_1,{\epsilon }_2,\dots ,{\epsilon }_n)}^T$ is a $n\times 1$ vector of random errors assumed to be jointly normally distributed random variables $\mathit{\epsilon }\sim N(\mathbf{0},\mathit{V})$, where V is a positive definite diagonal variance-covariance matrix.

Regression parameters $\mathit{\beta }$ are usually estimated using the WLS method,

(11)$\underset{\mathit{\beta }}{\mathrm{Minimize}}{\mathit{\epsilon }}^T\mathit{W}\mathit{\epsilon }=\underset{\mathit{\beta }}{\mathrm{Minimize}}{\left(\mathit{Y}-\mathit{X}\mathit{\beta }\right)}^T\mathit{W}\left(\mathit{Y}-\mathit{X}\mathit{\beta }\right),$

where $\mathit{W}={\mathit{V}}^{-1}$. For model (10), WLS coincides with the maximum likelihood (ML) estimation method. Please note that for homoscedastic models, W corresponds to the identity matrix, i.e., $w_{ii}=1/{\sigma }^2;i=1,\dots ,n$; $w_{ij}=0;i,j=1,\dots ,n$ and $i\ne j$, and (11) becomes the traditional least squares (LS) method. However, matrix W is included in the formulation so that regression formulas remain valid for the reweighting approach presented by [7]. Fitting results are [17]:

(12)$\begin{array}{ccc}\hfill \widehat{\mathit{\beta }}& =& {\left({\mathit{X}}^T\mathit{W}\mathit{X}\right)}^{-1}{\mathit{X}}^T\mathit{W}\mathit{Y}\hfill \end{array}$

(13)$\begin{array}{ccc}\hfill \widehat{\mathit{Y}}& =& \mathit{X}\widehat{\mathit{\beta }}=\mathit{P}\mathit{W}\mathit{Y}\hfill \end{array}$

where the hat (^) refers to estimates, and

(14)$\begin{array}{ccc}\hfill \mathit{P}& =& \mathit{X}{\left({\mathit{X}}^T\mathit{W}\mathit{X}\right)}^{-1}{\mathit{X}}^T\hfill \end{array}$

(15)$\begin{array}{ccc}\hfill \widehat{\mathit{\epsilon }}=\mathit{Y}-\widehat{\mathit{Y}}& =& (\mathit{I}-\mathit{P}\mathit{W})\mathit{Y}\hfill \end{array}$

where P is the projection matrix. The residual mean square estimator of ${\sigma }^2$ is

(16)${\displaystyle {\widehat{\sigma }}^2=\frac{{\mathit{\epsilon }}^T\mathit{W}\mathit{\epsilon }}{n-k}.}$

3.1. Differences Between Influential Observations and Outliers

Influential observations can be defined, according to [18], as those observations with larger and excessive impact on the calculated values of some estimates. There are numerous influence measures in the literature, which, according to [19], can be classified into five groups based on: (1) residuals, (2) the prediction matrix, (3) volume of confidence ellipsoids, (4) influence functions, and (5) partial influence. Outliers, in contrast, are data points that cannot be explained by the model, as they are generated under dynamics distinct from those of regular data. A more intuitive definition describes them as observations that deviate from or disrupt an established pattern. One can find outliers that are influential as well as outliers that are not. Some outliers present large residuals and, therefore, are easy to detect. However, it is important to realize that some outliers may have small residuals because they have a large influence on the parameter estimates; when outliers of this type appear in groups, they are often more difficult to detect even though they are very influential. Finally, there may be some outliers with small residuals that are not influential; these are also difficult to detect, but they are much less important.

3.2. Influence Measures

To assess the effect of outliers’ different influence measures can be used. Among those, internally studentized residuals are widely used for outlier detection, which are a scaled version of residuals, that is

(17)${\displaystyle z_i=\frac{\sqrt{w_{ii}}{\epsilon }_i}{\widehat{\sigma }\sqrt{1-w_{ii}{p}_{ii}}};i=1,\dots ,n.}$

For “outlier” identification purposes, an internally studentized residual corresponds to a suspected “bad” data with a $1-\alpha$ confidence level (e.g., $0.99$) if $\left|z_i\right|>{\Phi }^{-1}(1-\alpha /2)$.

Since the aim is to identify the minimum threshold that allows for an acceptable fit to the Pareto distribution or ensures that the MRLP aligns with a straight line, an appropriate threshold would result in the studentized residual of the first observation not being excessively large, indicating that it has limited influence. This ensures that the fit is not unduly impacted. If the threshold is smaller than optimal and the left-tail data exhibits a very high studentized residual, it would suggest that this point has a large influence on the regression and may not belong to the Pareto tail. This forms one of the key criteria for our threshold selection process. Large values of studentized residuals are typically related to outliers in the response-factor space, representing points not well fitted by the model. This is particularly concerning if the largest studentized residual corresponds to the minimum exceedance above the selected threshold, signaling that the threshold may be inappropriate. Outliers, especially in the left tail of the distribution, are expected to exert significant influence, and this can be effectively detected using the studentized residual as a diagnostic measure.

This method has also been used in the literature to identify outliers in wave hindcast databases related to hurricanes and typhoons, which are often misrepresented due to modeling limitations [20]. By detecting and managing these extreme events, this approach improves the accuracy of statistical characterizations for offshore and coastal design applications.

3.3. Bad Data Detection

If errors are Gaussian-distributed unbiased, i.e., centered with zero mean, and independent, weighted square sum of residuals ${\widehat{\mathit{\epsilon }}}^T\mathit{W}\widehat{\mathit{\epsilon }}$ follows a ${\chi }^2$ distribution with $n-2$ degrees of freedom [21]. Thus,

(18)$\mathrm{Prob}\left({\widehat{\mathit{\epsilon }}}^T\mathit{W}\widehat{\mathit{\epsilon }}\le {\chi }^2(1-\alpha ,n-2)\right)=1-\alpha ,$

where $\alpha$ is the significance level.

Therefore, for a given $\alpha$ (e.g., $0.01$), the value ${\chi }^2(1-\alpha ,n-2)$ can be computed and the ${\chi }^2$ test applied, i.e., if ${\widehat{\mathit{\epsilon }}}^T\mathit{W}\widehat{\mathit{\epsilon }}<{\chi }^2(1-\alpha ,n-2)$ there is no evidence against the null hypothesis of the non-existence of bad data (improper threshold selection) at the $1-\alpha$ confidence level, otherwise, there is. Further details can be found in [22].

In summary, the threshold is sequentially searched, starting with small values, until both the chi-square test and the outlier detection test for the minimum exceedance are satisfied.

4. Independent Peak Detection and Threshold Selection

In extreme-value analysis, particularly when applying the Peaks Over Threshold (POT) method, it is essential to ensure that the peaks considered are independent of one another. The assumption of independence is a key theoretical requirement for the validity of extreme-value theory (EVT), which relies on the assumption that exceedances over a given threshold belong to the Generalized Pareto Distribution (GPD) [1]. Dependent peaks can introduce biases in parameter estimation, leading to erroneous conclusions about the distribution of extremes. To address this issue, a minimum peak distance is imposed to prevent the selection of multiple peaks originating from a single event (e.g., a storm) or closely related phenomena. This approach enhances the likelihood that the exceedances used to fit the GPD model are independent, therefore improving the reliability of the statistical inference. Furthermore, a cross-correlation test can be conducted for lags between sequential peaks, which are not necessarily evenly distributed in time, using methods such as the Ljung–Box test. If statistically significant evidence confirms the absence of correlations, it provides additional support for the assumption of independence. However, it is important to note that the absence of correlation does not necessarily imply independence.

Moreover, when dealing with highly quantized datasets—where data are recorded with limited precision or in discrete intervals—it is essential to use unique peak values to reduce computational load. Quantization often results in multiple identical peaks, which should be considered in estimating exceedance statistics. This approach is particularly important when the recording instrument’s resolution is coarse relative to the data’s variation scale [4].

The process of extracting independent peaks involves identifying exceedances above a specified threshold while ensuring that the peaks are separated by a minimum distance to preserve their independence. This assumption is further reinforced by performing a correlation test at different time lags. If the assumption of independence is validated by the correlation test (which is not conclusive), the mean exceedance, variance, and weight are computed for each unique peak based on the number of exceedances above the peak and their variability. If the number of exceedances is insufficient (fewer than a predefined threshold, $n_0$), the process is halted to avoid introducing unreliable estimates. The computed weights are then used to assign greater importance to peaks with higher variability in the subsequent threshold selection process.

The algorithm works as follows:

1. Detect all peaks in the data that exceed the given threshold and enforce a minimum peak distance to ensure independence.

2. Sort the unique peaks in ascending order.

3. For each unique peak, calculate the mean exceedance, variance, and weight based on the ratio between the number of exceedances and the variance.

4. If the number of exceedances is less than the required minimum $n_0$, truncate the data to retain only valid exceedances. $n_0$ was chosen to be 10 by [7]. Although this serves as a minimum value, in practice, the data quantization means this limit is effectively nonrestrictive: each unique peak typically has more than 10 exceedances, except for the very last one, which can be skipped. Importantly, this specific value is not critical to the method since the primary focus lies on accurately modeling the left tail of the exceedance distribution. Additionally, the data weighting is applied proportionally to the sample size to account for variability in the number of exceedances.

5. Remove any invalid data (e.g., NaN values) to ensure that the subsequent analysis is robust and free from errors.

6. Perform a correlation test for a sequential lag of 1, using p-values derived from Pearson’s correlation with a Student’s t distribution applied to the transformed correlation coefficient or a Ljung–Box Test. If the correlation is found to be significantly different from zero, increase the minimum peak distance and return to step 1 to re-evaluate the process.

This procedure enhances the likelihood of meeting the independence assumption while also addressing challenges related to data quantization by focusing on unique values. As a result, it facilitates a more accurate and reliable threshold selection process. Upon completing this procedure, the time series of independent peaks, unique peaks, and their corresponding mean exceedances, variances, and weights are obtained, providing the essential inputs for applying various threshold selection methods.

In Figure 1, the method described above is applied to positive rain rates extracted from a 134-year record of daily rainfall observations from Australia, obtained from the NOAA NCDC rainfall database, also utilized by [7]. The full precipitation record is shown in gray, independent peaks are marked with black triangles, and the MRLP is plotted using the unique peaks and their corresponding exceedances. The correlations and p-values (in parentheses) for the first three sequential lags are $0.0117\left(0.4079\right)$, $0.0228\left(0.1079\right)$, and $0.0222\left(0.1165\right)$, respectively. These results demonstrate that the linear correlations are very small and not significantly different from zero, with p-values exceeding typical significance levels of 1%, 5%, or 10%, therefore supporting the independence hypothesis. In such cases, it is also helpful to create a scatter plot of peaks with different time lags to assess whether the linearity assumption holds, which is necessary for linking correlation and independence.

5. Updated Langousis Algorithm for Threshold Selection

The method proposed by [7] for threshold selection has been improved by incorporating an independence check and applying a cubic smoothing spline to identify local minima in the mean square error (MSE) function. The choice is motivated by its desirable mathematical properties, including the existence and continuity of the first two derivatives, which ensure smoothness and stability in the approximation. These local minima are then used as candidate thresholds for selection.

This algorithm works in conjunction with the independent peak detection process outlined in the “Independent Peak Detection and Threshold Selection”section. The independent peaks, along with the mean exceedances and variance for the data points exceeding those peaks—obtained using the Peaks Over Threshold (POT) method—form the basis for evaluating candidate thresholds.

5.1. Algorithm Steps

The updated Langousis algorithm proceeds as follows:

1. From the unique independent peaks identified during the independent peak detection and threshold selection process, sort the peaks in ascending order. These sorted peaks form the potential thresholds for the Langousis method.

2. For each potential threshold $u_i$, and the mean exceedances and weight for the data points exceeding this threshold. Use the weighted least squares (WLS) method to fit a linear regression model to the exceedances and calculate the corresponding mean square error (MSE).

3. Apply a cubic smoothing spline to the logarithm of the MSE values across the thresholds. The smoothing parameter p controls the degree of smoothness in the spline and can be adjusted to prevent overfitting. The default value of $p=0.9$ has been found to be effective for daily precipitation records.

4. Compute the first derivative of the spline to find the points where the derivative is zero, identifying potential local minima and maxima.

5. Evaluate the second derivative of the spline at these points to classify them. Local minima are detected where the second derivative is positive, and these local minima are selected as candidate thresholds.

6. All local minima, if available, are retained as the final candidate thresholds for further analysis.

5.2. Mathematical Basis

The mathematical foundation of the method lies in the WLS model used to fit the exceedance data. For each threshold $u_i$, the weights are computed based on the sample variance of the exceedances. This ensures that the selected thresholds minimize variability while retaining enough data points to fit the Generalized Pareto Distribution (GPD) effectively. The spline fitting smooths the MSE curve, and the detection of local minima provides a rigorous way of identifying potential thresholds. Please note that the selection of the smoothing parameter is intended to ensure that the spline closely follows the overall curve, avoiding spurious local minima in the MSE function for each peak value caused by data noise. The value $p=0.9$ has been tested and found appropriate for daily precipitation records; however, its suitability should be confirmed if the method is applied to other types of data. Figure 2 illustrates how the smoothing curve effectively follows the data, providing the desired type of fit when selecting the smoothing parameter.

5.3. Output of the Algorithm

The algorithm returns all detected local minima as candidates for threshold selection. These local minima provide a range of possible threshold values, allowing the user to choose the most appropriate one based on the characteristics of the data. While the first local minimum is often used as the default threshold, additional local minima can be evaluated to ensure robustness in the analysis.

5.4. Graphical Representation

An optional graphical output is provided to visualize the Langousis MSE data, the fitted spline, and the detected local minima. As shown in Figure 2, the local minima are marked with red triangles, and vertical dashed lines indicate their positions. This graphical representation aids in understanding the behavior of the MSE function across different thresholds.

5.5. Advantages of the Updated Algorithm

The updated version of the Langousis method provides the following key benefits:

1. Automation: The inclusion of the spline fitting and automatic detection of local minima eliminates the need for manual intervention, resulting in a more objective and consistent threshold selection process.

2. Statistical Soundness: Using weighted least squares combined with local minima detection, the algorithm provides a statistically sound basis for threshold selection. The use of multiple local minima ensures robustness in the threshold selection process, especially in cases where the MSE function is complex.

6. Automatic Threshold Selection Procedure

The threshold selection algorithm based on studentized residuals builds on the method proposed by [7] but incorporates the use of independent peaks and the outlier detection and identification techniques explained in Section 3. Given an initial sample, the updated algorithm works as follows:

• (Initialization): Select a significance level $\alpha$ (possible alternatives are $\alpha =0.01$, $\alpha =0.05$, or $\alpha =0.1$) and choose an initial minimum threshold $u_{\min }$ to start the process. A minimum peak distance to ensure independence between the peaks is defined.

• Independent Peak Detection): Extract independent peaks over the threshold $u_{\min }$ using the Peaks Over Threshold (POT) method and the minimum peak distance to enforce independence. This step is detailed in the “Independent Peak Detection and Threshold Selection” section.

• (Sorting): Sort unique peaks in ascending order, resulting in $x_{1:n},x_{2:n},\dots ,x_{n:n}$, where $x_{i:n}$ corresponds to the i-th order statistic.

• (Mean Exceedance and Variance Computation): For each order statistic $x_{i:n}$, set the threshold equal to the corresponding order statistic, i.e., $u_i=x_{i:n}$, and compute the sample mean of exceedances ${\overline{e}}_i$ above the threshold $u_i$ as well as the sample standard deviation $s_i$. Ensure that at least 10 exceedances are used to compute the conditional mean and standard deviation, even in the worst-case scenario. The choice of 10 exceedances, while seemingly arbitrary, ensures a minimum sample size for the stable computation of variance. Importantly, this specific value is not critical to the method since the primary focus lies on accurately modeling the left tail of the exceedance distribution. Additionally, the data weighting is applied proportionally to the sample size to account for variability in the number of exceedances.

• (Weighted Least Squares Fit): For each threshold $u_i$ starting from $u_{\min }$, apply the weighted least squares (WLS) method to fit a linear regression model to the points $(u_j,{\overline{e}}_j)$ for all $u_j\ge u_i$. The weights $w_j$ for each point are computed as the inverse of the variance, $w_j=n_j/s_j^2$, to account for the increasing variance with higher thresholds, where $n_j$ is the number of independent peaks above the threshold used to compute the mean.

• (Studentized Residuals Computation): Compute the internally studentized residuals $r_j$ for each threshold using the residuals of the WLS fit using (17). These residuals help identify outliers and ensure that the first data point (smallest exceedance) does not unduly influence the regression.

• (Stopping Criteria): Evaluate the goodness of fit using the weighted sum of squared residuals ${\widehat{\epsilon }}^{T}W\widehat{\epsilon }$. If the value is less than the critical value from a ${\chi }^2$ distribution with $n-2$ degrees of freedom, proceed to test the studentized residuals for outliers. If the studentized residual associated with the minimum exceedance (left tail of the exceedance distribution) is within acceptable bounds and does not indicate an outlier (e.g., $|r_1|<z_{1-\alpha /2}$, where $z_{1-\alpha /2}$ is the critical value for the normal distribution at the chosen significance level), the current threshold is selected as optimal.

• (New Threshold Selection): If the sum of squared residuals exceeds the critical value or if there is evidence of an outlier (e.g., $|r_1|>z_{1-\alpha /2}$), select a new candidate threshold. The new candidate is chosen as the closest local maximum of the studentized residual absolute values. This ensures that other possible thresholds that might better balance the fit and the residuals are explored.

• (Iterate and Finalize): Repeat the above steps iteratively, adjusting the threshold until both the Chi-square goodness-of-fit test and the studentized residuals test are satisfied. Once a threshold is selected, stop the procedure and report the optimal threshold $u_{\mathrm{opt}}$.

This procedure integrates both the detection of independent peaks and the use of residuals to ensure the robustness of the selected threshold.

7. Synthetic Numerical Experiment

To evaluate the performance of various threshold selection methods in extreme-value analysis, a synthetic numerical experiment was designed. This experiment aimed to assess the performance of the different threshold selection methods under controlled conditions with known characteristics. The objective was not to replicate the real behavior of precipitation records but rather to compare their performance when the optimal threshold is known. The experiment is based on the following assumptions:

1. The record length mimics daily precipitation records over 40 years (14,610 data points), simulating long-term environmental datasets.

2. The generated data combines Gaussian-distributed values with a Generalized Pareto Distribution (GPD) to model the heavy-tailed behavior of extremes. Please note that only the tail behavior approximates extreme precipitation patterns.

3. The central part of the distribution follows a normal distribution, which does not represent precipitation records; however, negative values were removed and replaced with zeros to mimic the effect of the absence of precipitation.

4. The threshold was established based on a quantile criterion. For values exceeding a specified quantile over the normally distributed sample, data were sampled from a Pareto distribution with varying scale and shape parameters, ensuring that the tail was governed by the Pareto distribution.

5. The data sampled is statistically independent, which aligns with the assumption typically made when selecting peaks in precipitation records.

It is important to emphasize that the goal of this experiment is not to replicate the dynamics of precipitation patterns. Instead, the experiment is designed to evaluate the performance of the methods in scenarios with some similarities, where the true solution is known and can serve as a benchmark.

The experiment systematically varied several parameters to assess the robustness of the methods under different conditions:

• Threshold quantiles (CDF thresholds): $\{0.99,0.995,0.999\}$, representing the proportion of data exceeding the threshold.

• GPD shape parameters ($\xi$): $\{-0.15,-0.1,-0.05,0,0.05,0.1,0.15\}$, controlling the tail heaviness of the extreme values.

• Scale parameters ($\sigma$): $\{0.5,1,2\}$, influencing the variability of the extreme values.

• Rounding levels: $\{1,2,3\}$, mimicking the effects of data quantization.

• Significance levels ($\alpha$): $\{0.01,0.05\}$, defining thresholds for the statistical tests applied.

• Runs per combination: 10 independent simulations to assess variability in the results.

The rest of the features of the experiment are:

• Data Generation: Gaussian-distributed daily precipitation data were generated with truncation at zero to simulate non-negative rainfall. For values exceeding the specified CDF threshold, a GPD transformation was applied using the given shape ($\xi$) and scale ($\sigma$) parameters.

• Peak Extraction: Unique peaks were identified with a minimum separation of one day to ensure independence. Since the data were synthetically generated, the independence assumption holds in practice. For each peak, exceedance statistics (mean, variance, and weights) were calculated, provided that at least 10 exceedances were available to ensure robust estimation.

• Threshold Selection Methods: The following four methods were applied to determine thresholds:

  1. Studentized Residuals (SR): Identifies thresholds based on statistical residuals.

  2. Langousis Method (L): Evaluates thresholds using local minima of a smoothed error curve.

  3. Anderson–Darling (AD): Selects thresholds based on goodness-of-fit statistics for the GPD.

  4. Cramer–Von Mises (CVM): Similar to AD, but uses a different fit statistic.

For each parameter combination and each method, the results included (i) the real threshold values, (ii) selected thresholds, (iii) computational times, (iv) errors between estimated and real thresholds, and (v) the percentage of cases that the absolute difference between real and estimated thresholds is below 1 unit.

The performance of different threshold selection methods is compared in Figure 3, which shows their median thresholds, variability, and robustness across significance levels. The Langousis methods (L1, L2, L3) demonstrate varying levels of reliability, with L1 showing the lowest median and narrower variability, but 8.9% of the results are NaN, and only 57.6% of the deviations are within an absolute error of 1. L2 achieves the highest percentage of deviations within 1 (87.1%) while maintaining a moderate NaN value of 28.6%, making it a balanced choice. However, L3 exhibits the highest variability and NaN proportion (54.4%), reducing its reliability. The Studentized Residuals (SR) methods perform consistently, with no NaN values and percentages of deviations within 1 of 56.9% (SR 1%) and 63.3% (SR 5%), indicating better performance at the 5% significance level. The Anderson–Darling (AD) and Cramer–Von Mises (CVM) methods, while robust (no NaN values), show lower percentages of deviations within 1, particularly AD 1% (39.3%) and CVM 1% (35.0%). Overall, SR 5% and L2 provide the most accurate thresholds, while CVM and AD methods appear less reliable in terms of precision. This highlights the trade-offs between numerical robustness and accuracy among the methods.

Comparison of Threshold Selection Methods Across Rounding Levels

To evaluate the performance of the threshold selection methods under different levels of rounding precision, three rounding levels ($r=1$, $r=2$, and $r=3$) were tested. The results are summarized in Table 1, where the percentage of NaN values and the percentage of deviations within 1 are presented for each combination of method and rounding level. These metrics provide insight into the robustness and accuracy of the methods under varying degrees of data quantization.

For rounding level $r=1$, the Langousis methods (L1, L2, and L3) display higher percentages of NaN values, particularly for L3 (80.3%). In contrast, SR 1%, SR 5%, AD 1%, AD 5%, CVM 1%, and CVM 5% methods exhibit no NaN values, demonstrating their robustness to coarse rounding. Regarding the percentage of deviations within 1 (<1%), the performance of L3 is notably low (4.8%), while SR 1% and SR 5% achieve higher percentages of 58.9% and 72.5%, respectively, followed by AD 1% (46.8%) and CVM 5% (51.6%). For rounding level $r=2$, the performance improves for the Langousis methods, with reduced NaN values. Specifically, L1, L2, and L3 show NaN percentages of 2.5%, 23.7%, and 55.9%, respectively. The percentage of deviations within 1 also increases significantly, with L1 achieving 62.2% and L2 72.2%. This indicates enhanced robustness for Langousis methods under moderate rounding. Similarly, the SR, AD, and CVM methods maintain consistent performance, with SR 1% achieving 58.4% and SR 5% 63.8%, while AD 5% and CVM 5% perform at 44.3% and 42.9%, respectively. For rounding level $r=3$, the percentage of NaN values is negligible for all methods, with L1, L2, and L3 reporting NaN percentages of 0.2%, 2.9%, and 27.0%, respectively. The percentage of deviations within 1 reaches its highest for L2 (78.6%), highlighting its improved stability at finer rounding levels. SR 5% remains consistent with 53.5%, while AD and CVM methods perform steadily, with AD 1% achieving 34.0% and CVM 5% at 39.8%. Overall, the results suggest that the methods become more stable under increased rounding precision. Langousis methods, in particular, benefit from finer rounding, exhibiting a notable reduction in NaN percentages and improved accuracy for deviations within 1. SR and AD methods consistently perform well across all rounding levels, showcasing their robustness to data quantization effects.

If the effect of the Pareto shape parameter is considered, the analysis of Table 2 reveals significant insights into the performance of threshold selection methods across different shape parameter values ($\xi$). Regarding the percentage of NaN values (NaN %), the Langousis methods (L1, L2, and L3) exhibit small dependence on $\xi$. The NaN percentages are highest for L3, with values reaching up to 57.8% for $\xi =-0.15$, indicating the small sensitivity of this method to tail behavior. Conversely, SR 1%, SR 5%, AD 1%, AD 5%, CVM 1%, and CVM 5% consistently show no NaN values, demonstrating their robustness to varying $\xi$ values. The percentage of deviations within 1 (<1%) highlights different trends. L2 achieves the highest performance among the Langousis methods, particularly for $\xi =-0.15$ (67.8%) and $\xi =0.15$ (59.3%), indicating its relative stability. Among the other methods, SR 5% consistently provides the highest values, with percentages exceeding 60% for most $\xi$ values and peaking at 65.9% for $\xi =0.00$. AD and CVM methods show moderate performance, with the percentages generally improving for lower values of $\xi$. The CVM 5% method achieves notable values of 53.7% for $\xi =-0.10$ and $\xi =-0.15$. Overall, the results indicate that SR 5% is the most robust across all $\xi$ levels, followed by the Langousis L2 method, which shows improved performance for moderate values of $\xi$. The AD and CVM methods also demonstrate consistent behavior, with slight variations depending on the tail parameter. This analysis underscores the importance of selecting methods that balance robustness and accuracy under different tail behaviors.

Finally, the computational time comparison across methods, shown in Figure 4, reveals notable differences in execution times. The Langousis method demonstrates significantly higher computational time compared to Studentized Residuals, Anderson–Darling, and Cramer–Von Mises. However, these differences are unlikely to impact practical applications on modern computational systems, where such time differences are negligible given the high performance of current processors.

In conclusion, for synthetic records, particularly those with highly quantized data, the Studentized Residuals method with a 5% significance level consistently exhibits the best overall performance, combining accuracy and reliability. While the Langousis method performs exceptionally well in certain scenarios, it also displays robustness issues, making it less consistent under varying conditions.

8. Illustrative Example

This section demonstrates the application of the automatic threshold selection method using studentized residuals on positive rainfall rates derived from a 134-year daily rainfall record from Australia, as shown in Figure 2. This dataset, also analyzed by [7], led to an optimal threshold choice of 2 mm/day.

The algorithm proceeds as follows when applied to the Australian rainfall data:

• Initialization: Set significance levels $\alpha =0.05$ and $\alpha =0.01$ (test both significant levels), and initialize with a minimum threshold $u_{\min }=0.0$ (the sample minimum value). The minimum peak separation is set to 2 days to enhance the likelihood of independence between peaks. The choice of a minimum peak distance of 2 days is based on practical hydrological considerations. For precipitation data, this threshold ensures that peaks are sufficiently separated to reduce the likelihood of dependence between consecutive exceedances caused by the same weather event. While this value is not universal, it is well suited for daily rainfall data with typical storm durations of 1–2 days. Researchers working with other types of data should adjust the minimum peak distance based on the expected temporal correlation structure of the phenomena under study.

• Independent Peak Detection: The initial sample includes 46,088 observations, with 9836 values exceeding $u_{\min }=0.0$. Applying the 2-day minimum separation, the number of potentially independent peaks reduces to 4987, of which 295 are unique values. This reduction is induced by the measurement precision and is closely related to the quantization level of the data. This sample forms the basis for $x_1,x_2,x_3,\dots ,x_n$ with $n=295$. The lag-1 correlation and p-value (in parentheses) are $0.0117\left(0.4079\right)$, which supports the adequacy of the 2-day minimum peak separation criterion.

• Sorting: Sort the unique peaks in ascending order, yielding $x_{1:n},x_{2:n},\dots ,x_{n:n}$, where $x_{i:n}$ denotes the i-th order statistic. Figure 2 presents the final list of independent peaks.

• Mean Exceedance and Variance Calculation: For each order statistic $x_{i:n}$, set the threshold to $u_i=x_{i:n}$, then compute the mean exceedances ${\overline{e}}_i$ and standard deviation $s_i$. Make sure that the number of elements to compute the mean and standard deviation is at least 10. Calculations are based on the 4987 peak records.

• Regression Fit: Perform the first regression fit with $u_i=0.1$, computing the internally studentized residuals shown in Figure 5a. The residual sum of squares is ${\widehat{\epsilon }}^{T}W\widehat{\epsilon }=115.4919$.

• Stopping Criteria: The $0.95$ and $0.99$ quantiles of the ${\chi }^2$ distribution with 293 degrees of freedom are $333.9219$ and $352.2375$, respectively, both larger than ${\widehat{\epsilon }}^{T}W\widehat{\epsilon }=115.4919$. Thus, accept the null hypothesis of no outliers. However, the first internally studentized residual absolute value, $|{\widehat{r}}_1|=5.1744$, exceeds the $z_{1-\alpha /2}$ quantiles ($1.96$ and $2.576$), indicating a high-impact outlier in the lower tail. Therefore, the hypothesis that no outliers are present is rejected, i.e., the current threshold does not belong to the upper tail of the distribution.

• New Threshold Selection: The next threshold is selected based on the largest absolute value among local maxima/minima of the normalized residual sensitivities above $u_{\min }$. This yields a new threshold of $u_{\min }=2$, as shown in Figure 5b.

• Further Regression Fit and Stopping Criteria: Perform a second regression fit with $u_i=2.0$. The sum of squares residuals is ${\widehat{\epsilon }}^{T}W\widehat{\epsilon }=62.0728$. Since this value is below both ${\chi }^2$ quantiles, accept the null hypothesis. However, ${\widehat{r}}_1=2.3115$ surpasses the $z_{0.95}=1.96$ quantile, so:

  1. For $\alpha =0.01$, the optimal threshold selection is ${\widehat{u}}_{\min }=2$ (${\widehat{r}}_1=2.3115<2.576$).

  2. For $\alpha =0.05$, outliers are indicated, prompting another threshold selection.

• Final Threshold Selection: The threshold yielding the largest absolute residual among local maxima/minima leads to $u_{\min }=10.6$ (see Figure 6b).

• Stopping Criteria for Final Fit: With ${\widehat{\epsilon }}^{T}W\widehat{\epsilon }=41.7095$ and first residual ${\widehat{r}}_1=0.7473$ well below the quantiles, the optimal threshold is ${\widehat{u}}_{\min }=10.6$ for $\alpha =0.05$ as shown in Figure 7.

Alternative threshold estimates from the Anderson–Darling and Cramer–Von Mises methods are $17.3$ and $12.7$, respectively. These findings align with [7] and indicate that using the asymptotic distributions of $A^2$ and $W^2$ in the failure-to-reject method may overestimate the threshold.

For model estimates, these are the results:

1. For $\alpha =0.05$: With $u^{*}=10.6$, ${\widehat{\beta }}_0=6.4603$, ${\widehat{\beta }}_1=0.2610$, yielding $\widehat{\sigma }=7.3171$, and $\widehat{\xi }=0.2070$.

2. For $\alpha =0.01$: With $u^{*}=2$, ${\widehat{\beta }}_0=6.8592$, ${\widehat{\beta }}_1=0.2377$, yielding $\widehat{\sigma }=5.9262$, and $\widehat{\xi }=0.192$.

These estimates closely match those in [7] with slight variations, implying that independence assumptions affect the results. Dependent peaks could reduce expected excesses for small thresholds but converge to independence at high thresholds, increasing the mean excess plot slope and, thus, the shape parameter estimates.

Figure 2 shows the modified Langousis’ method results, indicating similar threshold choices. Different local optima correspond to varying significance levels, allowing the Langousis’ method to be tuned accordingly.

9. Application of Threshold Selection Methods to the GHCN-Daily Dataset

In this section, several threshold selection methods are applied to the Global Historical Climatology Network-Daily (GHCN-Daily) dataset (Version 3.31). The GHCN-Daily dataset, which is a comprehensive record of daily climatic variables from multiple locations worldwide, is maintained by the NOAA National Climatic Data Center. For a detailed background on the dataset, reference can be made to [11,12].

9.1. Dataset Description

The GHCN-Daily dataset (Version 3.31) includes daily climate records, including rainfall data, from a wide array of global stations spanning over 100 years. Each station’s data can vary in terms of completeness and temporal coverage, necessitating the use of robust threshold selection techniques to analyze extreme precipitation events effectively.

There are 7795 precipitation records; however, many days have missing values. Figure 8 displays the distribution of records based on the percentage of available data. Only those records with more than 90% data availability are selected, resulting in a total of 1909 locations. In addition, Figure 9 shows the spatial distribution of selected stations, covering different regions of the globe with different densities.

9.2. Methodology

The same threshold selection methods as in the numerical experiment from Section 7 are utilized. The methodology followed is described below:

1. File Filtering: Only files with at least 90% data availability are retained for further processing.

2. Setting Significance Levels: Two significance levels, $\alpha =0.01$ and $\alpha =0.05$, are defined to evaluate the effectiveness of the threshold selection methods.

3. Loop through Filtered Files: For each file, the following procedures are executed:

   • (a). Data Preprocessing: The current file is read, NaN entries are removed, and precipitation data are stored for analysis.

   • (b). Independent Peak Extraction: Using a minimum peak distance of 2 days, independent peaks and their corresponding exceedances are extracted to avoid data dependence based on a predefined threshold.

   • (c). Applying Threshold Selection Methods: For each significance level, the script applies four threshold detection methods:

     • i.  . Studentized Residuals Method: Calculates the optimal threshold by fitting a weighted regression model and computing studentized residuals.

     • ii. . Langousis Method (Local Minima): Detects local minima in the mean square error to identify candidate thresholds, recording up to three values.

     • iii.. Anderson–Darling Test: Applies the Anderson–Darling goodness-of-fit test to identify an optimal threshold.

     • iv.. Cramer–Von Mises Test: Uses the Cramer–Von Mises criterion to evaluate potential thresholds.

(d). Execution Time Logging: The execution time for each method is recorded, and the total time is displayed.

(e). Result Storage: For each file and significance level, the results are stored, including optimal thresholds and execution times.

This automated procedure enables consistent and efficient threshold detection across multiple files, supporting robust analysis of precipitation data.

9.3. Results and Discussion

The application of each threshold selection method to various stations in the GHCN-Daily dataset provides valuable insights into their performance and characteristics. Table 3 presents a summary of key statistical metrics for threshold values identified by the different methods. Each row corresponds to a specific method: Langousis’ local minima (L1, L2, and L3), Studentized Residuals at significance levels of 1% and 5%, and the Anderson–Darling and Cramer–Von Mises (CVM) tests also at 1% and 5%. The summary includes the mean, median, standard deviation, minimum, and maximum values, as well as the first quartile (Q1), third quartile (Q3), and critical lower and upper bounds (Crit_lo and Crit_up) for each method, where $Q_1-1.5\mathrm{IQR}$ and $Q_3+1.5\mathrm{IQR}$ define these bounds. Additionally, the table incorporates a feature representing the percentage of cases where the threshold estimates fall within the range of 2 to 12 mm, as specified by [7]. This metric provides an important benchmark for evaluating the practical suitability of each method.

Results are also shown in Figure 10, where the boxplots associated with each result are displayed.

From the results provided in Table 3 and shown in Figure 10, considerable variation in threshold values across the methods is observed. The Langousis L1 method yields the lowest mean threshold of 10.96 mm/d, indicating that this criterion tends to select lower thresholds compared to others. In contrast, Langousis L3, along with the AD and CVM methods at 5% significance, results in higher mean thresholds, with Langousis L3 reaching a mean of 31.56 mm/d. This broad range suggests that the choice of method significantly influences threshold selection, potentially affecting the reliability and interpretation of extreme-value analysis.

L1 method results are used as reference values because they demonstrate the highest percentage of data falling within the bounds proposed by [7], instilling confidence in its reliability and suitability for threshold selection in precipitation data. Although the numerical experiments revealed instances where the L1 method produced NaN values, this is not an issue in the current analysis, as the percentage of NaN values for L1 is 0%. Even for the L2 and L3 methods, the occurrence of NaN values is minimal, at just 0.3% and 1.2%, respectively. Additionally, the best-performing method in the numerical experiment, the Studentized Residuals (SR) method, performs similarly here with respect to L1. Notably, the case with 1% significance aligns more closely with the results of [7], with only a slight difference observed compared to the 5% significance level. The remaining methods deviate significantly from these results, discarding a large number of peaks and leading to higher threshold estimates. These findings further highlight the advantages of the Langousis L1 and Studentized Residuals methods SR, particularly under strict significance levels, for obtaining accurate and reliable thresholds.

Figure 11 presents a boxplot of the computational time (in seconds) for the Langousis, Studentized Residuals (1% and 5%), Anderson–Darling (1% and 5%), and Cramer–Von Mises (1% and 5%) methods. Each method’s time distribution is displayed, with outliers marked to show variations in computational demands across datasets. The Langousis method consistently shows the shortest computational times, with minimal variance, indicating that it is the most computationally efficient method in this analysis. In contrast, the Anderson–Darling and Cramer–Von Mises methods, particularly at the 5% significance level, exhibit higher mean times and a larger spread, suggesting that these methods are more computationally intensive and may introduce additional variability depending on the data. Overall, the Studentized Residuals methods lie between these two extremes, providing a balanced trade-off between computational efficiency and rigorous threshold selection. In any case, method selection should not be influenced by computational time, as all methods are relatively fast.

The map in Figure 12 presents the spatial distribution of interpolated threshold values derived from the first local minima (L1) of Langousis’ method. Regions with higher thresholds generally correspond to areas of greater precipitation intensity, illustrating the spatial variability in extreme rainfall patterns. Notably, densely populated regions in North America and Europe display moderate thresholds, while areas in the southern hemisphere, such as parts of Africa, are characterized by relatively lower thresholds. It is important to note that the majority of locations with available data report thresholds below 15 mm, with isolated areas showing higher values, ranging between 15 mm and 35 mm (light blue). This highlights distinct geographic differences in rainfall extremes and the distribution of precipitation thresholds.

The maps and boxplots presented in Figure 13 provide an insightful comparison of threshold differences with respect to the first local minimum (L1) of Langousis’ method. In all the maps depicting the differences, a consistent color scale across all figures has been applied to facilitate direct comparison, even though some maps may not exhibit maximum or minimum values that reach the limits of the colormap. In subplot Figure 13a, the differences corresponding to the second local minimum (L2) are generally positive, indicating that the second minimum often results in higher thresholds, with notable regional variations. In both cases, the medians are very close. Subplots Figure 13b,c showcase the differences for the Studentized residuals method at 1% and 5% significance levels, respectively. The map for the 1% significance level shows a mix of very light positive and negative differences, suggesting that this stricter threshold often yields values closer to or below those of L1. For the 5% significance level, the differences tend to be more positive, reflecting the less restrictive nature of this criterion. Across all methods, the associated boxplots highlight the distribution and spread of threshold differences, emphasizing the variability in threshold selection depending on the method and significance level applied.

The maps and boxplots in Figure 14 illustrate the differences in thresholds relative to L1 using the Anderson–Darling (AD) method at 1% and 5% significance levels. The maps reveal spatial heterogeneity in threshold differences, with notable regions showing significant deviations in the southern hemisphere and parts of North America. Please note that the red colors are clearly darker with respect to studentized residual thresholds. The 1% significance level (Figure 14a) generally exhibits smaller differences compared to the 5% level (Figure 14b), reflecting the impact of stricter statistical criteria. The corresponding boxplots provide a statistical summary of the differences, indicating a wider range and higher variability for the 5% level. These results underscore the sensitivity of threshold selection to significance levels, with more lenient criteria (5%) producing larger deviations from L1. These results highlight the considerable increment in the threshold estimation exhibited by the Anderson–Darling (AD) method.

The maps and boxplots in Figure 15 illustrate the differences in thresholds relative to L1 using the Cramer–Von Mises (CVM) method at 1% and 5% significance levels. The maps reveal spatial heterogeneity in threshold differences, with notable regions showing significant deviations in the southern hemisphere and parts of North America, closely mirroring the patterns observed with the Anderson–Darling (AD) method. The 1% significance level (Figure 15a) generally exhibits smaller differences compared to the 5% level (Figure 15b), reflecting the impact of stricter statistical criteria. The corresponding boxplots provide a statistical summary of the differences, indicating a wider range and higher variability for the 5% level. These results underscore the sensitivity of threshold selection to significance levels, with more lenient criteria (5%) producing larger deviations from L1. The close resemblance between the CVM and AD methods highlights the consistent behavior of these goodness-of-fit-based approaches in estimating thresholds.

10. Conclusions

This study focuses on the improvement, development, and application of two threshold selection methods designed to enhance the extraction of extreme precipitation events for Generalized Pareto (GP) distribution fitting. The methods—internally studentized residuals and an adaptation of the Langousis approach—demonstrate robustness in identifying suitable thresholds, particularly for long-duration precipitation records.

The findings align with prior research, emphasizing the importance of selecting an appropriate threshold $u^{*}$ to ensure a reliable GP model fitting for extreme-value analysis. The choice of threshold significantly influences the balance between incorporating sufficient extreme data points and maintaining the assumptions required for accurate parameter estimation. The developed methods address challenges related to data quantization and dependence, which often introduce biases in threshold-based extreme-value analysis.

The proposed methods mitigate subjective bias and eliminate the need for extensive visual inspection, as required in traditional approaches such as the mean residual life plot (MRLP) method. The automation of threshold selection provides a statistically sound alternative, performing reliably even with quantized data. Additionally, the integration of independence checks enhances the reliability of GP modeling by removing potentially dependent peaks, resulting in more accurate estimates for return levels.

In terms of computational performance, the analysis demonstrates that the proposed methods are computationally efficient across datasets with varying levels of data availability. This efficiency is particularly advantageous for large datasets, such as the NOAA NCDC Daily Rainfall Database analyzed in this study.

The results highlight the potential of these methods for applications in hydrology and climate science, where accurate characterization of extreme precipitation events is essential for risk assessment and infrastructure planning. Future research could explore the adaptation of these methods to other types of environmental extremes and further refine the diagnostic criteria to accommodate varying data quality and quantization levels.

In summary, by addressing the limitations of existing threshold selection methods and proposing a solution tailored for high-frequency, quantized data, this study provides a novel toolset for extreme-value analysis, effectively meeting the demands of modern climatological and hydrological datasets.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The MATLAB code used to download and reproduce all the results presented in this manuscript is available in the following GitHub repository: https://github.com/rminguez/Matlab_Threshold_Selection_Project (accessed on 9 December 2024). Detailed information about the project content is provided in Appendix A.

Acknowledgments

I would like to express my gratitude to the three anonymous reviewers for their constructive comments and insightful suggestions, which have significantly improved the quality of this manuscript.

Conflicts of Interest

The author declares no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

GHCN	Global Historical Climatology Network
GoF	Goodness of Fit
AM	Annual Maxima
AMS	Annual Maxima Series
PDS	Partial Duration Series
POT	Peaks Over Threshold
EVT	Extreme-Value Theory
MRLP	Mean Residual Life Plot
MSE	Mean Square Error
AD	Anderson–Darling
CVM	Cramer–Von Mises
GPD	Generalized Pareto Distribution
GP	Generalized Pareto
CDF	Cumulative Distribution Function
MLE	Maximum Likelihood Estimation
RWLS	Robust Weighted Least Squares
SR	Studentized Residuals
L1	First Local Minimum (Langousis Method)
L2	Second Local Minimum (Langousis Method)
L3	Third Local Minimum (Langousis Method)
PRCP	Precipitation (Dataset Identifier)
GHCN-Daily	Daily Version of the Global Historical Climatology Network Dataset
NOAA	National Oceanic and Atmospheric Administration
NCDC	National Climatic Data Center
MCMC	Markov Chain Monte Carlo
PDF	Probability Density Function
SE	Standard Error

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Table 1

Performance of threshold selection methods across rounding levels. The table reports the percentage of NaN values (NaN %) and the percentage of deviations within 1 (<1%) for each method under different rounding levels.

Metric	Rounding	L1	L2	L3	SR 1%	SR 5%	AD 1%	AD 5%	CVM 1%	CVM 5%
NaN %	1	24.1	59.2	80.3	0.0	0.0	0.0	0.0	0.0	0.0
	2	2.5	23.7	55.9	0.0	0.0	0.0	0.0	0.0	0.0
	3	0.2	2.9	27.0	0.0	0.0	0.0	0.0	0.0	0.0
<1%	1	41.0	35.9	4.8	58.9	72.5	46.8	53.3	41.7	51.6
	2	62.2	72.2	23.3	58.4	63.8	37.1	44.3	31.9	42.9
	3	54.3	78.6	52.2	53.3	53.5	34.0	43.2	31.3	39.8

Table 2

Performance of threshold selection methods across shape parameter ($\xi$) levels. The table reports the percentage of NaN values (NaN %) and the percentage of deviations within 1 (<1%) for each method under different $\xi$ values.

Metric	$\mathit{\xi }$	L1	L2	L3	SR 1%	SR 5%	AD 1%	AD 5%	CVM 1%	CVM 5%
NaN %	0.15	10.7	28.5	48.5	0.0	0.0	0.0	0.0	0.0	0.0
	0.10	9.6	30.0	54.1	0.0	0.0	0.0	0.0	0.0	0.0
	0.05	10.0	30.0	55.6	0.0	0.0	0.0	0.0	0.0	0.0
	0.00	10.7	29.3	54.1	0.0	0.0	0.0	0.0	0.0	0.0
	−0.05	8.1	26.3	54.8	0.0	0.0	0.0	0.0	0.0	0.0
	−0.10	5.6	28.1	55.9	0.0	0.0	0.0	0.0	0.0	0.0
	−0.15	7.8	27.8	57.8	0.0	0.0	0.0	0.0	0.0	0.0
<1%	0.15	54.8	59.3	31.5	56.7	62.2	40.7	28.5	24.4	37.4
	0.10	54.1	61.9	26.3	55.2	59.6	36.7	27.4	24.4	33.3
	0.05	61.5	60.4	23.7	57.4	64.4	41.1	35.2	32.6	39.6
	0.00	49.6	60.7	25.9	59.6	65.9	51.9	42.6	37.4	50.0
	−0.05	53.0	60.0	27.8	58.9	64.8	41.1	36.3	36.3	45.6
	−0.10	50.0	65.6	26.3	57.8	64.8	48.9	42.2	42.2	53.7
	−0.15	44.4	67.8	25.9	52.6	61.1	56.7	51.5	47.4	53.7

Table 3

Summary Statistics for Threshold Selection Methods for stations in the GHCN-Daily dataset, including a percentage of values meeting specific criteria.

Method	Mean	Median	StdDev	Min	Max	Q1	Q3	Crit_lo	Crit_up	Percentage (%)
L1	10.96	7.69	10.67	0.28	76.39	3.16	14.92	0.00	32.56	53.54
L2	21.97	18.71	13.79	2.55	107.21	11.87	28.04	0.00	52.29	25.64
L3	31.56	27.71	16.23	4.71	139.04	20.27	39.67	0.00	68.77	4.98
SR 1%	8.96	7.10	9.81	0.10	109.50	2.00	12.40	0.00	28.00	49.55
SR 5%	11.88	8.10	12.16	0.10	109.50	4.30	16.30	0.00	34.30	46.05
AD 1%	27.38	25.10	13.60	6.10	191.00	18.08	33.65	0.00	57.01	4.40
AD 5%	30.43	27.70	16.17	6.20	191.00	19.60	37.30	0.00	63.85	2.72
CVM 1%	26.39	23.60	13.31	6.10	191.00	17.30	32.43	0.00	55.11	5.66
CVM 5%	29.33	26.70	15.56	6.10	191.00	19.18	35.80	0.00	60.74	3.46

Appendix A. Code Availability

The MATLAB code used to reproduce all the results presented in this manuscript is available in the following GitHub repository: https://github.com/rminguez/Matlab_Threshold_Selection_Project (accessed on 9 December 2024). The repository is open-source and distributed under the MIT License, allowing free use and modification. The code has been tested primarily in MATLAB 2016, ensuring compatibility with this version. However, it is expected to work with later versions, though minor adjustments may be necessary.

Repository Structure

The repository contains the following main components:

• Scripts: The main scripts for setting up the data processing workflow and running threshold selection analyses. The key scripts include:

  • –

    Script_Download_Rainfall_Database.m: Downloads and processes the GHCN-Daily precipitation data, including filtering stations with more than 110 years of data.

  • –

    Script_GHCN_DAILY_Database.m: Applies threshold selection methods to the GHCN-Daily dataset and organizes data availability results.

  • –

    Script_GHCN_DAILY_AnalysesResults.m: Analyzes threshold selection outcomes and generates statistical summaries and boxplots.

  • –

    Script_New_PRCP_ASN00021043.m: Illustrates the application of the threshold selection methods to a specific rainfall station (PRCP_ASN00021043).

  • –

    Script_NewSimulationExperiment.m: Conducts new simulation experiments to assess the performance of threshold selection methods under varying conditions, including synthetic datasets with different parameters.

  • –

    Script_NewSimExpAnalysisResults.m: Analyzes the results of the new simulation experiments, comparing the performance of threshold selection methods and generating detailed summaries, plots, and statistical evaluations.

• Functions: A collection of modular functions performing core analyses, including threshold selection methods:

  • –

    threshold_peak_extraction.m: Identifies independent peaks above specified thresholds in a time series, filtering based on a minimum peak distance.

  • –

    threshold_studentized_residuals.m: Implements threshold selection based on internally studentized residuals and weighted least squares regression.

  • –

    threshold_MSE.m: Uses Langousis’ method to identify optimal thresholds by locating local minima in the mean square error spline fit.

  • –

    threshold_AD.m: Selects thresholds using the Anderson–Darling test for the Generalized Pareto Distribution fit.

  • –

    threshold_CVM.m: Similar to threshold_AD.m, but uses the Cramer–Von Mises statistic to identify optimal thresholds.

  • –

    RWLSfit.m: A helper function that performs weighted least squares regression, used by several threshold selection methods.

• Data and Figures: Organized into data and figures folders to store raw data, processed data, and output visualizations, respectively.

  • –

    The data folder includes files for the GHCN-Daily inventory, data availability lists, and precipitation data for individual stations.

  • –

    The figures folder contains output figures, such as the mean residual life plots, boxplots comparing threshold selection methods, and any supplementary visualizations generated by the scripts.

This repository provides a complete and reproducible environment for researchers interested in threshold selection techniques for extreme-value analysis. For more details on specific function usage, refer to the comments within each script and function file.

References

1. Coles, S. An Introduction to Statistical Modeling of Extreme Values; Springer: Berlin/Heidelberg, Germany, 2001.

2. Pickands, J. Statistical inference using extreme order statistics. Ann. Stat.; 1975; 3, pp. 119-131.

3. Davison, A.; Smith, R. Models for exceedances over high thresholds. J. R. Stat. Soc. Ser. B (Methodol.); 1990; 52, pp. 393-442. [DOI: https://dx.doi.org/10.1111/j.2517-6161.1990.tb01796.x]

4. Grimshaw, S. Computing maximum likelihood estimates for the generalized Pareto distribution. Technometrics; 1993; 35, pp. 185-191. [DOI: https://dx.doi.org/10.1080/00401706.1993.10485040]

5. Beirlant, J.; de Wet, T.; Goegebeur, Y. A goodness-of-fit statistic for Pareto-type behaviour. J. Comput. Appl. Math.; 2006; 186, pp. 99-116. [DOI: https://dx.doi.org/10.1016/j.cam.2005.01.036]

6. Hosking, J.; Wallis, J. Parameter and quantile estimation for the generalized Pareto distribution. Technometrics; 1987; 29, pp. 339-349. [DOI: https://dx.doi.org/10.1080/00401706.1987.10488243]

7. Langousis, A.; Mamalakis, A.; Puliga, M.; Deidda, R. Threshold detection for the generalized Pareto distribution: Review of representative methods and application to the NOAA NCDC daily rainfall database. Water Resour. Res.; 2016; 52, pp. 2659-2681. [DOI: https://dx.doi.org/10.1002/2015WR018502]

8. Serinaldi, F.; Kilsby, C. Rainfall extremes: Toward reconciliation after the battle of distributions. Water Resour. Res.; 2014; 50, pp. 336-352. [DOI: https://dx.doi.org/10.1002/2013WR014211] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24829512]

9. Papalexiou, S.; Koutsoyiannis, D.; Makropoulos, C. How extreme is extreme? An assessment of daily rainfall distribution tails. Hydrol. Earth Syst. Sci.; 2013; 17, pp. 851-862. [DOI: https://dx.doi.org/10.5194/hess-17-851-2013]

10. Anderson, T.; Darling, D. Asymptotic theory of certain goodness-of-fit criteria based on stochastic processes. Ann. Math. Stat.; 1952; 23, pp. 193-212. [DOI: https://dx.doi.org/10.1214/aoms/1177729437]

11. Menne, M.; Durre, I.; Vose, R.; Gleason, B.; Houston, T. An overview of the Global Historical Climatology Network-Daily Database. J. Atmos. Ocean. Technol.; 2012; 29, pp. 897-910. [DOI: https://dx.doi.org/10.1175/JTECH-D-11-00103.1]

12. Menne, M.; Durre, I.; Korzeniewski, B.; McNeill, S.; Thomas, K.; Yin, X.; Anthony, S.; Ray, R.; Vose, R.; Gleason, B. et al. Global Historical Climatology Network-Daily (GHCN-Daily), Version 3.31; NOAA National Centers for Environmental Information: Lincoln, NE, USA, 2012; [DOI: https://dx.doi.org/10.7289/V5D21VHZ]

13. Mínguez, R.; Méndez, F.J.; Izaguirre, C.; Menéndez, M.; Losada, I.J. Pseudo-Optimal Parameter Selection of Non-Stationary Generalized Extreme Value Models for Environmental Variables. Environ. Model. Softw.; 2010; 25, pp. 1592-1607. [DOI: https://dx.doi.org/10.1016/j.envsoft.2010.05.008]

14. Mínguez, R.; Menéndez, M.; Méndez, F.J.; Losada, I.J. Sensitivity Analysis of Time-Dependent Generalized Extreme Value Models for Ocean Climate Variables. Adv. Water Resour.; 2010; 33, pp. 833-845. [DOI: https://dx.doi.org/10.1016/j.advwatres.2010.05.003]

15. Salas, J.D.; Anderson, M.L.; Papalexiou, S.M.; Francés, F. PMP and Climate Variability and Change: A Review. J. Hydrol. Eng.; 2020; 25, 03120002. [DOI: https://dx.doi.org/10.1061/(ASCE)HE.1943-5584.0002003]

16. Castillo, E.; Castillo, C.; Mínguez, R. Use of Extreme Value Theory in Engineering Design. Safety, Reliability and Risk Analysis: Theory, Methods and Applications, Proceedings of the European Safety and Reliability Conference 2008 (ESREL 2008), Valencia, Spain, 22–25 September 2008; Martorel, S.; Guedes Solares, C.; Barnett, J. Taylor & Francis Group: Abingdon, UK, 2008; Volume 3, pp. 2473-2488.

17. Draper, N.R.; Smith, H. Applied Regression Analysis (Wiley Series in Probability and Statistics); 2nd ed. John Wiley & Sons Inc.: Hoboken, NJ, USA, 1981.

18. Belsley, D.A.; Kuh, E.; Welsch, R.E. Regression Diagnostics: Identifying Influential Data and Sources of Multicollinearity; John Wiley & Sons: New York, NY, USA, 1980.

19. Chatterjee, S.; Hadi, A.S. Influential Observations, High Leverage Points, and Outliers in Linear Regression. Stat. Sci.; 1986; 1, pp. 379-393.

20. Mínguez, R.; Reguero, B.G.; Luceño, A.; Méndez, F.J. Regression Models for Outlier Identification (Hurricanes and Typhoons) in Wave Hindcast Databases. J. Atmos. Ocean. Technol.; 2012; 29, pp. 267-285. [DOI: https://dx.doi.org/10.1175/JTECH-D-11-00059.1]

21. Chatterjee, S.; Hadi, A.S.; Price, B. Regression Analysis by Example; John Wiley & Sons: New York, NY, USA, 1990.

22. Abur, A.; Expósito, A.G. Electric Power System State Estimation. Theory and Implementations; Marcel Dekker: New York, NY, USA, 2004.