Improved Shape Parameter Estimation for the Three-Parameter Log-Logistic Distribution

The log-logistic distribution is widely used in di ﬀ erent ﬁ elds of study such as survival analysis, hydrology, insurance, and economics. Recently, Ahsanullah and Alzaatreh studied the best linear unbiased estimators for the location and the scale parameters of the three-parameter log-logistic model. The same authors also propose a shift-invariant Hill estimator for the unknown shape parameter. In this work, we propose a new estimation method for the shape parameter. We derive its nondegenerate asymptotic behaviour and analyse its ﬁ nite sample performance through a Monte Carlo simulation study. To have precise estimates, we present a method for selecting the threshold. To illustrate the improvement achieved, e ﬃ ciency comparisons are also provided.


Introduction
Over the last decades, the need to solve problems in a diversity of applied areas, such as finance, hydrology, insurance, or survival analysis, gave rise to many new statistical distributions in the literature. A significant effort has been made toward the generalization of some classic distributions. A common technique consists in adding a parameter to a well-known distribution. It is often observed that the introduction of such an extra parameter brings more flexibility to the class of distribution functions. Here, we are interested in the three-parameter or shifted log-logistic distribution which is obtained from the classic log-logistic model with the addition of a location parameter. Thus, a random variable X follows a three-parameter log-logistic distribution if its probability density function and distribution function (d.f.) are, respectively, given by The corresponding quantile function is The constants α > 0, μ ∈ R, and σ > 0 are the shape, location, and scale parameters, respectively. Once there is no restriction between the values of the parameters, they can vary freely in the parameter space. If α = 1, X has a location and scale beta prime distribution. When α > 1, the probability density function f is unimodal with mode at the value μ + σððα − 1Þ/ðα + 1ÞÞ 1/α . Since the log-logistic model is heavy tailed, with a tail index equal to α, the i-th moment of X is finite only when i < α. In the following, the notation X~llogistðα, μ, σÞ will be used whenever the random variable X has the distribution function in (2). If X~llogistðα, μ, σÞ, then ðX − μÞ/σ, with σ > 0, has a standard log-logistic distribution, llogistðα, 0, 1Þ, with d.f. given by Figure 1 illustrates the probability density function and cumulative distribution function, for selected values of the shape parameter of the standard log-logistic distribution. As noted in [1], the standard log-logistic model is characterized by the relation The three-parameter log-logistic distribution is also known as a Pareto type III distribution (see Arnold [2]). Also, the two-parameter log-logistic distribution is a member of Burr's type XII family of distributions [3] and is known as the Fisk [4] distribution in the economic literature.
The log-logistic distribution in (2) is closely related with the logistic distribution. More precisely, if X~llogistðα, μ, σÞ, then is a logistic random variable with probability density function The half-logistic model, that is, the absolute value of the standard logistic model is another important model in the scientific literature. Extensions of the log-logistic or halflogistic models can be found in Cordeiro et al. [5,6], Alizadeh et al. [7], Mohammad [8], Lemonte [9], and Shakhatreh [10] among others.
Considerable attention has been paid to the estimation of the model parameters of the log-logistic distribution. Although more attention has been paid to the two-parameter case, several estimation procedures for the three-parameter log-logistic model are already available in the literature. Balakrishnan et al. [1] derived the best linear unbiased estimators (BLUE) for the location and scale parameters of a three-parameter log-logistic model, with a known shape parameter. In practice, it is unrealistic to assume that the shape parameter α is known, and it should be estimated. More recently, Ahsanullah and Alzaatreh [11] considered again the BLUE for the location and scale parameters of the log-logistic model. The authors propose the estimation of the shape parameter with the reciprocal of a Hill-type estimator applied to a tail sample fraction, shifted by the sample minimum. Moreover, Ahsanullah and Alzaatreh [11] proposed a sample fraction of 10%, if the sample size is greater than 100. However, such a simple suggestion has no theoretical or empirical support. The main objective of this paper is to improve the estimation of the shape parameter α of the log-logistic model. The motivation comes from the fact that the Hill estimator is biased and the sample fraction proposed in [11] is nonoptimal.
The paper is organized as follows. Section 2 describes the estimator proposed [11] and presents two alternative estimators and their asymptotic properties. It is shown that the estimators are asymptotically normal distributed and the choice of the sample fraction is discussed. In Section 3, we introduce a simple threshold selection method. The results of a Monte Carlo simulation that evaluates the mean value, the median, the standard deviation, and the root mean squared error of the estimators under study are reported in Section 4. Finally, Section 5 concludes the paper.

