Using Simple Fixed-Point Iterations to Estimate Generalized Pareto Distribution Parameters

Abstract

Estimating generalised Pareto distribution (GPD) parameters is a fundamental step in modelling the extremevalue distribution of random variables. It is generally done with the maximum likelihood method, but there are generally difficulties in estimating GPD parameters using this method as there is no closed-form solution for the first derivative of the GPD log-likelihood function. This makes the solution difficult to determine analytically. However, numerical methods can be used as an alternative. Therefore, this study estimates the solution numerically using a simple fixed-point iteration method that is intuitive for both practitioners and professionals. We obtained three fixed-point iterations when estimating GPD parameters that met the unbiased estimator and convergence criteria. The iterations allow practitioners and professionals to directly and efficiently estimate GPD parameters when modelling extreme-value distributions of random variables.

Headnote

Abstract-Estimating generalised Pareto distribution (GPD) parameters is a fundamental step in modelling the extremevalue distribution of random variables. It is generally done with the maximum likelihood method, but there are generally difficulties in estimating GPD parameters using this method as there is no closed-form solution for the first derivative of the GPD log-likelihood function. This makes the solution difficult to determine analytically. However, numerical methods can be used as an alternative. Therefore, this study estimates the solution numerically using a simple fixed-point iteration method that is intuitive for both practitioners and professionals. We obtained three fixed-point iterations when estimating GPD parameters that met the unbiased estimator and convergence criteria. The iterations allow practitioners and professionals to directly and efficiently estimate GPD parameters when modelling extreme-value distributions of random variables.

Index Terms-generalized Pareto distribution; parameter estimation; maximum likelihood method; fixed-point iteration; unbiasedness criteria; convergence criteria

(ProQuest: ... denotes formulae omitted.)

I. Introduction

XTREME value theory (EVT) is a unique approach to determining the probability of the extreme value of a random variable on the left or right tail of the probability function [1]. This approach is widely used in various fields, such as climatology, finance, and insurance. In general, studies that apply EVT are more likely to analyse the probability of extreme values of random variables in the right tail of the probability function since this is more interesting [2]. This includes studies on maximum rainfall risk [3], maximum temperature risk [4], [5], [6], maximum sea wave height [7], [8], maximum disaster loss risk [9], [10], and maximum earthquake magnitude risk [11], [12], [13].

There are two methods for measuring the extreme-value probability of random variables using the EVT approach: 1) the block maxima (BM) method; and 2) the peaks-over threshold (POT) method. The latter is the most modern method of the two, and it facilitates the design of the cumulative-distribution function (CDF) model for extreme values in which the 'father' CDF of a random variable is unknown [2]. The advantage of this method Fes in its abifity to define extreme values based on values that exceed a particular threshold [14]. This definition often makes the size of extreme values in the POT method more significant than those in the BM method, indicating that the POT method is more representative of the population and that the probabifity results will be closer to the actual probabifity.

One of the fundamental steps in modelfing the right-tailed probabifity function of a random variable using the POT method entails estimating generalized Pareto distribution (GPD) parameters. GPD is a two-parameter probabifity distribution that is crucial in modelfing. Many methods can be used to estimate GPD parameters, but the most popular is the maximum likefihood (ML) method [15]. This method is intuitive as the estimated parameter sought is a twodimensional vector that maximizes the GPD likefihood function. If the value of the GPD likefihood function reaches the maximum value, the sample that represents the population will very likely be observed. However, there are obstacles to determining the vector solution due to the first derivative of the GPD log-likefihood function (i.e., the natural logarithmic form of the GPD likefihood function) not having a closed form. Hence, it is difficult to determine the solution for the local-maximum value using ordinary analytical methods.

The difficulty in determining this solution can be overcome by numerical methods. Identifying a solution that maximizes the GPD log-likefihood function is similar to finding the root of its first derivative when it has a value equal to zero. Therefore, these numerical methods are categorized as rootfinding methods. Several numerical root-finding methods, such as Newton, secant, and fixed-point iteration, can be used, of which fixed-point iteration is the simplest basic method [16]. This method determines the root as the fixed point of a function [17], denoting that only sfight manipulation of the function's form is required and that further differentiation is not needed, as is the case in the Newton or secant methods [18].

