An efficient estimator of the parameters of the Generalized Lambda Distribution

TL;DR

This paper introduces a fast, two-step estimator for the four parameters of the Generalized Lambda Distribution, improving computational efficiency while maintaining good performance, enabling easier interval estimation via bootstrapping.

Contribution

A novel, simple two-step estimator for the Generalized Lambda Distribution parameters that is faster and performs well compared to existing methods.

Findings

Estimator is significantly faster than existing methods.

Performs well in simulations and real data applications.

Enables effective bootstrap-based interval estimation.

Abstract

Estimation of the four generalized lambda distribution parameters is not straightforward, and available estimators that perform best have large computation times. In this paper, we introduce a simple two-step estimator of the parameters that is comparatively very quick to compute and performs well when compared with other methods. This computational efficiency makes the use of bootstrapping to obtain interval estimators for the parameters possible. Simulations are used to assess the performance of the new estimators and applications to several data sets are included.

Tables

Table 1. Table 5: Parameters chosen for the Generalized Lambda Distribution using the pdQ method and quartile matching.

Distribution	${\lambda }_1$	${\lambda }_2$	${\lambda }_3$	${\lambda }_4$
Normal (0,1)	0	1.4420	0.1469	0.1469
Lognormal (0,1)	0.8038	1.8141	0.7589	-0.7082
${\chi }_5^2$	4.0559	0.4977	0.5167	-0.1470
Beta (2,3)	0.3770	5.3836	0.4958	0.2637

Table 2. Table 6: This table shows the coverage probability (cp) and mean width ( ω ¯ ¯ 𝜔 \bar{\omega} ) of bootstrap confidence intervals for the location ( λ 1 subscript 𝜆 1 \lambda_{1} ) of the closest GLD distribution at nominal level 95%. These values are obtained using 500 bootstrap re samples from and 500 iterations

n	F		pdQ		TM		PM
			perc	bca	perc	bca	perc	bca
100	Normal	cp	0.966	0.966	0.964	0.960	0.966	0.942
		$\overline{\omega }$	0.560	0.561	0.489	0.490	0.698	0.701
	${\chi }_5^2$	cp	0.974	0.968	0.946	0.956	0.978	0.918
		$\overline{\omega }$	1.735	1.734	1.500	1.453	2.072	2.035
	Lognormal	cp	0.988	0.960	0.950	0.950	0.968	0.882
		$\overline{\omega }$	0.578	0.572	0.493	0.471	0.682	0.660
	Beta	cp	0.968	0.964	0.948	0.960	0.956	0.940
		$\overline{\omega }$	0.142	0.142	0.123	0.121	0.174	0.173
250	Normal	cp	0.948	0.948	0.960	0.966	0.942	0.948
		$\overline{\omega }$	0.341	0.342	0.290	0.291	0.394	0.397
	${\chi }_5^2$	cp	0.974	0.966	0.978	0.974	0.960	0.944
		$\overline{\omega }$	1.010	1.011	0.822	0.824	1.362	1.350
	Lognormal	cp	0.972	0.964	0.956	0.950	0.962	0.928
		$\overline{\omega }$	0.347	0.363	0.282	0.270	0.498	0.484
	Beta	cp	0.950	0.934	0.934	0.940	0.944	0.942
		$\overline{\omega }$	0.083	0.083	0.067	0.067	0.107	0.107

Table 3. Table 7: This table shows the coverage probability (cp) and mean width ( ω ¯ ¯ 𝜔 \bar{\omega} ) of bootstrap confidence intervals for the λ 3 − λ 4 subscript 𝜆 3 subscript 𝜆 4 \lambda_{3}-\lambda_{4} of the closest GLD distribution at nominal level 95%. These values are obtained using 1000 and 2000 bootstrap re samples subsequently for n=500 and n=1000 and 500 iterations

n	F		pdQ		PM
			perc	bca	perc	bca
500	Normal	cp	0.972	0.964	0.952	0.918
		$\overline{\omega }$	0.251	0.253	0.457	0.465
	${\chi }_5^2$	cp	0.950	0.962	0.956	0.934
		$\overline{\omega }$	0.311	0.297	0.743	0.747
	Lognormal	cp	0.876	0.942	0.978	0.918
		$\overline{\omega }$	0.560	0.477	1.541	1.450
	Beta	cp	0.978	0.972	0.972	0.952
		$\overline{\omega }$	0.295	0.300	0.646	0.650