Estimator of the Shape Parameter and Their Properties
In what follows, we consider the estimation of the shape parameter of the log-logistic model in (2). We shall assume that ðX 1 , X 2 , ⋯, X n Þ is a sample of size n of independent and identically distributed random variables, with a common d.f. F, given in (2). All three parameters are assumed unknown. The corresponding sample of nondecreasing order statistics is denoted by ðX 1:n , X 2:n , ⋯, X n:n Þ.   [11] noted that the log-logistic distribution has a Pareto-type tail behaviour. More precisely, we can write with α the tail (or Pareto) index and lðxÞ = ð1 + x −α Þ −1 , x > 0, a slowly varying function at infinity. Moreover, since lðxÞ admits the Taylor expansion lðxÞ = 1 − x −α + oð−x α Þ, as x ⟶ ∞, the standard log-logistic model belongs to Hall's subclass (see equation (1) of [12]) of Pareto-type models with survival function, with ξ > 0, c > 0 is a first-order scale parameter, and ρ < 0 and β are second-order tail parameters. The second-order parameter ρ quantifies the deviation of the model to the Pareto distribution. For the standard log-logistic model, we have ξ = α −1 , ρ = −1, and c = β = 1. Therefore, F is in the max domain of attraction of the extreme value distribution with a positive shape parameter. This means that there exist  (i) Generate 5000 samples of size n from a log-logistic distribution, with n taking values between 50 and 2000, with step 50 (ii) Let b α j ðk, iÞ denote the estimates based on the i-th sample. For each sample size, compute b α j ðk, iÞ, k = 1, 2, ⋯, n − 2, i = 1, ⋯, 5000 (iii) Compute the empirical root mean squared error as a function of k (iv) Obtain the level k j ðnÞ that minimizes the empirical root mean squared error (v) Finally, perform a power regression with k j ðnÞ as the response variable and the sample size as the predictor variable. The regression coefficients are the vector values (a 1 , a 2 ).
3 Computational and Mathematical Methods sequences of normalizing constants a n > 0 and b n such that ðX n:n − b n Þ/a n converges in distribution to a nondegenerate random variable with distribution function with a positive shape parameter ξ. This parameter is the socalled extreme value index. In the statistic literature, one can find several estimators for ξ or, equivalently, for the shape parameter α. For a general overview of the available estimators, we refer to [13,14]. Whenever working with Pareto-type models, the Hill estimator [15] is frequently used to estimate the extreme value index. This estimator is defined as the average of the log excesses over the threshold u = X n−k:n > 0, where k represents the number of upper order statistics used in the estimation. For the strict Pareto model, with d.f. FðyÞ = 1 − ðy/cÞ −α , y > 0, (α > 0, c > 0), the Hill estimator is consistent and asymptotically unbiased for the estimation of α −1 . Moreover, if k = n − 1, HðkÞ is the maximum likelihood estimator of the reciprocal of the shape parameter α of the strict Pareto distribution. If F differs from the strict Pareto model, Also, the variance of HðkÞ decreases and the absolute bias increases, as k increases. Therefore, the choice of k leads to a trade-off between the bias and the variance of the estimator.
Many of the estimators that have been suggested to estimate a positive extreme value index, including the one in (10), are only scale-invariant. A change in the location can modify the asymptotic behaviour of the tail and the bias of a location-variant estimator (for more information, see the papers [16,17]). The properties of the Hill estimator and the fact that the d.f. in (2) has a location parameter, lead Ahsanullah and Alzaatreh [11] to propose the estimation Standard deviation Root mean squared error  Computational and Mathematical Methods of the shape parameter α with the following locationinvariant Hill-type estimator, The sample values are thus shifted by the sample minimum. In the following, we always add the notation * to any statistic based on shifted data. Note that the estimator H * ðkÞ in (12) is a member of the class of estimators in references [16,18,19]. Regarding the choice of the parameter k, Ahsanullah and Alzaatreh [11] proposed k = ½n/10, if n > 100, where ½x denotes the integer part of x. In order to improve the choice of the threshold k, it is crucial to gain information about the asymptotic behaviour of b α 1 . We begin with the following proposition that provides the distributional representation of b α 1 . The proof can be found in the appendix.