Thus, this study's research objective is to develop novel fixed-point iteration equations to estimate GPD parameters via the ML method. The fiterature review in Section II demonstrates that fixed-point iteration equations have not been studied in terms of GPD parameter estimation. Moreover, we also review the equations for accuracy and goodness of fit using case studies on the disaster loss data from the United States between 1977 and 2021. This study will allow professionals and practitioners to estimate GPD parameters in a way that the design of the extreme distribution ofa random variable can be earned out directly and easily using the POT method.

II. Literature Review

Several studies were earned out on estimating GPD parameters between 1987 and 2000. For example, Hosking and Wallis [19] estimated GPD parameters using the method of moments (MOM) and the probability-weighted moment (PMW) approaches. They also compared these methods using Monte Carlo simulations. Moharram et al. [20] used the least squares (LS) method to estimate GPD parameters, after which they compared the estimation results with MOM and PWM methods based on root mean square error. The comparison showed that the estimated results deduced from the LS method were generally smaller than those deduced using the MOM and PWM methods. Furthermore, Singh and Guo [21] implemented the principle of maximum entropy (POME) method to estimate GPD parameters and found that the estimated results were better than those for the MOM and PMW methods in several data ranges. In addition, Lin and Wang [22] estimated GPD parameters using the ML method on censored data.

Several studies on estimating GPD parameters were also earned out from 2000 to 2022, with de Zea Bermudez and Turkman [23] first estimating GPD parameters using the Bayesian method and showing that this method was suitable for small data sizes in simulations. Juarez and Schucany [24] applied the minimum density power divergence estimator method to estimate GPD parameters, which demonstrated that the coefficients for the robustness and the efficiency of the estimated GPD parameters can be controlled. Zhang [25] estimated GPD parameters using likelihood-moment estimation (LME) and showed the high asymptotic efficiency that resulted from the parameter estimates. Zhang and Stephens [26] used the ML and quasi-Bayesian methods to estimate the GPD parameters. Following this, de Zea Bermudez and Kotz [27] estimated GPD parameters by combining robust and Bayesian techniques. Moreover, Mackay et al. [7] estimated GPD parameters using LME, which has a low sensitivity for selecting a threshold value and is suitable for data sizes lower than 500. Song and Song [28] estimated GPD parameters using the nonlinear least square method. Wang and Chen [29] introduced a hybrid estimation method that minimized the statistical value of the Anderson-Darling (AD) test and maximized the GPD loglikelihood function. In addition, El-Sagheer et al. [30] assessed GPD parameters on progressive and censored failure data using the Bayesian method, the performance of which was determined by Monte Carlo simulations. Finally, Form and Ratnasingam [31] estimated GPD parameters using the elemental percentile method.

None of these studies estimated GPD parameters using the ML method in the form of fixed-point iterations. Thus, this study makes a novel contribution to the literature.

III. Estimating GPD Parameters Using the ML Method

The GPD was first introduced by Pickands [32] and Balkema and de Haan [33]. Smith [34], Davison [35], Smith [36], and van Montfort and Witter [3] then contributed to the development of this method through their research. For example, X is a random variable with a fat-tailed CDF, ц is the threshold value between the extreme and non-extreme values of X, and Y = X - jU with X > и denotes an extreme-excess random variable. Thus, if the value of ц is large, then the CDF of Y can be approximated by the CDF of the GPD, which is expressed as follows:

...

In this equation, ст > 0 and shape parameters, respectively. Moreover, the probability density function (PDF) of Y, approximated by the PDF of the GPD, is stated as follows:

...

when < 0 [37]. Based on its shape parameters, Salvadori et al. [2] stated that there are three types of GPD (see Table I). The shape parameter values in this study are assumed to be non-zero to limit the scope of the study.

The most popular method for estimating GPD parameters is the ML method. Herein, Xj with y=l,2,...,n are independent and identically distributed random variables for which the PDF tails are fat. Moreover Y = X¡-u, where X¡ > и Xi > u, i = i2,...,m, represents the extreme-excess random variables that are independent and identically distributed in the GPD. Thus, the likelihood function of Y, is expressed as follows:

...

The first derivative of (4) with respect to each element of 0 is expressed as follows:

...

Thus, the solution for (6) is difficult to estimate analytically since it does not have a closed form, but a numerical method is a fitting alternative for solving this problem. Therefore, this study uses the fixed-point iteration method to do so.

IV. Main Results

A. Designing Possible Fixed-Point Iterations in GPD Parameter Estimation using the ML Method