1000	Normal	cp	0.940	0.942	0.958	0.934
		$\overline{\omega }$	0.167	0.167	0.305	0.304
	${\chi }_5^2$	cp	0.944	0.956	0.964	0.940
		$\overline{\omega }$	0.194	0.193	0.415	0.415
	Lognormal	cp	0.894	0.934	0.956	0.932
		$\overline{\omega }$	0.304	0.293	1.012	0.996
	Beta	cp	0.976	0.968	0.948	0.940
		$\overline{\omega }$	0.172	0.172	0.383	0.383

Table 4. Table 8: The KS test statistic for Spanish household data

pdQ	TL	PM	ML	MPS	TM	SM	DLA
0.0069	0.0043	0.0091	0.0326	0.0326	0.0326	0.0183	0.0155

Table 5. Table 9: Running times of each method for spanish household data

Method	Elapsed time (seconds)	Relative to pdQ
PM	0.02	0.05
pdQ	0.41	1
TL	2.13	5.20
ML	11.62	28.34
TM	17.19	41.93
MPS	17.22	42
DLA	21.03	51.29
SM	37.39	91.20

Table 6. Table 10: The KS test statistic and p-value for twin study data

	pdQ	TL	PM	ML	MPS	TM	SM	DLA
D	0.0465	0.0447	0.0583	0.0519	0.0487	0.0487	0.0417	0.0494
p-value	0.9501	0.9657	0.7973	0.8993	0.9344	0.9284	0.9791	0.9236

Table 7. Table 11: Point and interval estimates for the difference of location for earnings of males and females from 1992 to 1998 in US

		pdQ	TM	PM
Est.		2.276	2.282	2.288
Percentile	CI	(2.026, 2.674)	(2.011, 2.544)	(1.989, 2.736)
	Time(mins)	0.945	213	0.405
BCa	CI	( 1.968, 2.604 )	(2.021, 2.545)	(1.929, 2.631)
	Time(mins)	59.17	6451.2	45.03

Table 8. Table 12: Point and interval estimates for the difference of l3 and l4 for earnings of males from 1992 to 1998 in US

	pdQ	TM	PM
Est.	0.383	0.351	0.417
CI-Perc	(0.365, 0.398)	(0.328, 0.370)	(0.385, 0.459)
CI-BCa	(0.345, 0.418)	(0.329, 0.373)	(0.353, 0.496)

Table 9. Table 13: Point and interval estimates for the difference of l3 and l4 for earnings of females from 1992 to 1998 in US

	pdQ	TM	PM
Est.	0.486	0.410	0.524
CI-Perc	(0.452, 0.529)	(0.385, 0.436)	(0.455, 0.618)
CI-BCa	(0.447, 0.524)	(0.385, 0.437)	(0.440, 0.602)

Equations

$$\hat{q}_{n}(u)=\sum_{i=1}^{n}X_{(i)}\bigg{\{}k_{b}\bigg{(}u-\frac{(i-1)}{n}\bigg{)}-k_{b}\bigg{(}u-\frac{i}{n}\bigg{)}\bigg{\}}$$

$${\hat{f}}_Q^{\ast }(u_j)=\frac{1}{\hat{\kappa }{\hat{q}}_n(u_j)}\text{where}\hat{\kappa }=\frac{1}{J}\sum _{j=1}^J\frac{1}{{\hat{q}}_n(u_j)}.$$

$$Q(u)={\lambda }_1+\frac{1}{{\lambda }_2}\left[\frac{u^{{\lambda }_3}-1}{{\lambda }_3}-\frac{(1-u{)}^{{\lambda }_4}-1}{{\lambda }_4}\right]$$

$$f_Q(u)=\frac{{\lambda }_2}{u^{{\lambda }_3-1}+(1-u{)}^{{\lambda }_4-1}}.$$

$$({\hat{\lambda }}_3,{\hat{\lambda }}_4)=\underset{{\lambda }_3,{\lambda }_4}{\mathrm{argmin}}\sum _{j=1}^J{\left[{\hat{f}}_Q^{\ast }(u_j)-f_Q^{\ast }(u_j;{\lambda }_3,{\lambda }_4)\right]}^2.$$

