# A Doubling Method for the Generalized Lambda Distribution

Authors: Todd C. Headrick; Mohan D. Pant
Journal: ISRN Applied Mathematics
 (2012)

---

## Abstract

This paper introduces a new family of generalized lambda distributions (GLDs) based on a method ofdoubling symmetric GLDs. The focus of the development is in the context of L-moments and L-correlation theory. As such, included is the development of a procedure for specifying double GLDs with controlled degrees of L-skew, L-kurtosis, and L-correlations. The procedure can be applied in a variety of settings such as modeling events and Monte Carlo or simulation studies. Further, it is demonstrated that estimates of L-skew, L-kurtosis, and L-correlation are substantially superior to conventional product-moment estimates of skew, kurtosis, and Pearson correlation in terms of both relative bias and efficiency when heavy tailed distributions are of concern.


---

## Body

## 1. Introduction

The conventional moment-based family of generalized lambda distributions (GLDs) is often used in various applied mathematics contexts to model and describe data by a single functional form [1, page 5]. Some examples include modeling non-log-normal security price distributions [2], biological and physical phenomena [3], and solar radiation data [4]. The GLD is also a popular tool for generating random variables for Monte Carlo or simulation studies. Some examples include studies in such areas as operations research [5], microarray research [6], and structural equation modeling [7].The family of GLDs is based on the transformation [1, page 21], [9, 11], [10, page 127](1.1)q(u)=λ1+uλ3-(1-u)λ4λ2,
where u is uniformly distributed on the interval (0,1). The parameters λ1 and λ2 are location and scale parameters, while λ3 and λ4 are shape parameters that determine the skew and kurtosis of a GLD. The pdf and cdf associated with (1.1) can be expressed as [10, page 127](1.2)fq(u)(q(u))=f̅(u)=(q(u),1q′(u)),(1.3)Fq(u)(q(u))=F̅(u)=(q(u),u),
where f̅:ℜ↦ℜ2 and F̅:ℜ↦ℜ2 are parametric forms of the pdf and cdf with the mappings u↦(x,y) and u↦(x,v) with x=q(u), y=1/q′(u), v=u, and where 1 and u are the regular uniform pdf and cdf, respectively. It is assumed that q′(u)>0 in (1.2) to ensure a valid pdf that is, the transformation in (1.1) is strictly increasing. An essential requirement for a valid pdf is that λ2, λ3, λ4 in (1.1) all have the same sign [1, page 24]. For more specific details on the parameter space and conditions related to valid GLDs, see Karian and Dudewicz [1, pages 21–47]. Provided in Figure 1 is an example of a valid GLD pdf based on (1.1) and (1.2).Figure 1
A graph of a standardized GLD distribution (α1=0, α2=1) based on (1.1) and (1.2). The GLD has skew of α3=3 and kurtosis of α4=65. The solutions for λ1,…,4 are based on solving the equations in [1, page 57].Symmetric GLDs are produced for the case whereλ3=λ4 in (1.1) and where the mean (α1), variance (α22), skew (α3), and kurtosis (α4) can be determined from [8](1.4)α1=λ1,(1.5)α22=2/(1+2λ3)-2β[1+λ3,1+λ3]λ22,(1.6)α3=0,(1.7)α4=6β[1+2λ3,1+2λ3]-4β[1+λ3,1+3λ3]-4β[1+3λ3,1+λ3]+2/(1+4λ3)λ24.
Numerical solutions for λ2 and λ3 in (1.5) and (1.7) can be found in [1, Appendix B], which are associated with standardized GLDs (i.e., α1=λ1=0, α22=1). Note that the term β in (1.5) and (1.7) represents the complete beta function where the arguments cannot be negative. As such, for the kth moment to exist then we must have λ3>-1/k [8, 9, 11]. Thus, the condition λ3>-1/4 ensures that the first four moments exist.We propose a new family of asymmetric GLDs based on a technique referred to herein asdoubling symmetric GLDs. More specifically, a family of double GLDs can be created by selecting a pair of constants (λℒ, λℛ) and transforming separately for u≤1/2 and u≥1/2 as follows: (1.8)q(u)={(uλL-(1-u)λL)(δλL22-λL),for0<u≤12,(uλR-(1-u)λR)(δλR22-λR),for12≤u<1,
where λℒ (λℛ) is the nonzero shape parameter for the left (right) side of a distribution and δ is the (positive) parameter that determines the height of the double GLD at u=1/2. Provided in Figure 2 is an example of a double GLD pdf based on (1.2) and (1.8). Inspection of Figures 1 and 2 clearly indicates that these two GLDs are markedly different even though both distributions have the same values of skew (α3=3) and kurtosis (α4=65). Note that the values of λℒ and λℛ in Figure 2 were determined based on the equations for α3 and α4 given in the appendix.Figure 2
A graph of a double GLD distribution based on (1.2) and (1.8). The values of λℒ and λℛ were determined based on the equations for α3 and α4 in the appendix. The values of τ3 and τ4 were determined based on (2.10) and (2.11) in Section 2. The estimates (α̂3,4; τ̂3,4) and bootstrap confidence intervals (C.I.s) were based on resampling 25,000 statistics. Each sample statistic was based on a sample size of n=250.Conventional moment-based estimators (e.g.,α̂3,4) have unfavorable attributes insofar as they can be substantially biased, have high variance, or can be influenced by outliers. For example, inspection of Figure 2 indicates, on average, that the estimates of α̂3,4 are only 67.70% and 14.70% of their associated population parameters. Note that each estimate of α̂3,4 in Figure 2 was calculated based on samples of size n=250 and the formulae currently used by most commercial software packages such as SAS, SPSS, and Minitab, for computing skew and kurtosis.L-moment-based estimators, such as L-skew (τ3) and L-kurtosis (τ4), have been introduced to address the limitations associated with conventional moment-based estimators [12, 13]. Specifically, some of the advantages that L-moments have over conventional moments are that they (a) exist whenever the mean of the distribution exists, (b) are nearly unbiased for all sample sizes and distributions, and (c) are more robust in the presence of outliers. For example, the estimates of τ̂3,4 in Figure 2 are relatively much closer to their respective parameters τ3,4 with much smaller standard errors than their corresponding conventional moment-based analogs α̂3,4. More specifically, the estimates of τ̂3,4 that were simulated are, on average, 98.89% and 99.23% of their parameters.In the context of multivariate data generation, the methodology has been developed for simulating symmetric (or asymmetric) GLDs based on (1.1) with specified Pearson correlation structures [2, 14]. This methodology is based on conventional product moments and the popular NORTA (NORmal To Anything, [15]) approach, which begins with generating multivariate standard normal deviates. However, the NORTA approach is not without its limitations. Specifically, one limitation arises because the Pearson correlation is not invariant under nonlinear strictly increasing transformations such as (1.1). As such, the NORTA approach must begin with the computation of an intermediate  correlation (IC) matrix, which is different than the specified correlation matrix between the GLDs. The purpose of the IC matrix is to adjust for the effect of the transformation in (1.1) such that the resulting GLDs have their specified skew, kurtosis, and Pearson correlation matrix.Two additional limitations associated with the NORTA approach in this context are that solutions to an IC matrix may neither (a) exist in the range of [–1, +1] as the absolute values of the ICs must be greater than (or equal to) their specified Pearson correlations nor (b) yield a positive definite IC matrix albeit the specified Pearson correlation matrix is positive definite [16]. Further, these two problems can be exacerbated when heavy tailed distributions are involved in the computation of ICs as functions performing numerical integration can more frequently either fail to converge or yield incorrect solutions. In contradistinction, it has been demonstrated in the context of the L-correlation that the limitations associated with the NORTA approach are less pronounced because the solution values of an IC matrix are in closer proximity to their specified (positive definite) L-correlation matrix [17].In view of the above, the present aim is to derive the double GLD family of distributions based on (1.8) in the contexts of L-moment and L-correlation theory. Specifically, the purpose of this paper is to develop the methodology and a procedure for simulating double GLDs with specified L-moments and L-correlations. The primary advantages of the proposed procedure are that estimates of L-skew, L-kurtosis, and L-correlation are less biased and more efficient.The remainder of this paper is organized into four sections. The next section provides a summary of univariateL-moment theory and the derivations of the system of equations and boundary conditions for generating double GLDs with specified values of L-skew and L-kurtosis. The section thereafter introduces the coefficient of L-correlation and the equation is subsequently derived for determining ICs for specified L-correlations between double GLDs. The steps for implementing the proposed L-moment procedure are subsequently described. A numerical example and results of a simulation are also provided to confirm the derivations and compare the new L-moment-based procedure with the conventional product-moment-based procedure. In the last section, the results of the simulation are discussed.


## 2.  Methodology

### 2.1.  Preliminaries

LetX1,…,Xj,…,Xn be iid random variables each with continuous pdf f(x), cdf F(x), order statistics denoted as X1:n≤⋯≤Xj:n≤⋯≤Xn:n, and L-moments defined in terms of either linear combinations of (a) expectations of order statistics or (b) probability weighted moments (βi). For the purposes considered herein, the first four L-moments associated with Xj:n are expressed as [13, pages 20–22](2.1)Λ1=E[X1:1]=β0,(2.2)Λ2=12E[X2:2-X1:2]=2β1-β0,(2.3)Λ3=13E[X3:3-2X2:3+X1:3]=6β2-6β1+β0,(2.4)Λ4=14E[X4:4-3X3:4+3X2:4-X1:4]=20β3-30β2+12β1-β0,
where the βi are determined from(2.5)βi=∫x{F(x)}if(x)dx,
where i=0,…,3. The coefficients associated with βi in (2.5) are obtained from shifted orthogonal Legendre polynomials and are computed as shown in [13, page 20].TheL-moments Λ1 and Λ2 in (2.1) and (2.2) are measures of location and scale and are the arithmetic mean and one-half the coefficient of mean difference, respectively. Higher-order L-moments are transformed to dimensionless quantities referred to as L-moment ratios defined as τr=Λr/Λ2 for r≥3, and where τ3 and τ4 are the analogs to the conventional measures of skew and kurtosis. In general, L-moment ratios are bounded in the interval -1<τr<1 as is the index of L-skew (τ3) where a symmetric distribution implies that all L-moment ratios with odd subscripts are zero. Other smaller boundaries can be found for more specific cases. For example, the index of L-kurtosis (τ4) has the boundary condition for continuous distributions of [18](2.6)5τ32-14<τ4<1.


### 2.2.L-Moments for Double GLDs

The family of double GLDs, associated with (1.8), based on L-moments is less restrictive than the family based on conventional moments insofar as we may consider the nonzero parameters of λℒ,λℛ>-1 for any distribution with finite k order L-moments rather than λℒ,λℛ>-1/k for the kth-order conventional moment to exist. This advantage is attributed to Hosking’s Theorem 1 [12] which states that if the mean (Λ1) exists, then all other L-moments will have finite expectations.As such, the family of double GLDs can be derived in the context ofL-moments by defining the probability weighted moments based on (2.5) in terms of q(u) in (1.1) and the regular uniform pdf and cdf as(2.7)βi=∫01/2q(u,λL,δ)uidu+∫1/21q(u,λR,δ)uidu.
Integrating (2.7) for i=0,1,2,3 and using (2.1)–(2.4) yields (2.8)Λ1=1-2λLδ(1+λL)4λL+2λR-1δ(1+λR)4λR,(2.9)Λ2=2λL-2δ(2+3λL+λL2)+2λR-2δ(2+3λR+λR2),(2.10)τ3={((11)2λL+1-49+(3)2λR+2+(2λL+1-5)λR(6+λR))12(2λL-2λR)+6λL2(5+2λR+1(λR-1)+λR)+λL3(5+2λR+1(λR-1)+λR)+λL(49-(11)2λR+1+(11)2λR+1λR-6λR2-λR3-2λL+1(1+λR)(2+λR)(3+λR))+λR((11)2λL+1-49+(3)2λR+2+(2λL+1-5)λR(6+λR))}/{2(3+λL)(3+λR)(2λR(1+λL)(2+λL)+2λL(1+λR)(2+λR))},(2.11)τ4={(2λL(λL-2)(λL-1)(1+λL)(2+λL)(3+λL)(4+λL)+2λR(λR-2)(λR-1)(1+λR)(2+λR)(3+λR)(4+λR))(1+λL)(2+λL)(1+λR)(2+λR)×(2λL(λL-2)(λL-1)(1+λL)(2+λL)(3+λL)(4+λL)+2λR(λR-2)(λR-1)(1+λR)(2+λR)(3+λR)(4+λR))}/{2λR(1+λL)(2+λL)+2λL(1+λR)(2+λR)}.
Thus, given specified values of τ3 and τ4, (2.10) and (2.11) can be numerically solved for the corresponding values of λℒ and λℛ. Note that the values of L-skew (τ3) and L-kurtosis (τ4) in (2.10) and (2.11) are independent of the value of δ selected in (1.8). Further, inspection of (2.10) and (2.11) indicates that interchanging values for the parameters λℒ and λℛ reverses the direction of τ3 and has no effect on τ4.For the special case ofλℒ=λℛ the double GLD is symmetric where τ3=0 in (2.10) and τ4 in (2.11) will simplify to the expression(2.12)τ4=(λL-2)(λL-1)(λL+3)(λL+4).
Differentiating (2.12) with respect to λℒ and equating the resulting expression to zero and solving for λℒ yields a minimum value of L-kurtosis of(2.13)min⁡⁡(τ4)=12-5612+56=-0.010205…,
where λℒ=λℛ=6-1. As such, using (2.10) and (2.11) with λℒ=6-1 and λℛ∈(-1,6-1), provided in Figure 3 is a graph of the region for feasible combinations of τ3 and τ4 for double GLDs. Feasible combinations of τ3 and τ4 for (2.10) and (2.11) will lie in the region above the curve graphed in the |τ3|,τ4 plane of Figure 3.Figure 3
Graph of the region for feasible combinations of (absolute value)L-skew |τ3| and L-kurtosis τ4. A double GLD will lie in the area above the curve graphed in the |τ3| and τ4 plane.Provided in Figure4 are some examples of various double GLDs. These distributions are used in the simulation portion of this study in Section 4. The next section begins with an introduction to the L-correlation.Four asymmetric double GLDs with their conventional andL-moment parameters of skew (α3) and L-skew (τ3), kurtosis (α4), and L-kurtosis (τ4), and corresponding shape parameters for (1.8). Note that the height of each distribution is δ=1/2π.
(a)(b)(c)(d)


## 2.1.  Preliminaries

LetX1,…,Xj,…,Xn be iid random variables each with continuous pdf f(x), cdf F(x), order statistics denoted as X1:n≤⋯≤Xj:n≤⋯≤Xn:n, and L-moments defined in terms of either linear combinations of (a) expectations of order statistics or (b) probability weighted moments (βi). For the purposes considered herein, the first four L-moments associated with Xj:n are expressed as [13, pages 20–22](2.1)Λ1=E[X1:1]=β0,(2.2)Λ2=12E[X2:2-X1:2]=2β1-β0,(2.3)Λ3=13E[X3:3-2X2:3+X1:3]=6β2-6β1+β0,(2.4)Λ4=14E[X4:4-3X3:4+3X2:4-X1:4]=20β3-30β2+12β1-β0,
where the βi are determined from(2.5)βi=∫x{F(x)}if(x)dx,
where i=0,…,3. The coefficients associated with βi in (2.5) are obtained from shifted orthogonal Legendre polynomials and are computed as shown in [13, page 20].TheL-moments Λ1 and Λ2 in (2.1) and (2.2) are measures of location and scale and are the arithmetic mean and one-half the coefficient of mean difference, respectively. Higher-order L-moments are transformed to dimensionless quantities referred to as L-moment ratios defined as τr=Λr/Λ2 for r≥3, and where τ3 and τ4 are the analogs to the conventional measures of skew and kurtosis. In general, L-moment ratios are bounded in the interval -1<τr<1 as is the index of L-skew (τ3) where a symmetric distribution implies that all L-moment ratios with odd subscripts are zero. Other smaller boundaries can be found for more specific cases. For example, the index of L-kurtosis (τ4) has the boundary condition for continuous distributions of [18](2.6)5τ32-14<τ4<1.


## 2.2.L-Moments for Double GLDs

The family of double GLDs, associated with (1.8), based on L-moments is less restrictive than the family based on conventional moments insofar as we may consider the nonzero parameters of λℒ,λℛ>-1 for any distribution with finite k order L-moments rather than λℒ,λℛ>-1/k for the kth-order conventional moment to exist. This advantage is attributed to Hosking’s Theorem 1 [12] which states that if the mean (Λ1) exists, then all other L-moments will have finite expectations.As such, the family of double GLDs can be derived in the context ofL-moments by defining the probability weighted moments based on (2.5) in terms of q(u) in (1.1) and the regular uniform pdf and cdf as(2.7)βi=∫01/2q(u,λL,δ)uidu+∫1/21q(u,λR,δ)uidu.
Integrating (2.7) for i=0,1,2,3 and using (2.1)–(2.4) yields (2.8)Λ1=1-2λLδ(1+λL)4λL+2λR-1δ(1+λR)4λR,(2.9)Λ2=2λL-2δ(2+3λL+λL2)+2λR-2δ(2+3λR+λR2),(2.10)τ3={((11)2λL+1-49+(3)2λR+2+(2λL+1-5)λR(6+λR))12(2λL-2λR)+6λL2(5+2λR+1(λR-1)+λR)+λL3(5+2λR+1(λR-1)+λR)+λL(49-(11)2λR+1+(11)2λR+1λR-6λR2-λR3-2λL+1(1+λR)(2+λR)(3+λR))+λR((11)2λL+1-49+(3)2λR+2+(2λL+1-5)λR(6+λR))}/{2(3+λL)(3+λR)(2λR(1+λL)(2+λL)+2λL(1+λR)(2+λR))},(2.11)τ4={(2λL(λL-2)(λL-1)(1+λL)(2+λL)(3+λL)(4+λL)+2λR(λR-2)(λR-1)(1+λR)(2+λR)(3+λR)(4+λR))(1+λL)(2+λL)(1+λR)(2+λR)×(2λL(λL-2)(λL-1)(1+λL)(2+λL)(3+λL)(4+λL)+2λR(λR-2)(λR-1)(1+λR)(2+λR)(3+λR)(4+λR))}/{2λR(1+λL)(2+λL)+2λL(1+λR)(2+λR)}.
Thus, given specified values of τ3 and τ4, (2.10) and (2.11) can be numerically solved for the corresponding values of λℒ and λℛ. Note that the values of L-skew (τ3) and L-kurtosis (τ4) in (2.10) and (2.11) are independent of the value of δ selected in (1.8). Further, inspection of (2.10) and (2.11) indicates that interchanging values for the parameters λℒ and λℛ reverses the direction of τ3 and has no effect on τ4.For the special case ofλℒ=λℛ the double GLD is symmetric where τ3=0 in (2.10) and τ4 in (2.11) will simplify to the expression(2.12)τ4=(λL-2)(λL-1)(λL+3)(λL+4).
Differentiating (2.12) with respect to λℒ and equating the resulting expression to zero and solving for λℒ yields a minimum value of L-kurtosis of(2.13)min⁡⁡(τ4)=12-5612+56=-0.010205…,
where λℒ=λℛ=6-1. As such, using (2.10) and (2.11) with λℒ=6-1 and λℛ∈(-1,6-1), provided in Figure 3 is a graph of the region for feasible combinations of τ3 and τ4 for double GLDs. Feasible combinations of τ3 and τ4 for (2.10) and (2.11) will lie in the region above the curve graphed in the |τ3|,τ4 plane of Figure 3.Figure 3
Graph of the region for feasible combinations of (absolute value)L-skew |τ3| and L-kurtosis τ4. A double GLD will lie in the area above the curve graphed in the |τ3| and τ4 plane.Provided in Figure4 are some examples of various double GLDs. These distributions are used in the simulation portion of this study in Section 4. The next section begins with an introduction to the L-correlation.Four asymmetric double GLDs with their conventional andL-moment parameters of skew (α3) and L-skew (τ3), kurtosis (α4), and L-kurtosis (τ4), and corresponding shape parameters for (1.8). Note that the height of each distribution is δ=1/2π.
(a)(b)(c)(d)


## 3.  TheL-Correlation for Double GLDs

The coefficient ofL-correlation (see [19]) is introduced by considering two random variables Yj and Yk with continuous distribution functions F(Yj) and F(Yk), respectively. The second L-moments of Yj and Yk can alternatively be expressed as(3.1)Λ2(Yj)=2Cov⁡(Yj,F(Yj)),(3.2)Λ2(Yk)=2Cov⁡(Yk,F(Yk)).
The second L-comoments of Yj toward Yk and Yk toward Yj are(3.3)Λ2(Yj,Yk)=2Cov⁡(Yj,F(Yk)),(3.4)Λ2(Yk,Yj)=2Cov⁡(Yk,F(Yj)).
As such, the L-correlations of Yj toward Yk and Yk toward Yj are expressed as(3.5)ηjk=Λ2(Yj,Yk)Λ2(Yj),(3.6)ηkj=Λ2(Yk,Yj)Λ2(Yk).
The L-correlation in (3.5) or (3.6) is bounded such that -1≤ηjk≤1 where a value of ηjk=1 (ηjk=-1) indicates a strictly increasing (decreasing) monotone relationship between the two variables. In general, we would also note that ηjk≠ηkj.In the context of theL-moment-based double GLD, suppose it is desired to simulate T distributions of the form in (1.8) with a specified L-correlation matrix and where each distribution has its own specified values of τ3 and τ4. Let Z1,…,ZT denote standard normal variables where the distribution functions and bivariate density function associated with Zj and Zk are expressed as
(3.7)Φ(zj)=Pr⁡{Zj≤zj}=∫-∞zj(2π)-1/2exp⁡⁡{-wj22}dwj,(3.8)Φ(zk)=Pr⁡{Zk≤zk}=∫-∞zk(2π)-1/2exp⁡{-wk22}⁡dwk,(3.9)fjk=(2π(1-ρjk2)1/2)-1exp⁡{-(2(1-ρjk2))-1(zj2+zk2-2ρjkzjzk)}.
Using (3.7), it follows that the jth double GLD associated with (1.8) can be expressed as qj(Φ(zj)) because Φ(zj)~U(0,1). As such, using (3.5), the L-correlation of qj(Φ(zj)) toward qk(Φ(zk)) can be evaluated using solved values of λℒj and λℛj for qj(Φ(zj)), a specified intermediate correlation (IC) ρjk in (3.9), and the following integral expressed as(3.10)ηjk=2π∬-∞+∞xj(qj(Φ(zj),λLj,λRj))Φ(zk)fjkdzjdzk.
The double GLD in (3.10) is standardized by a linear transformation such that it has a mean of zero and one-half the coefficient of mean difference equal to that of the unit-normal distribution as(3.11)xj(qj(Φ(zj),λLj,λRj))=ξ(qj(Φ(zj),λLj,λRj)-Λ1),
where Λ1 is the mean from (2.8) and ξ is a constant that scales Λ2 in (2.9) and in the denominator of (3.5) to 1/π as(3.12)ξ=δ4(1+λLj)(2+λLj)(1+λRj)(1+λRj)2λRj(1+λLj)(2+λLj)π+2λLj(1+λRj)(1+λRj)π.
Analogously, the L-correlation of qk(Φ(zk)) toward qj(Φ(zj)) is expressed as(3.13)ηkj=2π∬-∞+∞xk(qk(Φ(zk),λLk,λRk))Φ(zj)fjkdzkdzj.
Note for the special case that if qj(Φ(zj)) in (3.10) and qk(Φ(zk)) in (3.13) have the same parameters that is, λℒj=λℒk and λℛj=λℛk, then ηjk=ηkj. Provided in Algorithm 1 is source code written in Mathematica [20] that implements the computation of an IC (ρjk) based on (3.10). The details for simulating double GLDs with specifiedL-correlations are described in the next section.Algorithm 1:Mathematica source code for computing intermediate correlations for specified L-correlations. The example is for distribution j=3 towards distribution k=4 (η34) in Figure 4. See also Table 3(b).
Φj=CDF[NormalDistribution[0,1],Zj];Φk=CDF[NormalDistribution[0,1], Zk];λℒ=0.2779733598283232;λℛ=-0.08329479143312762;qℒ=(Φjλℒ-(1-Φj)λℒ)/(λℒ×2(3/2)-λℒ/π);qℛ=(Φjλℛ-(1-Φj)λℛ)/(λℛ×2(3/2)-λℛ/π);(* Standardizing constantsΛ1 from Eq. (2.8) and  ξ from Eq. (3.12) *)Xℒ=ξ(qℒ-Λ1);Xℛ=ξ(qℛ-Λ1);(* Intermediate Correlation *)ρjk=0.395685;Needs[“MultivariateStatistics`”]fjk=PDF[MultinormalDistribution[{0,0},{{1,ρjk},{ρjk,1}}],{Zj, Zk}];(* Compute the specifiedL-correlation *)ηjk=2π*NIntegrate[Piecewise[{{Xℒ, Φj≤0.5}, {Xℛ, Φj>0.5}}]×Φk×fjk,{Zj,-10, 10}, {Zk, −10, 10}, Method→MultiDimensional]0.40


## 4.  The Procedure and Simulation Study

To implement the procedure for simulating double GLDs with specifiedL-moments and L-correlations we suggest the following six steps. (1)
Specify theL-moments for T transformations of the form in (1.8), that is, q1(Φ(z1)),…,qT(Φ(zT))and obtain the solutions for the parameters of λℒj and λℛjby solving (2.10) and (2.11) using the specified values of L-skew (τ3) and L-kurtosis (τ4) for each distribution. Specify a T×T matrix of L-correlations (ηjk) for qj(Φ(zj)) toward qk(Φ(zk)), where j<k∈{1,2,…,T}.(2)
Compute the (Pearson) intermediate correlations (ICs)ρjk by substituting the solutions of λℒj and λℛj from Step (1) into (3.10) and then numerically integrate to solve for ρjk (see Algorithm 1 for an example). Repeat this step separately for all T(T-1)/2 pairwise combinations of correlations.(3)
Assemble the ICs into aT×T matrix and decompose this matrix using a Cholesky factorization. Note that this step requires the IC matrix to be positive definite.(4)
Use the results of the Cholesky factorization from Step (3) to generateT standard normal variables (Z1,…,ZT) correlated at the intermediate levels as follows:
(4.1)Z1=a11V1Z2=a12V1+a22V2⋮Zj=a1jV1+a2jV2+⋯+aijVi+⋯+ajjVj⋮ZT=a1TV1+a2TV2+⋯+aiTVi+⋯+ajTVj+⋯+aTTVT,
where V1,…,VT are independent standard normal random variables and where aij represents the element in the ith row and the jth column of the matrix associated with the Cholesky factorization performed in Step 3.(5)
SubstituteZ1,…,ZT from Step (4) into the following Taylor series-based expansion for the standard normal cdf [21]:
(4.2)Φ(Zj)=(12)+ϕ(Zj){Zj+Zj33+Zj5(3⋅5)+Zj7(3⋅5⋅7)+⋯},
where ϕ(Zj) denotes the standard normal pdf and where the absolute error associated with (4.2) is less than 8×10-16.(6)
Substitute the regular uniform deviates,Φ(Zj), generated from Step (5) into the T equations of the form in (1.8), as noted in Step (1), to generate the double GLDs with the specified L-moments and L-correlations.To demonstrate the steps above and evaluate the proposed procedure, a comparison between the newL-moment and conventional product moment-based procedures is subsequently described. Specifically, the distributions in Figure 4 are used as a basis for a comparison using the specified correlation matrices in Table 1 where both strong and moderate levels of correlation are considered. Tables 2 and 3 give the solved IC matrices for the conventional moment- and L-moment-based procedures, respectively. See Algorithm 2 for the algorithm and an example for computing ICs for the conventional procedure. Tables 4 and 5 give the results of the Cholesky decompositions on the IC matrices, which are then used to create Z1,…,Z4 with the specified ICs by making use of the formulae given in (4.1) of Step (4) with T=4. The values of Z1,…,Z4 are subsequently transformed to Φ(Z1),…,Φ(Z4) using (4.2) and then substituted into equations of the form in (1.8) to produce q1(Φ(Z1)),…,q4(Φ(Z4)) for both procedures.Table 1
Specified correlation matrices for the distributions in Figure4.
(a)12341120.70130.700.70140.850.700.701(b)12341120.40130.500.40140.600.500.401Table 2
Intermediate correlations for the conventional moment procedure.
(a)12341120.815469130.7285370.762929140.8990990.7419350.7095971(b)12341120.441594130.5276600.429990140.6363150.5270470.4067861Table 3
Intermediate correlations for theL-moment procedure.
(a)12341120.686456130.6864560.694980140.8410130.6949800.6954271(b)12341120.387335130.4858400.395135140.5854680.4946090.3956851Table 4
Cholesky decomposition for the conventional moment procedure.
(a)a11=1a12=0.815469a13=0.728537a14=0.8990990a22=0.578801a23=0.291689a24=0.01513300a33=0.619800a34=0.080932000a44=0.429933(b)a11=1a12=0.441594a13=0.527660a14=0.6363150a22=0.897215a23=0.219544a24=0.27424200a33=0.820594a34=0.013185000a44=0.720917Table 5
Cholesky decomposition for theL-moment procedure.
(a)a11=1a12=0.686456a13=0.686456a14=0.8410130a22=0.727171a23=0.307710a24=0.16180700a33=0.658857a34=0.103693000a44=0.505730(b)a11=1a12=0.387335a13=0.485840a14=0.5854680a22=0.921939a23=0.224475a24=0.29051500a33=0.844731a34=0.054488000a44=0.754890Algorithm 2:Mathematica source code for computing intermediate correlations for specified conventional Pearson correlations. The example is for distributions j=3 and k=4 (ρ34*) in Figure 4. See also Table 2(b).
Φj=CDF[NormalDistribution[0,1], Zj];Φk=CDF[NormalDistribution[0,1], Zk];λℒj=0.2779733598283232;λℛj=-0.08329479143312762;λℒk=0.3662599069298994;λℛk=0.1297771263080004;qℒj=(Φjλℒj-(1-Φj)λℒj)/(λℒj×2(3/2)-λℒj/π);qℛj=(Φjλℛj-(1-Φj)λℛj)/(λℛj×2(3/2)-λℛj/π);qℒk=(Φkλℒk-(1-Φk)λℒk)/(λℒk×2(3/2)-λℒk/π);qℛk=(Φkλℛk-(1-Φk)λℛk)/(λℛk×2(3/2)-λℛk/π);(* Standardizing constantsα1 and α2 from Equation (A.2) in the Appendix *)Xℒj=(qℒj-α1j)/α2j;Xℛj=(qℛj-α1j)/α2j;Xℒk=(qℒk-α1k)/α2k;Xℛk=(qℛk-α1k)/α2k;(* Intermediate Correlation *)ρjk=0.406786;Needs[“MultivariateStatistics`”]fjk=PDF[MultinormalDistribution[{0,0},{{1,ρjk},{ρjk,1}}],{Zj, Zk}];(* Compute the specified conventional Pearson correlation *)ρjk*=NIntegrate[Piecewise[{{Xℒj, Φj≤0.5}, {Xℛj, Φj>0.5}}]×Piecewise[{{Xℒk, Φk≤0.5}, {Xℛk, Φk>0.5}}]×fjk, {Zj, −10, 10}, {Zk, −10, 10}, Method→MultiDimensional]0.40In terms of the simulation, a Fortran algorithm was written for both procedures to generate 25,000 independent sample estimates for the specified parameters of (a) conventional skew (α3), kurtosis (α4), and Pearson correlation (ρjk*); (b) L-skew (τ3), L-kurtosis (τ4), and L-correlation (ηjk). The estimates of τ3,4 were based on samples of size n=250 and the estimates of ηjk were based on sample sizes of n=25 and n=1000. The estimates for α3,4 were based on Fisher’s k-statistics, that is, the formulae currently used by most commercial software packages such as SAS, SPSS, and Minitab, for computing indices of skew and kurtosis (where α3,4=0 for the standard normal distribution). The formulae used for computing estimates for τ3,4 were Headrick’s equations (2.4) and (2.6) [22]. The estimate for ρjk* was based on the usual formula for the Pearson product-moment of correlation statistic and the estimate for ηjk was computed based on (3.5) using the empirical forms of the cdfs in (3.1) and (3.3). The estimates for ρjk* and ηjk were both transformed using Fisher’s z′ transformation. Bias-corrected accelerated bootstrapped average estimates, confidence intervals (C.I.s), and standard errors were subsequently obtained for the estimates associated with the parameters (α3,4, τ3,4, zρjk*′, zηjk′) using 10,000 resamples via the commercial software package Spotfire S+ [23]. The bootstrap results for the estimates of zρjk*′ and zηjk′ were transformed back to their original metrics. Further, if a parameter (P) was outside its associated bootstrap C.I., then an index of relative bias (RB) was computed for the estimate (E) as RB=((E-P)/P)×100. The results of the simulation are reported in Tables 6, 7, 8, 9, 10, 11 and are discussed in the next section.Table 6
Skew (α3) and Kurtosis (α4) results for the conventional moment procedure. Results are based on a sample size of n=250.
DistParameterEstimate95% Bootstrap C.I.Standard errorRelative bias %1α3=2.701.731.72, 1.750.0078−35.92α4=87.48.938.75, 9.120.0939−89.782α3=-1.38−1.23−1.24, −1.220.0039−10.87α4=5.673.783.71, 3.850.0360−33.333α3=0.8830.8240.819, 0.8290.0026−6.68α4=2.552.042.01, 2.080.0185−20.004α3=0.2950.2920.291, 0.2940.0008−1.02α4=-0.232−0.227−0.231, −0.2240.0018−2.16Table 7
L-skew (τ3) and L-kurtosis (τ4) results. Results are based on a sample size of n=250.
DistParameterEstimate95% Bootstrap C.I.Standard errorRelative bias %1τ3=0.150.14790.1473, 0.14840.0003−1.40τ4=0.200.19870.1982, 0.19910.0002−0.652τ3=-0.15−0.1487−0.1492, −0.14830.0002−0.87τ4=0.150.14940.1490, 0.14970.0002−0.403τ3=0.100.09920.0988, 0.09960.0002−0.80τ4=0.150.14970.1494, 0.15010.0002—4τ3=0.050.04980.0495, 0.05010.0002—τ4=0.100.10020.0999, 0.10040.0001—Table 8
Correlation (strong) results for the conventional moment procedure.
(a)
n=25
ParameterEstimate95% Bootstrap C.I.Standard errorRelative bias %ρ12*=0.700.76060.7597, 0.76150.001118.66ρ13*=0.700.72770.7263, 0.72910.001463.96ρ14*=0.850.88540.8847, 0.88580.001304.16ρ23*=0.700.72900.7279, 0.73000.001134.14ρ24*=0.700.72090.7198, 0.72200.001162.99ρ34*=0.700.71200.7108, 0.71340.001361.71(b)
n=1000
ParameterEstimate95% Bootstrap C.I.Standard errorRelative bias %ρ12*=0.700.70670.7063, 0.70710.000380.96ρ13*=0.700.70360.7033, 0.70390.000300.51ρ14*=0.850.85620.8560, 0.85650.000540.73ρ23*=0.700.70140.7012, 0.70160.000200.20ρ24*=0.700.70100.7008, 0.70120.000200.14ρ34*=0.700.70050.7002, 0.70060.000200.07Table 9
Correlation (strong) results for theL-moment procedure.
(a)
n=25
ParameterEstimate95% Bootstrap C.I.Standard errorRelative bias %η12=0.700.70820.7068, 0.70970.001481.17η13=0.700.70780.7064, 0.70930.001491.11η14=0.850.85540.8546, 0.85620.001510.64η23=0.700.70910.7077, 0.71050.001461.30η24=0.700.70900.7076, 0.71040.001471.29η34=0.700.70970.7083, 0.71110.001481.39(b)
n=1000
ParameterEstimate95% Bootstrap C.I.Standard errorRelative bias %η12=0.700.70000.6998, 0.70020.00023—η13=0.700.70020.7000, 0.70040.00023—η14=0.850.85020.8500, 0.85030.00022—η23=0.700.69990.6997, 0.70020.00022—η24=0.700.70050.6998, 0.70030.00022—η34=0.700.70010.6999, 0.70030.00021—Table 10
Correlation (moderate) results for the conventional moment procedure.
(a)
n=25
ParameterEstimate95% Bootstrap C.I.Standard errorRelative bias %ρ12*=0.400.42070.4187, 0.42270.001235.18ρ13*=0.500.52650.5244, 0.52850.001425.30ρ14*=0.600.63080.6292, 0.63240.001355.13ρ23*=0.400.41780.4156, 0.41970.001264.45ρ24*=0.500.51650.5146, 0.51830.001273.30ρ34*=0.400.41030.4081, 0.41250.001342.58(b)
n=1000
ParameterEstimate95% Bootstrap C.I.Standard errorRelative bias %ρ12*=0.400.39060.3903, 0.39100.00020−2.35ρ13*=0.500.50290.5026, 0.50330.000250.58ρ14*=0.600.60350.6032, 0.60390.000260.58ρ23*=0.400.40050.4001, 0.40070.000180.13ρ24*=0.500.50060.5003, 0.50090.000190.12ρ34*=0.400.40030.3999, 0.40060.000200.08Table 11
Correlation (moderate) results for theL-moment procedure.
(a)
n=25
ParameterEstimate95% Bootstrap C.I.Standard errorRelative bias %η12=0.400.40840.4061, 0.41090.001482.10η13=0.500.50730.5051, 0.50950.001491.46η14=0.600.60800.6061, 0.60980.001501.33η23=0.400.40830.4060, 0.41070.001452.08η24=0.500.50940.5073, 0.51160.001471.88η34=0.400.40890.4066, 0.41120.001422.23(b)
n=1000
ParameterEstimate95% Bootstrap C.I.Standard errorRelative bias %η12=0.400.40000.3996, 0.40040.00023—η13=0.500.50020.4999, 0.50050.00023—η14=0.600.60020.5999, 0.60050.00022—η23=0.400.39970.3994, 0.40010.00022—η24=0.500.50010.4998, 0.50040.00022—η34=0.400.40000.3996, 0.40030.00021—


## 5.  Discussion and Conclusion

One of the advantages thatL-moment ratios have over conventional moment-based estimators is that they can be far less biased when sampling is from distributions with heavy tails [13, 19]. Inspection of the simulation results in Tables 6 and 7 clearly indicates that this is the case. Specifically, the superiority that estimates of L-moment ratios (τ3,τ4) have over their corresponding conventional moment-based counterparts (α3,α4) is obvious. For example, with samples of size n=25 the estimates of skew and kurtosis for Distribution 1 were, on average, only 64.07% and 10.22% of their associated population parameters, whereas the estimates of L-skew and L-kurtosis were 98.6% and 99.35% of their respective parameters. It is also evident from comparing Tables 6 and 7 that L-skew and L-kurtosis are more efficient estimators as their standard errors are substantially smaller and more stable than the conventional moment-based estimators of skew and kurtosis.Presented in Tables8, 9, 10, 11 are the results associated with the conventional Pearson and L-correlations. Overall inspection of these tables indicates that the L-correlation is superior to the Pearson correlation in terms of relative bias. For example, for strong correlations (n=25) the relative bias for the two heavy tailed distributions (i.e., distributions 1 and 2) was 8.66% for the Pearson correlation compared to only 1.17% for the L-correlation. Further, for large sample sizes (n=1000), the L-correlation bootstrap C.I.s contained all of the population parameters, whereas the Pearson correlation C.I.s contained none of the parameters. It is also noted that the variability of the L-correlation appears to be more stable than that of the Pearson correlation both within and across the different conditions.In summary, the newL-moment-based procedure is an attractive alternative to the traditional conventional moment-based procedure. In particular, the L-moment-based double GLD procedure has distinct advantages when distributions with heavy tails are used. Finally, we note that Mathematica Version 8.0.1 [20] source code is available from the authors for implementing both the conventional and new L-moment-based procedures.



---

*Source: 725754-2012-05-07.xml*