Fixed-point iteration is a method used to numerically determine the root of a function. This method can also be used to determine the maximum solution of a function, which is the root of its first derivative. Assuming that the function for which the maximum solution is sought is (4), using the fixed-point iteration method, (6) is first modified to form the equation 0 = g(0|y). This transformation results in the following fixed-point iteration equation:

...

If et is less than e , where e is the tolerance error, the iteration is stopped. To start the iteration for the GPD parameter estimation, (7) must first be designed. It is important that the equation satisfies the unbiased estimator and convergence requirements. Firstly, the equations for the possible iterations of at and z were designed separately. These possible iterations are presented in Table II, which demonstrates that there were eight possible schemes of combination for at and z. Secondly, each scheme's unbiased estimators and convergence were examined.

B. Determining the Unbiased Estimators c/ the GPD Parameter Estimation

The unbiasedness of a as the estimator of 0 was investigated by determining each element's unbiased estimator, a is an unbiased estimator of 0 if

...

where Bias (ст, cr) is the unbiased estimator of both à and er as well as where Bias5,5 is the unbiased estimator of ý. If a p t with /> = 1,2 is the iteration dt of type p in Table II, and if 5 with q = 1,2,3,4 represents the iteration 5t of type q in Table II, the unbiased estimator of dXt is determined first. Thus, the unbiased estimator of dXt must be zero.

...

If Biased,d) = 0 , d is the unbiased estimator of cr, the unbiased estimator of £1( must be reviewed. Herein, this unbiased estimator of must be zero.

...

As Bias5,5 = 0, i is the unbiased estimator of ý, the processes of determining the unbiased estimators of d2t, 52t,53t, and ¿41 are generally the same as those for cq t and 1. The results of these checks demonstrate that only 53t is biased since it results in Bias 5,5 = -25 · 0 * Therefore, 5Ъ t was omitted from the schemes, resulting in the iteration in which /7 = 1,2 and q· = 1,2,4, are the unbiased estimators of 6 .

C. Determining the Convergence c/ the GPD Parameters ' Iterations

The convergence of the iteration equation ØgJ was determined using Theorem 1.

Theorem 1. Let g(0>) and g'(0>) be continuous for a two-dimensional closed set T> e R2, and let

...

Moreover, if

...

then there is a unique solution a for >9 = g(i9|y), and the iteration 9t will converge to a for any initial estimators 90eP.

Proof. (16) shows that there is at least one solution a for ø = g(ø|y) in D. Thus, the solution for a is unique. Contrastingly, if there are two solutions for 0 = g(M' namely a, and a2, using the mean value theorem, it can be seen that

...

where c e cq xa2. If the absolute value of (18) is taken and (17) is then used, we obtain the following:

...

Note that

...

Since /<1, aj is the same as a2, the solution for a is unique. Following this, the iteration 9t will converge to a for any initial estimators 60eP. Based on (16), it is inductively clear that 90 e P,9Z e P . This is similar to (18), which demonstrates that

...

where d e a x 9Z. Thus, equation (21 ) can be inductively generalized as follows:

...

Since / < 1, /z ļļa - 901] -> 0 as t -> oo , this demonstrates that 9Z -> a as í -> « .

In practice, the convergence of iteration 9Z using Theorem 1 is rarely determined. This is because it is not easy to find the closed set P e K2 that satisfies (16). Therefore, Theorem 1 can be practically used by applying Corollary 1.

Corollary 1. If g(0|y) and g'(0[y) are continuous for a two-dimensional closed set P e K2, in which the solution a of 9 = g(0|y) is included, and if

...

there is a closed neighbourhood K, for a , thereby proving the results of Theorem 1.

Proof. Using assumption (23), we can create a neighbourhood for /C with a radius e > 0 around a, ( 21 1 * )C= keR'a-k<fļ, which satisfies the following conditions:

...

It must be shown that g(k|y)e/CVke/C. If any element k e /C is selected, it must be noted that

...

where seaxk. Equation (25) indicates that g(k|y)e/C. Finally, the iteration 9Z converging to a for any initial estimator 90 e /С can be justified in the same ways as the methods used in Theorem 1, except by replacing P with K.. □

In this study, the convergence for both cr and £ „ was P" q J determined using Corollary 1. To facilitate this determination, we used the data on the disaster losses in the United States from 1977 to 2021 in billions of dollars, with a data size of 280 being obtained. This information was accessed at https://www.emdat.be (accessed 8 November 2022). The data were confirmed to have extreme values since the probability function tail was fat and sloped more to the right. The fat tail of the probability function can be seen in the data's Kurtosis (80.1725), which was >3, and the slope can be seen in the data skewness (8.0095), which was positive.

To obtain the vector у , the threshold value ц was first determined. Herein, the threshold value was determined using the Kurtosis method. The threshold value selection algorithm is presented in Fig. 1.

The data were first partitioned into 280 sub-data sets. The first sub-data set was the data itself, and the second sub-data set was the former sub-data set without the most prominent datum. The threshold value was the most prominent datum of the largest sub-data set that had a Kurtosis of <3. The threshold value was determined using the Kurtosis method in Scilab v. 6.1.1. The threshold value obtained was 4.4881, and the most prominent datum was the 196th sub-data set, which was the largest sub-data set, with a Kurtosis of <3, namely 2.9399. The vector of у is presented as follows:

...

The next stage comprised running each iteration pair. The algorithm for determining the convergence with the data is presented in Fig. 2. The initial vector estimator and the tolerance error were determined in advance. According to Hosking and Walks [19], the initial vector can be determined using the following equation:

...

Alternatively, the initial vector can also be determined using the following equation [15]:

...

where s1 represents the variance of y. This study determined the initial vector estimator using equation (27). Following this, the selected tolerance error was e = I x 10 6 .

...

D.Fixed-Point Iterations in GPD Parameter Estimation using the ML Method

These results indicate that three iterations meet the unbiased estimator and convergence conditions, namely with q· = 1,2,4 . Expanding on this, the iteration

for each q· is expressed as follows:

...

Therefore, professionals and practitioners can use Equations (30), (31), and (32) to estimate GPD parameters using the ML method.

V. Discussion

A. The Unbiased Estimator Is Not Necessarily Convergent

The determination of the unbiased estimators of the GPD parameters in Section IV (B) shows that iteration with q· =1,2,4 is an unbiased estimate of 0, and it is justified by each unbiased estimator that is equal to zero. However, after determining the convergence of iteration in Section IV (C), only iteration to not meet the convergence criterion. Thus, an unbiased estimator of the GPD parameter is not necessarily convergent.

B. The Effect cf the |k(a|y)|| Value on Convergence Speed

As can be seen in Table III, the value of |k(a|y)|| has a positive relationship with the number of iterations required to converge to the solution. The closer to zero the value of M«iy)ii is, the fewer iterations it takes to converge, and vice versa.