$${\hat{x}}_p={\lambda }_1+\frac{1}{{\lambda }_2}c(p,{\hat{\lambda }}_3,{\hat{\lambda }}_4)$$

$$c(p,{\hat{\lambda }}_3,{\hat{\lambda }}_4)=\left[\frac{p^{{\hat{\lambda }}_3}-1}{{\hat{\lambda }}_3}-\frac{(1-p{)}^{{\hat{\lambda }}_4}-1}{{\hat{\lambda }}_4}\right]$$

$${\hat{\lambda }}_2=\frac{c(0.75,{\hat{\lambda }}_3,{\hat{\lambda }}_4)-c(0.25,{\hat{\lambda }}_3,{\hat{\lambda }}_4)}{{\hat{x}}_{0.75}-{\hat{x}}_{0.25}},{\hat{\lambda }}_1={\hat{x}}_{0.5}-\frac{c(0.5,{\hat{\lambda }}_3,{\hat{\lambda }}_4)}{{\hat{\lambda }}_2}$$

Full text

An efficient estimator of the parameters of the Generalized Lambda Distribution

Dilanka S. Dedduwakumara a, Luke A. Prendergast a and Robert G. Staudte a

aDepartment of Mathematics and Statistics, La Trobe University

ABSTRACT. Estimation of the four generalized lambda distribution parameters is not straightforward, and available estimators that perform best have large computation times. In this paper, we introduce a simple two-step estimator of the parameters that is comparatively very quick to compute and performs well when compared with other methods. This computational efficiency makes the use of bootstrapping to obtain interval estimators for the parameters possible. Simulations are used to assess the performance of the new estimators and applications to several data sets are included.

Key words: bootstrap interval estimator, generalized lambda distribution, probability density quantile

1 Introduction

The generalized lambda distribution (GLD, Ramberg & Schmeiser,, 1974) is a flexible, four parameter distribution that can approximate a large variety of distributions of varying shapes. Due to this flexibility, the GLD has become a popular distribution to model data in many fields, including economics and finance (e.g. Pfaff,, 2016). Given that it has four parameters, location, scale and two shape parameters, estimation is not a trivial task, and estimation methods continue to attract attention in the literature. Dedduwakumara *et al. *, (2019) introduced a method for choosing optimal parameters for generalized distributions to approximate other distributions. The method uses the probability density quantile function (pdQ, Staudte,, 2017) to first find the optimal shape parameters and is computationally quick, simple to implement and often outperforms other methods. Motivated by this, we introduce estimators of the GLD shape parameters arising from the estimated pdQ. The estimators, which compare favourably to other methods, is computationally efficient, providing estimates in just a fraction of the time required for other methods that are good estimators of the parameters. This efficiency makes obtaining interval estimators for the GLD parameters via bootstrapping possible, which is a considerable advantage given the lack of intervals in the literature.

We begin by providing important definitions and notations in Section 2 before reviewing several popular estimators already available in Section 3. The new estimator is introduced in Section 4 and performance is assessed by simulations in Section 5. Several data applications are given in Section 6 before we provide concluding remarks in Section 7.

2 Definitions

Throughout, let $Q(u)=F^{-1}(u)$ denote the quantile function associated with distribution function $F$ for $u\in(0,\ 1)$, and let $f=F^{\prime}$ denote the probability density function where $f(x)>0$ for all $x$ in its domain.

2.1 Probability density quantile functions

The quantile density function (Parzen,, 1979), also called the sparsity index by Tukey, (1965), is denoted $q(u)=Q^{\prime}(u)=1/f\left[Q(u)\right]$. The quantile density function is mainly used in non-parametric modeling and inference. Parzen, (1979) called the reciprocal of this quantile density function the density quantile function which we denote here by $f_{Q}(u)=f[Q(u)].$ In a recent paper, Staudte, (2017) introduced the probability density quantile (pdQ) denoted $f^{*}_{Q}(u)=f_{Q}(u)/\kappa$ where $\kappa=E[f_{Q}(U)]$ and $U\sim\text{Unif}(0,1)$. This pdQ function is free from location and scale and therefore can easily be used to examine shape behaviours of a distribution; see also Staudte & Xia, (2018) who provide other insights for the pdQ. Moreover, it is defined on the finite domain $[0,1]$ for all lattice distributions and continuous distributions having square-integrable densities.