Proposition 1.
Assume that k is an intermediate sequence of integers satisfying (11). Then, the following distributional representation A typical approach in the literature is to choose the threshold through the minimization of the asymptotic mean squared error (AMSE) of b α 1 ðkÞ. Because such a choice depends on asymptotic arguments, it may only be reliable when the sample size becomes large. Alternative methods for selecting the threshold can be found in references [20][21][22][23]. From (13), it follows that the AMSE of Then, the level k 0 that minimizes the AMSE in Equation (15) is asymptotically equivalent to Furthermore, the level k ).

Alternative
Estimators. Despite its wide use to estimate the extreme value index, the Hill estimator is difficult to apply in real data problems due to the substantial bias. This problem motivated several researchers to reduce the bias of such an estimator and to construct alternative estimators. Reduced bias estimators of the extreme value index usually have a stable sample path, close to the target value, which makes them less sensitive to the choice of k. We mention the first reduced bias estimators of the extreme value index in references [24][25][26][27] (see also the papers [28,29] for a general overview on bias reduction). Bias reduction of tail parameter estimators typically requires the estimation of tail second-order parameters. However, in the present paper, we are assuming the parametric log-logistic distribution, a model within the class of heavy tailed models in (8), with ðρ, βÞ = ð−1, 1Þ. Therefore, in this paper, we shall consider two alternative estimators, which provide a reduction of the bias whenever ðρ, βÞ = ð−1, 1Þ. We shall consider the following alternative estimators: are the moments of order j of the shifted log-excesses. The statistic WH * is defined as Both statistics in (18) and (20) are data-shifted versions of statistics already considered in the literature. The nonshifted version of (18) was first introduced in [30], a particular case of a generalized Jackknife statistic based on the Hill and the moment ratio (MR) [31] estimators. Moreover, both HðkÞ and MRðkÞ are members of Lehmer mean-of-order-p class of extreme value index estimators studied in [32,33]. The nonshifted version of the estimator G * 1 ðkÞ was also independently introduced in [34] (see also [35]). The statistic in (20) is a location invariant modification of the class of weighted Hill estimators introduced in [36]. Next, we provide the asymptotic representation of the alternative estimators in (17). We do not provide proofs, since the arguments required to obtain the results are similar to the ones provided in the proof of Proposition 1 and in references [30,36].

Proposition 2.
Assume the conditions of Proposition 1. Then, the following distributional representation holds: where Z k is an asymptotically standard normal variable. Moreover, if ffiffi ffi k p ðk/2nÞ ⟶ c, then Regarding the asymptotic variances, note that we have Vðb α 1 ðkÞÞ < Vðb α 3 ðkÞÞ < Vðb α 2 ðkÞÞ. The absolute asymptotic bias of either b α 2 ðkÞ or b α 3 ðkÞ is null, while the absolute asymptotic bias of b α 1 ðkÞ is |Biasðb α 1 ðkÞÞ | = k/2n. Moreover, the results in (21) and (22) do not allow to evaluate the asymptotic optimal value of k. To find such an optimal value, it is necessary to obtain higher-order terms of the asymptotic expansion of bias. The major drawback of such an approach is that it can lead to inaccurate results for small sample sizes. For that reason, we shall address the choice of k in the next section.

A Threshold Selection Method
In practice, the threshold is fundamental to obtain accurate estimates, and it must be chosen before applying any of the aforementioned estimators of the parameter α. The naive method of selecting the top 10% of the data is extremely simple but has some drawbacks. Since such a threshold is not an intermediate sequence, it is not well supported by the theory. To have precise estimates for small sample sizes, we now propose a method for selecting the threshold that combines theoretical and empirical results. Hall [12] noted that the optimal performance of the Hill estimator is achieved if k = Oðδ 1 n δ 2 Þ for some δ 1 > 0 and δ 2 ∈ ½0, 1. Such a result also holds for other estimators of α, with a similar limit distribution. In the following, let b α j ðkÞ denote any of the aforementioned estimators of the shape parameter. We propose the following threshold selection procedure: with a 1 > 0 and a 2 ∈ ð0, 1Þ as suitable real numbers, chosen independently for each estimator. Additionally, it is worth noting that, at least for the estimator b α 1 ðkÞ, the constants a 1 and a 2 are asymptotically independent of any parameter of the log-logistic model in (2). We propose Algorithm 1 to obtain empirically both constants (a 1 , a 2 ) in (24).

Computational and Mathematical Methods
Remark 3. The constants in the first step of Algorithm 1 can be changed to modify the precision of the estimates of (a 1 , a 2 ).