C. The Estimated GPD Parameters ' Accuracy and Goodness ć/ Fit

The estimated GPD parameters' accuracy can be determined through various measures. Herein, the mean absolute percentage error (MAPE) was used. The accuracy of the GPD parameter estimation was considered 'good' if the MAPE value was between 10%-20%, and it was considered 'very good' if the MAPE values were <10% [41], [42]. The MAPE of the estimated GPD parameters was calculated using the following equation:

...

where Fe represented the empirical distribution function and Fy 1 (*) represented the inverse function of (1). The MAPE of the estimated GPD parameters from each iteration I e is presented in Table V, which shows that this MAPE value was between 10%-20%. Therefore, the GPD parameter estimation for each iteration classified as 'good'.

Furthermore, the GPD parameters' goodness of fit can be determined visually or formally. Visually, determining goodness of fit can be conducted with probability-probability plots (P-P plots), where set = 1,2,...,mļ is placed into a Cartesian diagram. The estimated GPD parameters are considered visually fit for describing the data if the data for set A are scattered around a line with a gradient of 1 [43].

The P-P plots for each estimated GPD parameter of iteration are presented in Fig. 3, which demonstrates that the scatterplot of set A for each iteration around a red Une with a gradient of one. Thus, the GPD parameters estimated by iteration describe the data distribution. Following this, the goodness of ht can be formally reviewed using the KolmogorovSmirnov (KS) and AD tests. The statistical value of the KS test is determined using the following equation [44], [45], [46], [47]:

...

If the statistical value of the KS test is greater than the critical value (ra) , the estimated GPD parameter is not fit for describing the data distribution, and vice versa. The statistical value of the AD test can be determined using the following equation [48]:

...

If the statistical value of the AD test is greater than the critical value ( Aa ), the estimated GPD parameter is not fit for describing the data distribution, and vice versa. With a significance level of 0.05, the statistical values of the KS and AD tests and their respective critical values for the estimated GPD parameters for each iteration Table VI, which shows that the statistical values of the KS and AD tests for each iteration e were smaller than their respective critical values. Therefore, the GPD parameters of each iteration