2.2 Empirical pdQ for continuous distributions

Let $X_{1},\ldots,X_{n}$ denote a random sample of size $n$ from $F$ and let $X_{(1)}\leq X_{(2)}\leq\ldots\leq X_{(n)}$ denote the ordered sample. In this section we introduce the empirical pdQ earlier discussed by Staudte, (2017). For $k_{b}(\textperiodcentered)=k(\cdot/b)/b$ denoting a kernel function and $b$ a bandwidth, we start by estimating $q(u)$ using the quantile density estimator

[TABLE]

which consists of a linear combination of order statistics. This kernel density estimator has been studied extensively, e.g. see Jones, (1992), Falk, (1986) and Welsh, (1988) for some notable works. The choice of bandwidth is important and we choose our bandwidth to be $b(u)=(15/n)^{1/5}\{q(u)/q^{{}^{\prime\prime}}(u)\}^{2/5}$ since it minimizes the asymptotic mean squared error of $\hat{q}(u)$. Prendergast & Staudte, (2016) call $q(u)/q^{{}^{\prime\prime}}(u)$ the quantile optimality ratio (QOR) and estimate it to obtain a suitable bandwidth.

By using a discrete set of $u$s defined by $\{u_{j}=(j-1/2)/J\}^{J}_{j=1}$ and for some integer $J$, the empirical pdQ can be defined as

[TABLE]

2.3 The Generalized Lambda Distribution

The generalized lambda distribution (GLD) is a very flexible distribution that can approximate, or is equal to exactly, many other distributions for appropriately chosen parameters. Although several definitions for the GLD are available, we prefer the parameterization by Freimer *et al. *, (1988) since it is defined for all parameter value choices, with the exception that the scale parameter must be positive. The distribution is defined in terms of its quantile function which is

[TABLE]

where $\lambda_{1}$ is a location parameter, $\lambda_{2}>0$ an inverse scale parameter and $\lambda_{3},\lambda_{4}$ are shape parameters. It is easy to see that the quantile density function for the GLD is $q(u)=Q^{\prime}(u)=\lambda_{2}^{-1}\left[u^{\lambda_{3}-1}+(1-u)^{\lambda_{4}-1}\right]$ so that the density quantile function is

[TABLE]

In general, no closed-form solution for the integral $\kappa=\int^{1}_{0}f_{Q}(u)du$ exists, but it can be evaluated computationally and quite efficiently since integration occurs only between the finite bounds zero and one.

3 Existing Methods

Having four parameters, GLD estimation is not straightforward, and several estimation methods are available in the literature. For comparison with the pdQ method, we consider the more common estimation methods that are believed to provide good estimates of the true parameter values. These methods are available in either or both of the gld (King *et al. *,, 2016) and bda packages (Wang,, 2015) in the R statistical software (R Core Team,, 2017). We briefly describe those estimation methods and an overview with more technical details can be found in Dean, (2013).

The trimmed L-moments method (TL) is discussed in Asquith, (2007) and Dean, (2013) and is a robust version of the matching L-moments method (Karvanen *et al. *,, 2002). L-moments are estimated by linear combinations of order statistics and are related to conventional moments. Trimmed L-moments (see: Elamir & Seheult, (2003) for more details) are generalizations of these L-moments which applies zero weight to extreme observations. GLD parameters are estimated through this method by minimizing the difference between the sample trimmed L-moments of the data and the trimmed L-moments of the fitted GLD distribution.

The percentile matching method (PM) equates a selected number of empirical percentiles with their GLD counterparts to obtain a set of non-linear equations which are then solved to obtain the optimal GLD parameters. The PM method is discussed by Karian & Dudewicz, (1999) and Tarsitano, (2005). Further, Karian & Dudewicz, (2003) demonstrate the superiority of the PM method compared to the method of matching moments and the method of L-moments.

Maximizing the likelihood (ML) to obtain the GLD parameters has also been a popular method and is considered by Su *et al. *, (2007) and Su, (2007). First, initial estimated parameter values using the method of moments or the PM method are chosen and used as a starting point to seek the values that maximize the numerical log likelihood.