Numeric Results
In this section, we use Monte Carlo simulation to study and compare the finite-sample properties of the estimators of the shape parameter α, in (12) and (17). For comparison purposes, we also include the reciprocal of the Hill estimator, b α H = 1/HðkÞ, with HðkÞ in (10). All the computations were done using R [37] software. We generated 5000 samples of size n from the log-logistic model with the same set of parameters considered in paper [11], that is, For the sample size n, we considered values between 50 and 5000. Let b α · ðkÞ, with ·∈f1, 2, 3, Hg, denote any of the aforementioned estimators under study. For each sample, the estimates of α are first computed for every k. We then and the simulated previous characteristics at the optimal level in (26). Furthermore, since the exact minimizer of the RMSE is not known, we also study the estimators with the level proposed in Section 3. This problem occurs only in a few generated samples. Based on the simulation results, b α 2 ðkÞ provides globally a good bias and RMSE patterns as a function of k. Therefore it constitutes an important estimation procedure for the shape parameter α.

4.2.
Results for the Optimal Level of k. We now compare the performance of the estimators at their simulated optimal level. The exact optimal value of k provides a benchmark of the best possible performance obtainable with each estimator of α. Notice that in practice, such an optimal level may not be achieved. For the sake of simplicity, we excluded the estimators b α H ðkÞ and b α 3 ðkÞ, the first due to the mixed performance, the second because it is less efficient than b α 2 ðkÞ. Thus, we only provide results for the estimator b α 1 ðkÞ, suggested in [11], and b α 2 ðkÞ, the estimator that provided a good performance for the four sets of the model parameters. The simulated values of the optimal sample fraction (the optimal level divided by the sample size), mean value, and RMSE, both computed at the simulated optimal threshold, are given in Table 1. Values associated to a smaller absolute bias and smaller RMSE are presented in bold.
Notice that the b α 2 ðkÞ has always an optimal level much higher than b α 1 ðkÞ. At their corresponding optimal level, b α 1 gives better results for small sample sizes (n < 500 if α = 0:5or 1.5, n < 100 if α = 2:5, and n < 50 if α = 4). For large values of n, b α 1 is less biased and provides a smaller RMSE.

Results for the Proposed Level k.
In this section, we study the distributional properties of the estimators, based on the tail data defined by the threshold proposed in (26). Again, we restrict ourselves to the estimators b α 1 ðkÞ and b α 2 ðkÞ in (12) and (17), respectively. We applied algorithm 3 to samples of the log-logistic model with parameters ðα, μ, σÞ = ð0:5,1, 1Þ. The algorithm provided the thresholds k 1 = 1:3n 0:66 Â Ã for b α 1 k ð Þ and k 2 = 1:1n 0:922 Â Ã for b α 2 k ð Þ: ð27Þ Figure 6 presents the empirical optimal sample fraction and the corresponding regression curve for both estimators. The overall agreement between empirical and fitted values is quite good.
Next, we provide numerical results of the finite sample performance of the estimators with the thresholds proposed in [11] and (27). More precisely, the following estimators are considered: From the simulation results obtained in Section 4.1, for every value of k, we obtained the mean value and the root mean squared error of the estimators under study. Results are presented in Table 2.
When looking at the empirical bias, b α * 0 achieves the smallest absolute bias if α < 4 and the sample size has a moderate size (usually between 50 and 500). Otherwise, b α * 0 has usually the smallest absolute bias. In terms of RMSE, the threshold approach in (27) compared favourably to the threshold proposed in [11]. If the sample size is small, b α * 1 provides the smallest RMSE. Otherwise, b α * 2 has the smallest RMSE. For all sets of parameters, the values of RMSE of α * 1 and α * 2 are very close to the corresponding optimal level in Table 1.

Conclusion
It is known that the Hill estimator can be seriously biased. To deal with this problem, we introduced two reduced bias estimators for the shape parameter of the log-logistic model and established its asymptotic normality. Additionally, we improved earlier guidelines for the choice of the threshold. Also, the presented simulation study shows improvements in the estimation of the shape parameter when compared to the classical estimation method in Ahsanullah and Alzaatreh [11]. These improvements are more pronounced if the sample size is large. We conclude by noting that the approach taken in this paper could be applicable to other bias-reduced estimators. Therefore, further research concerning alternative estimators of the shape parameter will be taken in the future. 10 Computational and Mathematical Methods