D. The Necessary Conditions for Reaching the Maximum Solution for the Estimated GPD Parameters

The estimated GPD parameters for each iteration must meet the necessary conditions to reach the maximum solution of the log-likelihood function (4). The maximum local solution of the GPD log-likelihood function (4) is reached at the critical point a if the Hessian matrix of (4), denoted by H(6|y), at a is a negative definite, namely [49]:

...

The results from reviewing the necessary conditions for reaching the maximum solution for the estimated GPD parameters for each iteration „ are presented in Table VII, which delineates that the aTH(a | y)a value of each was negative for each a eK2-{0}. Therefore, every a in each iteration maximum solution of the log-likelihood function (4).

VI. Conclusion

This study presented the design of fixed-point iteration equations for estimating GPD parameters using the ML method. The fixed-point iteration method was used because it is the simplest numerical method for determining a function's solution. This method was also used to determine the solution for the first derivative of the GPD log-likelihood function, which does not have a closed form.

In the initial stage of designing the iterations, eight fixedpoint iteration equations were obtained (see Table II). After this, the eight fixed-point iteration equations were reviewed to determine their unbiased estimators and convergence, with three of them satisfying these conditions, as presented in equations (30), (31), and (32). Thus, numerical proof was obtained for the assertion that an unbiased iteration equation of the GPD parameter is not necessarily convergent to the solution. Moreover, the closer to zero the value of the first derivative of the iteration at the solution, the fewer iterations it takes to converge. Moreover, it was found that the estimated GPD parameters obtained through the simulation have good accuracy and goodness of fit. The accuracy was demonstrated by the MAPE value of 17%, while the goodness of fit was determined using KS and AD tests, which resulted in values less than their respective critical values. Finally, the estimated GPD parameters met the necessary conditions for the maximum solution of the GPD log-likelihood function.

Therefore, the fixed-point iteration equations designed in this study can facilitate the estimation of GPD parameters when modelling the extreme distribution of random variables. By using these fixed-point iteration equations, professionals and practitioners no longer need to design fixed-point iteration equations when estimating GPD parameters and when reviewing unbiased estimators and convergence.

Sidebar

Manuscript received 13 February 2023; revised 26 November 2023.

This work was supported in part by the Rector and the Directorate of Research and Community Service (DRPM) Universitas Padjadjaran.

References

References

[1] S. Thendo, N. Nethengwe, C. Sigauke and H. Chikoore, "Modelling of extremely high rainfall in Limpopo province of South Africa," Climate, vol. 10, no. 3,p. 33, 2022.

[2] G. Salvadori, C. D. Michele, N. T. Kottegoda and R. Rosso, Extremes in Nature: An Approach Using Copulas. Dordrecht: Springer, 2007.

[3] M. A. J. van Montfort and J. V. Witter, "Testing exponentiality against generalised Pareto distribution," J. Hydrol., vol. 78, no. 3-4, pp. 305- 315,1985.

[4] J. Ferrez, A. C. Davison and M. Rebetez, "Extreme temperature analysis under forest cover compared to an open field," Agric. For. Meteorok, vol. 151, no. 7, pp. 992-1001, 2011.

[5] X. L. Wang, B. Trewin, Y. Feng and D. Jones, "Historical changes in Australian temperature extremes as inferred from extreme value distribution analysis," Geophys. Res. Lett., vol. 40, no. 3, pp. 573-578, 2013.

[6] F. Ashkar and C. N. Tatsambon, "Revisiting some estimation methods for the generalized Pareto distribution," Jour. Hydrol., vol. 346, no. 34,pp. 136-143,2007.

[7] E. B. L. Mackay, P. G. Challenor and A. S. Bahaj, "A comparison of estimators for the generalized Pareto distribution," Ocean Eng., vol. 38, no. ll-12,pp. 1338-1346,2011.

[8] P. Stansell, "Distributions of extreme wave, crest and trough heights measured in the North Sea," Ocean Eng., vol. 32, no. 8-9, pp. 1015- 1036,2005.

[9] L. Wei, L. Liu and J. Hou, "Pricing hybrid-triggered catastrophe bonds based on copula-EVT model," Quant. Fin. and Econ., vol. 6, no. 2, pp. 223-243,2022.