There are instances where the Maximum Likelihood method fails as the support depends upon the parameters to be estimated. In such a scenario, the Maximum product of spacings (MPS) was introduced by Cheng & Amin, (1983); Ranneby, (1984) and used by Chalabi *et al. *, (2012) to estimate GLD parameters. Here spacings refer to the differences between the cumulative distribution function at neighbouring data points and the parameters are estimated by maximizing the geometric mean of these spacings.

Rather than the spacing between transformed data points as in MPS, Titterington, (1985) suggested spacing between transformed, adjacently averaged data points. This approach is called the Titterington Method (TM), and more details with the GLD estimation can be found in Dean, (2013).

Owen, (1988) introduces the Starship Method (SM), a computer-intensive method that focuses on the fit to the base distribution of the inversely transformed data. King & MacGillivray, (1999) developed this concept to be used for generalized lambda distribution where for a given data set the distribution function of GLD is obtained numerically, and the parameter values are chosen such that it minimizes the goodness-of-fit to the uniform distribution. One major drawback of this method is its slow computation time, especially with large sample sizes.

The method of distributional least absolutes (DLA) obtains the optimal GLD parameters by minimizing the sum of absolute deviations between the order statistics and the corresponding medians. See Dean, (2013) for more details.

4 Point and interval estimators using the pdQ

Obtaining optimal GLD parameters for a specific distribution using the pdQ method has been discussed in Dedduwakumara *et al. *, (2019) and we adopt this idea to model empirical data as detailed in the steps below. Throughout this section the estimated $p^{th}$ quantile, the empirical pdQ of the data and the corresponding GLD pdQ are denoted by $\widehat{x}_{p}$, $\hat{f}_{Q}^{*}$ and $f_{Q}^{*}$ respectively.

4.1 Point estimators

Step 1: Estimating the shape parameters

Since the pdQ is free of location and scale, we can use it to focus only on estimating the two shape parameters, $\lambda_{3}$ and $\lambda_{4}$. Using a discrete set of probability values given as $\{u_{j}=(j-1/2)/J\}^{J}_{j=1}$, the estimated shape parameters are those that minimize the sum of squared distances between the empircal pfdQ and GLD pdQ. That is,

[TABLE]

To find the shape parameters that minimize the sum of squared distances, it is simple to use a computational optimizer. We use the R function nlminb in stats package for this, which we find to be both simple to use and quick to compute. As starting values for the optimizer we chose to apply the objective function in (5) to a grid of values for $\lambda_{3}$ and $\lambda_{4}$ consisting of each paired combination of $\lambda_{3},\lambda_{4}$ from $−0.9,−0.5,−0.1,0,0.1,0.2,0.4,0.8,1,1.5$ and choose those values that result in the smallest value. This grid of values are used in Dean, (2013) and considered as the default choice in gld (King *et al. *,, 2016) package for other GLD estimation methods. This covers a wide span of the $\lambda_{3},\lambda_{4}$ values and also we can expect uniformity between methods for the comparisons to follow. Our simulations, to be summarized in the next section, reveal that $J=50$ is more than adequate to obtain good estimates for moderate to larger sample sizes and there is little value in increasing $J$ to be more than this. Therefore, in what follows, we use $J=50$ for sample sizes greater than 200. We found that $J=25$ worked well for smaller sample sizes.

Step 2: Estimating the location and scale parameters

Given the estimated shape parameters from Step 1, we then match sample quartiles $\widehat{x}_{p}$, using $p\in\{0.25,0.5,0.75\}$, to their theoretical GLD counterparts to generate the linear equations

[TABLE]

where,

[TABLE]

By solving the the system of three linear equations from (6) using $p=0.25,0.5$ and $0.75$, we obtain our estimated inverse scale and location parameters as

[TABLE]

respectively.

4.2 Interval estimators

In order to calculate confidence interval estimates, we consider two commonly used bootstrap procedures; the percentile method (Efron & Tibshirani,, 1994) and the bias-corrected and accelerated (BCa) bootstrap intervals (Efron,, 1987). As the pdQ method requires little computational time, the computational cost to calculate the bootstrap confidence intervals is not prohibitive.

The percentile bootstrap method uses the upper and lower $\alpha/2$ percentiles of the GLD sample estimators as the interval bounds, whereas the BCa method also estimates a bias-corrected parameter and an acceleration parameter. The bias-correction parameter is related to the proportion of bootstrap estimates that are less than the observed statistic, and the acceleration parameter is proportional to the skewness of the bootstrap distribution. Unlike the percentile method, the BCa method corrects for both bias and skewness of the estimators.