[10] W. Chao and H. Zou, "Multiple-event catastrophe bond pricing based on CIR-Copula-POT model," Discrete Dyn. Nat. Soc., vol. 2018, 2018, Art. no. 5068480.

[11] A. A. Zimbidis, N. E. Frangos and A. A. Pantelous, "Modeling earthquake risk via extreme value theory and pricing the respective catastrophe bonds," ASI IN Bulk, vol. 37,no. l,pp. 163-183,2007.

[12] J. Shao, A. Pantelous and A. D. Papaioannou, "Catastrophe risk bonds with applications to earthquakes," Eur. Actuar. J., vol. 5, no. 1, pp. 113-138,2015.

[13] L. Huyse, R. Chen and J. A. Stamatakos, "Application of generalized Pareto distribution to constrain uncertainty in peak ground accelerations," Bulk Seismok Soc. Am., vol. 100, no. 1, pp. 87-101, 2010.

[14] R. A. Ibrahim, S. Sukono and R. Riaman, "Estimation of the extreme distribution model of economic losses due to outbreaks using the POT method with Newton Raphson iteration," Int. J. Quant. Res. Mod., vol. 2, no. 1, pp. 37 45, 2021.

[15] P. de Zea Bermudez and S. Kotz, "Parameter estimation of the generalized Pareto distribution-Part I," J. Stat. Plan. Irference, vol. 140, no. 6, pp. 1353-1373, 2010.

[16] K. E. Atkinson, An Introduction to Numerical Analysis, 2nd ed. New Jersey: John Wiley & Sons, 1989.

[17] K. Atkinson and W. Han, Elementary Numerical Analysis, 3rd ed. New Jersey: John Wiley & Sons, 2004.

[18] S. Purwani, A. F. Ridwan, R. A. Hidayana and S. Sukono, "Secant method with Aitken extrapolation outperform Newton-Raphson method in estimating stock implied volatility," 1AENG Int. J. Comput. Sci,, vol. 50, no. 2, pp. 368-374, 2023.

[19] J. R. M. Hosking and J. R. Wallis, "Parameter and quantile estimators for generalized Pareto distribution," Technometrics, vol. 29, no. 3, pp. 339-349, 1987.

[20] S. H. Moharram, A. K. Gosain and P. N. Kapoor, "A comparative study for the estimators of the generalized Pareto distribution," J. Hydrol., vol. 150, no. l,pp. 169-185, 1993.

[21] V. P. Singh and H. Guo, "Parameter estimation for 2-parameter generalized Pareto distribution by POME," Stoch. Hydrol. Hydrauk, vol. 9, pp. 81-93, 1995.

[22] W. T. Lin and W. Y. Wang, "Estimation for the generalized Pareto distribution with censored data," Commun. Stat. Sim. Comput., vol. 29, no. 4,pp. 1183-1213,2000.

[23] P. de Zea Bermudez and M. A. A. Turkman, "Bayesian approach to parameter estimation of the generalized Pareto distribution," Test, vol. 12, no. l,pp. 259-277,2003.

[24] S. F. Juarez and W. R. Schucany, "Robust and efficient estimation for the generalized Pareto distribution," Extremes, vol. 7, no. 3, pp. 237251,2004.

[25] J. Zhang, "Likelihood moment estimation for the generalized Pareto distribution," Aust. N. Z. J, Stat., vol. 49, no. 1, pp. 69-77, 2007.

[26] J. Zhang and M. Stephens, "A new and efficient estimation method for the generalized Pareto distribution," Technometrics, vol. 51, no. 3, pp. 316-325,2009.

[27] P. de Zea Bermudez and S. Kotz, "Parameter estimation of the generalized Pareto distribution-Part II," J. Stat. Plan. Ir.ference, vol. 140, no. 6, pp. 1374-1388, 2010.

[28] J. Song and S. Song, "A quantile estimation for massive data with generalized Pareto distribution," Comput. Stat. Data Anal., vol. 56, no. l,pp. 143-150,2012.

[29] C. Wang and G. Chen, "A new hybrid estimation method for the generalized Pareto distribution," Comm. Stat. Theory Methods, vol. 45, no. 14, pp. 4285-4294, 2016.

[30] R. M. El-Sagheer, T. M. Jawa and N. Sayed-Ahmed, "Inferences for generalized Pareto distribution based on progressive first-failure censoring scheme," Complexity, vol. 2021, 2021, Art. No. 9325928.

[31] S. G. From and S. Ratnasingam, "Some efficient closed-form estimators of the parameters of the generalized Pareto distribution," Environ. Ecok Stat.,vol. 29, pp. 827-847, 2022.

[32] J. Pickands, "Statistical inference using extreme order statistics," The Ann. Stat., vol. 3, no. l,pp. 119-131, 1975.

[33] A. Balkema and L. de Haan, "Residual life time at great age," Ann. Probab., vol. 2, no. 5, pp. 792-804, 1974.

[34] R. L. Smith, "Threshold Methods for Sample Extremes," in Statistical Extremes and Applications, J. T. de Oliveira (Ed.), Dordrecht: Springer, 1984, pp. 621-638.

[35] A. C. Davison, "Modelling Excesses over High Thresholds, with an Application," in Statistical Extremes and Applications, J. T. de Oliveira (Ed.), Dordrecht: Springer, 1984, pp. 461-482.

[36] R. L. Smith, "Maximum likelihood estimation in a class of nonregular cases," Biometrika, vol. 72, no. 1, pp. 67-90, 1985.

[37] H. Nagatsuka and N. Balakrishnan, "Efficient likelihood-based inference for generalized Pareto distribution," Ann. Instit. Stat. Math., vol. 73, no. 6, pp. 1153-1185, 2021.

[38] P. Worapun and A. Kangtunyakarn, "An approximation method for solving fixed points of general system of variational inequalities with convergence theorem and application," 1AENG Int. J. Appt Math., vol. 51, no. 3, pp. 751-756, Sep. 2021.

[39] Mohammed Barrada, Reda Benkhouya, Mohammed Lahmer, and Idriss Chana, "On the Global Convergence of a New Super Halley's Family for Solving Nonlinear Equations," Engineering Letters, vol. 29, no.4, ppl407-1413,2021.

[40] D. Shah, and M. Sahni, "DMS way of finding the optimum number of iterations for fixed point iteration method," Lecture Notes in Engineering and Computer Science: Proceedings cf the World Congress on Engineering 2018, 4-6 July, 2018, London, U.K., pp. 8789.

[41] P. C. Chang, Y. W. Wang and C. H. Liu, "The development of a weighted evolving fuzzy neural network for PCB sales forecasting," Expert Syst. Appk, vol. 32, no. 1, pp. 86-96, 2007.

[42] R. A. Ibrahim, Sukono, H. Napitupulu, R. I. Ibrahim, M. D. Johansyah and J. Saputra, "Estimating flood catastrophe bond prices using approximation method of the loss aggregate distribution: Evidence from Indonesia," Decis. Sei. Lett., vol. 12, no. 2, pp. 179-190, Mar. 2023.

[43] M. Gerson, "The techniques and uses of probability plotting," J. R. Stat. Soc., vol. 24, no. 4, pp. 235-257, 1975.

[44] A. Luceño, "Fitting the generalized Pareto distribution to data using maximum goodness-of-fit estimators," Comput. Stat. Data Anal., vol. 51,no.2,pp. 904-917,2005.

[45] M. Raschke, "Opportunities of the minimum Anderson-Darling estimator as a variant of the maximum likelihood method," Commun. Stat. Simul. Comput, vol. 46, no. 9, pp. 6879-6888, 2017.

[46] Sukono, R.A. Ibrahim, M. P. A. Saputra, Y. Hidayat, H. Juahir, I. G. Prihanto and N.B.A. Halim, "Modeling multiple-event catastrophe bond prices involving the trigger event correlation, interest, and inflation rates," Mathematics, vol. 10, no. 24, pp. 4685, Dec. 2022.

[47] R. A. Ibrahim, Sukono and H. Napitupulu, "Multiple-trigger catastrophe bond pricing model and its simulation using numerical methods," Mathematics, vol. 10, no. 9, pp. 1363, Apr. 2022.

[48] Z. G. Ma and C.Q. Ma, "Pricing catastrophe risk bonds: A mixed approximation method," Insur. Math. Econ., vol. 52, no. 2, pp. 243254,2013.

[49] R. A. Ibrahim, Sukono, H. Napitupulu and R. I. Ibrahim, "How to price catastrophe bond for sustainable earthquake funding? A systematic review of the pricing framework," Sustainability, vol. 15, no. 9,pp. 7705, May. 2023.