In the following section to come, we look at the confidence intervals for the location ($\lambda_{1}$) and the skewness ($\lambda_{3}-\lambda_{4}$) of the underlying GLD distribution by using these bootstrap interval estimators. Note that the difference $\lambda_{3}-\lambda_{4}$ may be of interest since $\lambda_{3}-\lambda_{4}=0$ indicates that the distribution is symmetric.

5 Simulation results

In this section, we compare our pdQ estimation method with the existing methods presented in Section 3. To do so, we summarize the results for 500 simulated data sets by reporting the mean standard error and the mean absolute bias for estimates of the four GLD parameters compared to their true values. We adopt four representative FMKL GLD settings used in Corlu & Meterelliyoz, (2016) to evaluate the performances of each method for these different shapes of GLD. The shapes of the probability density functions for these GLD settings are shown below in Figure 1.

Results are evaluated for a range of sample sizes: 100, 250, 500 and 1000. In the pdQ method, for sample size 100, we take $J=25$, and for other sample sizes, we use $J=50$.

As can be seen in Table 1, the pdQ method shows comparatively low mean squared error and absolute bias for the $\lambda_{1}=0,\lambda_{2}=1,\lambda_{3}=1.5,\lambda_{4}=1.5$ setting for all the sample sizes. The MPS and TM methods also perform well, being second best or third best in most instances. Similar results can be found in Table 2 for $\lambda_{1}=0,\lambda_{2}=1,\lambda_{3}=2.5,\lambda_{4}=1.5$ as well, where the pdQ method has the minimum MSE and bias values in almost all cases. The DLA and MPS methods also perform well for this setting. Although the pdQ is not the best method always for the last two GLD settings, $\lambda_{1}=0,\lambda_{2}=1,\lambda_{3}=2,\lambda_{4}=0.5$ and $\lambda_{1}=0,\lambda_{2}=1,\lambda_{3}=0.5,\lambda_{4}=0.6$, it produces competitive results compared to the other methods. The MPS and TM methods again provide very good results for these two settings, displaying low MSE and absolute bias values (Table 3 and 4). However, their good performance comes with much greater computational cost which we will consider below.

In Figure 2, we depict the computational time for each method as the sample size increases. We consider 100 replications of data from GLD(0,1,0.5,0.6) for each sample size and calculate the average computational time for each estimation method. The run times were measured on a computer with an Intel(R) Core(TM) i7-6700 processor running at 3.40GHz using 32 GB of RAM, running Windows version 10. The graph shows that the PM and pdQ methods have much lower computational times compared to all other methods. The PM takes just a few milliseconds to estimate the parameters from a sample even as large as 100,000 observations. However, the PM method was typically not a good estimator when compared to the pdQ, MPS and TM methods. On the other hand, the pdQ method requires about a second to compute even for 100,000 observations. As expected, SM and ML are the slowest among all the methods, spending almost 4 minutes each to calculate the above results. The TM method has a computation time of around 2.5 minutes for this setting. Such large computation times make bootstrapping onerous.

We now consider the performance of bootstrap intervals for some common distributions. The closest GLD parameters for these considered distributions are adopted from Dedduwakumara *et al. *, (2019) and are presented in Table 5. Here one sample is generated from the above closest GLD parameters and bootstrapped from those observations.

In Table 6 and Table 7 we provide the coverage probability results for bootstrap confidence intervals for both location ($\lambda_{1}$) and the difference of shape parameters ($\lambda_{3}-\lambda_{4}$), where $\lambda_{3}-\lambda_{4}$ is an indication of the skewness of the GLD distribution. We only selected the TM, PM methods together with PdQ methods for the comparison as they are the close competitors in performance wise and computational time wise.

From Table 6, all the methods display coverage probabilities close to the nominal value for the location parameter. The TM method has the narrower width, with the pdQ method the next narrowest. We found that larger sample sizes are required to obtain reliable intervals for the difference in shape parameters. The TM method is not considered here as BCa bootstrap confidence intervals need a higher number of re-samples and when the sample size increases to 1000, and it takes roughly 33.98 minutes to obtain a single confidence interval for TM. The pdQ and PM methods, on the other hand, take only 47.65 seconds and 32.02 seconds subsequently to provide a single confidence interval for $n=1000$. Compared to the PM method, the pdQ method shows better coverage in the vicinity of the nominal level for skewness with a narrower width. So it can be concluded that the pdQ is favourable compared to these other two methods.

6 Applications

6.1 Application 1: Total income data of private households of Spain in 1980

In this example, we consider total income data of Spanish households (Figure 3), which is from the 1980 Spanish Family Expenditure Survey (FES) described in Alonso-Colmenares *et al. *, (1994). This data set consists of 23,972 observations, and total income is recorded with household characteristics and expenditure on several categories. This data set is readily available in the Ecdat package (Croissant,, 2016) under the name ‘BudgetFood’.

In Table 8, we assess the goodness of fit of each method to this data by using the Kolmogorov Smirnov test statistic. As can be seen, both TL-moments and pdQ method provide the best fit to the data. Further, this is elaborated in Figure 4 quantile plots where, unlike the other methods, estimated quantiles from the pdQ method and TL-moments method are clearly aligned with the sample quantiles.

In Table 9, we provide the running times for parameter estimation by each method. The PM and pdQ methods are clearly the fastest. The pdQ is about five times faster than the TL-moments method and 91 times faster than SM for this particular data set. This further suggests the suitability of the pdQ method with large sample sizes.

6.2 Application 2: Twin Study data

We now consider a much smaller data set known as the Indiana Twin Study data, which consists of the birth weights of a set of 123 twins. The data originally comes from the PhD thesis of Dr Cynthia Moore, Department of Medical and Molecular Genetics, Indiana University School of Medicine and also has been used in Karian & Dudewicz, (2003) and Karian & Dudewicz, (2000), for GLD parameter estimation.

Table 10 presents the Kolmogorov Smirnov test statistic for the GLD fit to the data by each method and the corresponding p-value for the test. According to the p-values, all the methods suggest a better fit to the data, but in particular, the SM method, TL-moments method and pdQ method are the best fits, as they display smaller KS statistics compared to the other methods.

6.3 Application 3: Earnings data

We now look at obtaining 95% bootstrap confidence intervals for the difference of the location parameter for two groups and also confidence intervals for the skewness. The data set we consider here is also available in the Ecdat package (Croissant,, 2016) and is named ‘CPSch3’. The data set includes the hourly earnings of males ($n=5956$) and females ($n=5174$) in the US from 1992 to 1998.

In Tables 11,12 and 13 we display the bootstrap confidence intervals calculated using three methods and the time efficiency of each method. For percentile bootstrap intervals, we used 500 bootstrap samples, but for BCa confidence intervals, the number of re-samples has to be much higher than the sample size. Therefore, we used 15,000 bootstrap re-samples to obtain the BCa confidence intervals as the sample sizes here are large.

All the methods provide similar results for the confidence intervals with narrower width while declaring there is a difference in the location for both males and females in earnings. Additionally, there is some skewness present in both the male and female earnings distributions. When the width of the confidence interval is considered, pdQ has a narrower width after the TM method and outperforms the PM method. When the time efficiency is considered the pdQ and PM take only 0.945 minutes and 0.405 minutes, respectively, to obtain the percentile confidence interval whereas the TM method spends around 3.5 hours to obtain the same results. For the BCa confidence intervals, TM method is much more inefficient in this case as it takes roughly 4.5 days to obtain the BCa confidence interval for this data whereas the pdQ method and PM method only takes 59.17 minutes and 45.03 minutes respectively.

7 Discussion

In this paper we have introduced a simple two-step estimation method for the GLD parameters. The first step uses the empirical probability density quantile function to find estimates of the shape parameters. The second step can then simply obtain estimates for the location and inverse scale. How simulations show that the method performs very well, often beating existing methods. An additional advantage of the approach is the very small computation time compared to the best of the methods. This small computation time makes bootstrapping possible which is important given the lack of interval estimators available for the GLD parameters. For the more time consuming estimation methods, bootstrap intervals can take hours, even days to compute, as our examples show.

While we focused on the GLD distribution in this paper, the method can be used for other distributions too. An obvious extension would be to consider other generalized distributions, like the generalized Beta distribution, since estimation for such distributions is not straightforward.