# A Method for Simulating Nonnormal Distributions with SpecifiedL-Skew, L-Kurtosis, and L-Correlation

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

---

## Abstract

This paper introduces two families of distributions referred to as the symmetricκ and asymmetric κℒ-κℛ distributions. The families are based on transformations of standard logistic pseudo-random deviates. The primary focus of the theoretical development is in the contexts of L-moments and the L-correlation. Also included is the development of a method for specifying distributions with controlled degrees of L-skew, L-kurtosis, and L-correlation. The method can be applied in a variety of settings such as Monte Carlo studies, simulation, or modeling events. It is also demonstrated that estimates of L-skew, L-kurtosis, and L-correlation are superior to conventional product-moment estimates of skew, kurtosis, and Pearson correlation in terms of both relative bias and efficiency when moderate-to-heavy-tailed distributions are of concern.


---

## Body

## 1. Introduction

Monte Carlo investigations often require the need for univariate or multivariate nonnormal distributions with specified conventional product moments. For example, it is common practice for methodologists to investigate the Type I error and power properties associated with inferential statistics under various degrees of nonnormality. In many cases, these investigations may only require an elementary transformation on standard normal or zero-one uniform deviates to produce nonnormal distributions with specified values of conventional skew and kurtosis. More specifically, the popular power method (PM) [1–3] and the generalized lambda distribution (GLD) [4–6] are two basic transformations that could be used for this purpose.However, conventional-moment-based estimators of skew and kurtosis have unfavorable attributes insofar as they can be substantially biased, have high variance, or can be influenced by outliers [7, 8]. In view of this, L-moment-based estimators such as L-skew and L-kurtosis were introduced to address the limitations associated with conventional moment-based estimators. 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 [7–10]. In view of these advantages, both the PM and GLD transformations have been characterized in the context of L-moments to enable the simulation of univariate nonnormal distributions with specified values of L-skew and L-kurtosis [6, 11–13].In addition to issues of a transformation’s simplicity and ease of execution, another important criterion is its ability to simulate nonnormal distributions with specified correlation structures. To meet this criterion, methodologists have extended the conventional-moment-based PM and GLD transformations from univariate to multivariate data generation in the context of the Pearson correlation [2, 3, 14–17]. However, and although these product-moment-based methods are commonly used [3, pages 2-3], they would not make feasible extensions of the univariate L-moment-based PM or GLD transformations because the second L-moment and the L-correlation are based on the coefficient of mean difference and not the variance or standard deviation [7, 8, 18].Thus, the primary goal of this paper is to derive two families of distributions based on transformations that produce continuous symmetric and asymmetric nonnormal distributions with specified values ofL-skew, L-kurtosis, and L-correlation. These two families of distributions are referred to herein as the symmetric κ and asymmetric κℒ-κℛ distributions. The proposed transformations are computationally efficient because they only require the knowledge of one or two parameters (κ, or κℒ, κℛ) and an algorithm that generates standard logistic pseudo-random deviates. We would note that the proposed transformations are analogs to the L-moment-based Tukey h and hh distributions [19]. However, the κ and κℒ-κℛ transformations can generate more leptokurtic and elongated distributions as they are based on the logistic distribution as opposed to the standard normal based h or hh distributions.The remainder of the paper is outlined as follows. In Section2, the parametric forms of the probability density function, cumulative distribution function, and other properties associated with the κ and κℒ-κℛ distributions are derived. In Section 3, a summary of univariate L-moment theory is first provided. The derivations of the systems of equations for specifying values of L-skew and L-kurtosis for the κ and κℒ-κℛ distributions are subsequently provided. The boundary region for feasible combinations of L-skew and L-kurtosis is also provided. In Section 4, the coefficient of L-correlation is first introduced and then the equations are developed for determining intermediate correlations for specified L-correlations between κ and (or) κℒ-κℛ distributions. In Section 5, the steps for implementing the new method are described. A numerical example and the results of a simulation are provided to confirm the derivations and compare the new method with its conventional-moment-based counterpart. In Section 6, the results of the simulation are discussed.


## 2. Methodology

The derivation of the probability density function (pdf) and the cumulative distribution function (cdf) associated with theκ and κℒ-κℛ families of distributions begins with the following definitions.Definition 2.1.
LetX be a random variable that has a standard logistic distribution with pdf and cdf expressed as
(2.1)fX(x)=e-x(1+e-x)2(2.2)FX(x)=Pr(X≤x)=1(1+e-x),-∞<x<+∞.
Let x=(u,v) be the auxiliary variable that maps the parametric curves of (2.1) and (2.2) as
(2.3)f:x↦ℜ2≔fX(x)=fX(u,v)=fX(x,fX(x)),F:x↦ℜ2≔FX(x)=FX(u,v)=FX(x,FX(x)).Definition 2.2.
The analytical and empirical forms of the quantile functions for symmetricκ distributions are defined as
(2.4)q(x)=qκ(x)=xeκ|x|,(2.5)q(X)=qκ(X)=Xeκ|X|,
where κ is a real-valued parameter that controls the tail weight of a distribution.Definition 2.3.
The analytical and empirical forms of the quantile functions for asymmetricκℒ-κℛ distributions are defined as
(2.6)q(x)=qκℒ,κℛ(x)={xeκℒ|x|,forx≤0,xeκℛ|x|,forx≥0,(2.7)q(X)=qκℒ,κℛ(X)={Xeκℒ|X|,forX≤0,Xeκℛ|X|,forX≥0,
where κℒ (κℛ) is the real-valued parameter that controls the left (right) tail of a κℒ-κℛ distribution. Note that if κℒ=κℛ, then (2.6) and (2.7) are equivalent to the symmetric κ distribution in (2.4) and (2.5).The explicit forms of the derivatives associated with (2.4) and (2.6) are
(2.8)q′(x)=qκ′(x)=eκ|x|+xκeκ|x|d|x|dx,q′(x)=qκℒ,κℛ′(x)={eκℒ|x|-xκℒeκℒ|x|,forx<0,eκℛ|x|+xκℛeκℛ|x|,forx>0,
where it is assumed that q′(x)>0, that is, the transformations in (2.4) and (2.6) are strictly increasing monotone functions and where q′(x)=1 at x=0 in (2.8).Proposition 2.4.
If the compositionsf∘q and F∘q map the parametric curves of fq(X)(q(x)) and Fq(X)(q(x)) where q(x)=q(u,v) as
(2.9)f∘q:q(x)↦ℜ2≔fq(X)(q(x))=fq(X)(q(u,v))=fq(X)(q(x),fX(x)q′(x)),(2.10)F∘q:q(x)↦ℜ2≔Fq(X)(q(x))=Fq(X)(q(u,v))=Fq(X)(q(x),FX(x)),
then (2.9) and (2.10) are the parametric forms of the pdf and cdf associated with the quantile functions in q(X) in (2.5) and (2.7), respectively.Proof.
It is first shown thatfq(X)(q(x),fX(x)/q′(x)) in (2.9) has the properties:Property 1.
∫
-
∞
+
∞
f
q
(
X
)
(
q
(
x
)
,
f
X
(
x
)
/
q
′
(
x
)
)
d
x
=
1.Property 2.
f
q
(
X
)
(
q
(
x
)
,
f
X
(
x
)
/
q
′
(
x
)
)
≥
0
,
for
-
∞
<
x
<
+
∞.
To prove Property1, let v=f(u) be a function where ∫-∞+∞f(u)du=∫-∞+∞vdu. Thus, given that u=q(x) and v=fX(x)/q′(x) in fq(X)(q(u,v)) in (2.9) we have
(2.11)∫-∞+∞fq(X)(q(x),fX(x)q′(x))dx=∫-∞+∞vdu=∫-∞+∞(fX(x)q′(x))dq(x)=∫-∞+∞(fX(x)q′(x))q′(x)d(x)=∫-∞+∞fX(x)d(x)=1,
which integrates to one because fX(x) is the standard logistic pdf in (2.1). To show Property 2, it is given by definition that fX(x)≥0 and it is assumed that  q′(x)>0. Hence, fq(X)(q(x),fX(x)/q′(x))≥0 because fX(x)/q′(x) will be nonnegative on the support of x for all x∈(-∞,+∞) and where the limx→±∞fq(X)(q(x),fX(x)/q′(x))=0 provided that κ≥0 for symmetric distributions or for the case of asymmetric distributions we must have κℒ≥0 and κℛ≥0.A corollary to Proposition2.4 is stated as follows.Corollary 2.5.
The derivative of the cdfFq(X)(q(x),FX(x)) in (2.10) is the pdf fq(X)(q(x),fX(x)/q′(x)) in (2.9).Proof.
It follows fromu=q(x) and v=FX(x) in Fq(X)(q(u,v)) in (2.10) that du=q′(x)dx and dv=fX(x)dx. Hence, using the parametric form of the derivative we have v=dv/du=fX(x)/q′(x) in (2.9). Whence, Fq(X)′(q(x),FX(x))=Fq(X)′(q(u,dv/du))=fq(X)(q(u,v))=fq(X)(q(x),fX(x)/q′(x)). Thus, fq(X)(q(x),fX(x)/q′(x)) in (2.9) and Fq(X)(q(x),FX(x)) in (2.10) are the parametric forms of the pdf and cdf associated with the quantile functions in (2.5) and (2.7).Remark 2.6.
Inspection of (2.1), (2.4), (2.6), and (2.9) indicates that the height of any symmetric κ or asymmetric κℒ-κℛ distribution at x=0 will be fX(x)/q′(x)=1/4.Definition 2.7.
If the monotonicity assumption (q′(x)>0) holds in (2.9) for all x∈(-∞,+∞), which requires κℒ,ℛ≥0, then (2.9) is defined as a global pdf.Remark 2.8.
The mode associated with a global pdf in (2.9) is located at fq(X)(q(x~),fX(x~)/q′(x~)), where x=x~ is the critical number that solves dv/du=d(fX(x)/q′(x))/du=0 and globally maximizes v=fX(x~)/q′(x~) at x=q(x~). It is noted that the pdf in (2.9) will have a global maximum because the standard logistic pdf in (2.1) has a global maximum and the transformation q(x) is assumed to be a strictly increasing monotone function.Remark 2.9.
The median associated with a global pdf in (2.9) is located at q(x=0)=0. This can be shown by letting u0.50=q(x) and v0.50=0.50=FX(x)=Pr{X≤x} denote the coordinates in the cdf in (2.10) that are associated with the 50th percentile. In general, we must have x=0 such that v0.50=0.50=FX(0)=Pr{X≤0} holds in (2.10) for the standard logistic distribution. As such, the limit of the quantile function q(x) locates the median at limx→0q(x)=0.Remark 2.10.
The monotonicity assumption (q′(x)>0) also holds in (2.9) for cases where κℒ,ℛ<0 and (1/κℒ,ℛ<x<0 or 0≤x<-1/κℒ,ℛ). This leads to Definition 2.11.Definition 2.11.
A symmetricκ or asymmetric κℒ-κℛ distribution has a local pdf if the monotonicity assumption (q′(x)>0) holds in accordance to Remark (2.5) and from (2.10) 1-Pr{q(x≤|1/κℒ,ℛ|)}≤γ where γ is a specified threshold probability (e.g. γ=0.001).Figures1 and 2 give examples of symmetric κ and asymmetric κℒ-κℛ distributions based on (2.9) and (2.10). Figures 1(a) and 1(b) are examples of symmetric global pdfs. Figure 2(a) is an example of an asymmetric global pdf and Figure 2(b) is an example of an asymmetric local pdf. The point where monotonicity fails for the local pdf in Figure 2(b) is located at q(x=1/-0.0497)=-7.402 where making use of (2.10) yields 1-Pr{q(x)≤|-7.402|}=0.0006, which would not pose a serious limitation for most purposes. In the next section the L-moments for κ and κℒ-κℛ distributions are derived and other distributional properties are discussed after a preliminary discussion of univariate L-moments.Two symmetricκ global pdfs and their associated cdfs. The percentiles were computed based on (2.10). The parameters of kurtosis (α4) and L-kurtosis (τ4) are based on (A.2) in the appendix and (3.11), respectively.
(a)
(b)Two asymmetricκℒ-κℛ distributions pdfs and their associated cdfs. The parameters of skew (α3) and kurtosis (α4) are based on (A.2) in the appendix. The parameters of L-skew (τ3) and L-kurtosis (τ4) are based on (3.15) and (3.16). The distribution in (b) is a local pdf in the range of q(x)∈(-7.402,+∞).
(a)
(b)


## 3.L-Moments for κ and κℒ-κℛ Distributions

### 3.1. Preliminaries: UnivariateL-Moments

LetY1,…,Yj,…,Yn be independent and identically distributed random variables each with continuous pdf fY(y), cdf FY(y), order statistics denoted as Y1:n≤⋯≤Yj:n≤⋯≤Yn: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 Yj:n are expressed as [8, pages 20–22]
(3.1)λ1=E[Y1:1]=β0,(3.2)λ2=12E[Y2:2-Y1:2]=2β1-β0,(3.3)λ3=13E[Y3:3-2Y2:3+Y1:3]=6β2-6β1+β0,(3.4)λ4=14E[Y4:4-3Y3:4+3Y2:4-Y1:4]=20β3-30β2+12β1-β0,
where the βi are determined from
(3.5)βi=∫y{FY(y)}ifY(y)dy,
where i=0,…,3. The coefficients associated with βi in (3.5) are obtained from shifted orthogonal Legendre polynomials and are computed as shown in [8, page 20] or in [11].TheL-moments λ1 and λ2 in (3.1) and (3.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 [20]:
(3.6)5τ32-14<τ4<1.


### 3.2.L-Moments for Symmetric κ and Asymmetric κℒ-κℛ Distributions

The derivation of the first fourL-moments associated with symmetric κ distributions begins by defining the probability weighted moments based on (3.5) in terms of q(x) in (2.4) and the standard logistic pdf and cdf in (2.1) and (2.2) as
(3.7)βi=∫-∞+∞q(x,κ){FX(x)}ifX(x)dx.
Integrating (3.7) for i=0,1,2,3 and using (3.1)–(3.4) gives λ1, λ2, τ3, and τ4 as
(3.8)λ1=0,(3.9)λ2=12(2+κ(-4A+4B+κ(C-D))),(3.10)τ3=0,(3.11)τ4=16(1+5κ2+10κ(1+2κ+2κ2(H[-(κ/2)]-H[-(1/2)-(κ/2)]))2λ2),
where A=ψ[0,(1/2)-(κ/2)],B=ψ[0,1-(κ/2)],C=ψ[1,(1/2)-(κ/2)], and D=ψ[1,1-(κ/2)]. The notations ψ[a,b] and H[·] in (3.9) and (3.11) are the polygamma and harmonic number functions, respectively. The argument for b must be positive in ψ[a,b] and thus a κ-distribution with a global pdf and defined values of λ1, λ2, τ3, τ4 requires that 0≤κ<1.The derivation of theL-moments for asymmetric κℒ-κℛ distributions associated with (2.6) begins with determining the probability-weighted moments βi in (3.5) by separately evaluating and summing two integrals as
(3.12)βi=∫-∞0q(x,κℒ){FX(x)}ifX(x)dx+∫0+∞q(x,κℛ){FX(x)}ifX(x)dx.As such, using (3.12) to obtain β0,…, β3 and subsequently substituting these terms into (3.1)–(3.4) we have
(3.13)λ1=14(2A-2B-2C+2D-κℒE+κℒF+κℛG-κℛH),(3.14)λ2=14(4+κℒ(-4A+4B+κℒ(E-F))+4+κℛ(-4C+4D+κℛ(G-H))),(3.15)τ3=(4κℛ-4κℒ+κℒ2(6A-6B+κℒ(F-E))+κℛ2(-6C+6D+κℛ(G-H)))(4λ2),(3.16)τ4=𝔚(24λ2),
where 𝔚  denotes {-4(κℛ+10κℛ3)C+5κℒ((4+40κℒ2)B+κℒ(1+5κℒ2)(E-F)+κℛ((4+40κℛ2)D+κℛ(1+5κℛ2)(G-H)))2(2+5κℒ(1+3κℒ)+5κℛ(1+3κℛ))-4(κℒ+10κℒ3)A-4(κℛ+10κℛ3)C+5κℒ((4+40κℒ2)B+κℒ(1+5κℒ2)(E-F)+κℛ((4+40κℛ2)D+κℛ(1+5κℛ2)(G-H)))4(κℛ+10κℛ3)C+5κℒ((4+40κℒ2)B+κℒ(1+5κℒ2)(E-F)+κℛ((4+40κℛ2)D+κℛ(1+5κℛ2)(G-H)))},  where the polygamma functions are A=ψ[0,(1/2)-(κℒ/2)], B=ψ[0,1-(κℒ/2)], C=ψ[0,(1/2)-(κℛ/2)], D=ψ[0,1-(κℛ/2)], E=ψ[1,(1/2)-(κℒ/2)],  F=ψ[1,1-(κℒ/2)], G=ψ[1,(1/2)-(κℛ/2)], and H=ψ[1,1-(κℛ/2)]. Thus, given specified values of τ3 and τ4, (3.15) and (3.16) can be numerically solved to obtain the corresponding values of κℒ and κℛ. Analogous to κ-distributions, an asymmetric κℒ-κℛ distribution with a global pdf and defined values of λ1, λ2, τ3, τ4 will require both 0≤κℒ<1 and 0≤κℛ<1. Note that careful inspection of (3.15) and (3.16) reveals that interchanging values for the parameters κℒ and κℛ reverse the direction of τ3 and has no effect on τ4. Further, if κℒ=κℛ, then equations (3.13)–(3.16) simplify to (3.8)–(3.11) associated with the symmetric κ family of distributions.In terms of boundary conditions, the minimum value ofL-kurtosis for global pdfs is min(τ4)=1/6 where κℒ=κℛ=κ=0, which is associated with the standard logistic pdf in (2.1). As such, using (3.15) and (3.16) with κℒ=0 and κℛ∈[0,1], provided in Figure 3 is a graph of the region for feasible combinations of |τ3| and τ4 for asymmetric κℒ-κℛ global pdfs. Feasible combinations of L-skew and L-kurtosis for (3.15) and (3.16) 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 for global pdfs. An asymmetric κℒ-κℛ distributions will lie in the area above the curve graphed in the |τ3| and τ4 plane.The conventional-moment-based system forκℒ-κℛ distributions is given in the appendix. This system was used to determine the values of skew and kurtosis associated with the distributions given in Figures 1 and 2. It is worthy to point out that the conventional moment based system has a disadvantage in terms of moment existence. That is, the integral in (A.1) of the appendix reveals that for the rth moment to exist we must have 0≤κℒ<1/r and 0≤κℛ<1/r to yield appropriate polygamma functions. More specifically, if the mean, variance, skew, and kurtosis exist, then the parameters must be bounded in the range of 0≤κℒ<0.25 and 0≤κℛ<0.25. The advantage that the L-moment system has in this context is attributed to Hosking’s Theorem 1 [7] which states that if the mean (λ1) exists, then all other L-moments will have finite expectations.Figure4 provides some additional examples of κ and κℒ-κℛ distributions. These distributions will be used in the simulation portion of this study in Section 5. Note that the distributions associated with Figures 4(a) and 4(b) are local pdfs with probabilities at the point where monotonicity fails of 1-Pr{q(x)≤|-8.2084|}=0.00027 and 1-Pr{q(x)≤|12.9285|}=2.43×10-6, which do not pose any limitation for the purposes considered herein.  In the next section we first introduce the topic of the L-correlation and subsequently develop the methodology for simulating κℒ-κℛ (or κ) distributions with specified L-correlations.Two asymmetricκℒ-κℛ distributions and two symmetric κ distribution with their L-moment and conventional parameters of L-skew (τ3) and skew (α3), L-kurtosis (τ4) and kurtosis (α4), and corresponding shape parameters for (2.5) and (2.7). Note that in Tables 1–11 (a), (b), (c), (d) are denoted as 1, 2, 3, 4, respectively.
(a)
(b)
(c)
(d)


## 3.1. Preliminaries: UnivariateL-Moments

LetY1,…,Yj,…,Yn be independent and identically distributed random variables each with continuous pdf fY(y), cdf FY(y), order statistics denoted as Y1:n≤⋯≤Yj:n≤⋯≤Yn: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 Yj:n are expressed as [8, pages 20–22]
(3.1)λ1=E[Y1:1]=β0,(3.2)λ2=12E[Y2:2-Y1:2]=2β1-β0,(3.3)λ3=13E[Y3:3-2Y2:3+Y1:3]=6β2-6β1+β0,(3.4)λ4=14E[Y4:4-3Y3:4+3Y2:4-Y1:4]=20β3-30β2+12β1-β0,
where the βi are determined from
(3.5)βi=∫y{FY(y)}ifY(y)dy,
where i=0,…,3. The coefficients associated with βi in (3.5) are obtained from shifted orthogonal Legendre polynomials and are computed as shown in [8, page 20] or in [11].TheL-moments λ1 and λ2 in (3.1) and (3.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 [20]:
(3.6)5τ32-14<τ4<1.


## 3.2.L-Moments for Symmetric κ and Asymmetric κℒ-κℛ Distributions

The derivation of the first fourL-moments associated with symmetric κ distributions begins by defining the probability weighted moments based on (3.5) in terms of q(x) in (2.4) and the standard logistic pdf and cdf in (2.1) and (2.2) as
(3.7)βi=∫-∞+∞q(x,κ){FX(x)}ifX(x)dx.
Integrating (3.7) for i=0,1,2,3 and using (3.1)–(3.4) gives λ1, λ2, τ3, and τ4 as
(3.8)λ1=0,(3.9)λ2=12(2+κ(-4A+4B+κ(C-D))),(3.10)τ3=0,(3.11)τ4=16(1+5κ2+10κ(1+2κ+2κ2(H[-(κ/2)]-H[-(1/2)-(κ/2)]))2λ2),
where A=ψ[0,(1/2)-(κ/2)],B=ψ[0,1-(κ/2)],C=ψ[1,(1/2)-(κ/2)], and D=ψ[1,1-(κ/2)]. The notations ψ[a,b] and H[·] in (3.9) and (3.11) are the polygamma and harmonic number functions, respectively. The argument for b must be positive in ψ[a,b] and thus a κ-distribution with a global pdf and defined values of λ1, λ2, τ3, τ4 requires that 0≤κ<1.The derivation of theL-moments for asymmetric κℒ-κℛ distributions associated with (2.6) begins with determining the probability-weighted moments βi in (3.5) by separately evaluating and summing two integrals as
(3.12)βi=∫-∞0q(x,κℒ){FX(x)}ifX(x)dx+∫0+∞q(x,κℛ){FX(x)}ifX(x)dx.As such, using (3.12) to obtain β0,…, β3 and subsequently substituting these terms into (3.1)–(3.4) we have
(3.13)λ1=14(2A-2B-2C+2D-κℒE+κℒF+κℛG-κℛH),(3.14)λ2=14(4+κℒ(-4A+4B+κℒ(E-F))+4+κℛ(-4C+4D+κℛ(G-H))),(3.15)τ3=(4κℛ-4κℒ+κℒ2(6A-6B+κℒ(F-E))+κℛ2(-6C+6D+κℛ(G-H)))(4λ2),(3.16)τ4=𝔚(24λ2),
where 𝔚  denotes {-4(κℛ+10κℛ3)C+5κℒ((4+40κℒ2)B+κℒ(1+5κℒ2)(E-F)+κℛ((4+40κℛ2)D+κℛ(1+5κℛ2)(G-H)))2(2+5κℒ(1+3κℒ)+5κℛ(1+3κℛ))-4(κℒ+10κℒ3)A-4(κℛ+10κℛ3)C+5κℒ((4+40κℒ2)B+κℒ(1+5κℒ2)(E-F)+κℛ((4+40κℛ2)D+κℛ(1+5κℛ2)(G-H)))4(κℛ+10κℛ3)C+5κℒ((4+40κℒ2)B+κℒ(1+5κℒ2)(E-F)+κℛ((4+40κℛ2)D+κℛ(1+5κℛ2)(G-H)))},  where the polygamma functions are A=ψ[0,(1/2)-(κℒ/2)], B=ψ[0,1-(κℒ/2)], C=ψ[0,(1/2)-(κℛ/2)], D=ψ[0,1-(κℛ/2)], E=ψ[1,(1/2)-(κℒ/2)],  F=ψ[1,1-(κℒ/2)], G=ψ[1,(1/2)-(κℛ/2)], and H=ψ[1,1-(κℛ/2)]. Thus, given specified values of τ3 and τ4, (3.15) and (3.16) can be numerically solved to obtain the corresponding values of κℒ and κℛ. Analogous to κ-distributions, an asymmetric κℒ-κℛ distribution with a global pdf and defined values of λ1, λ2, τ3, τ4 will require both 0≤κℒ<1 and 0≤κℛ<1. Note that careful inspection of (3.15) and (3.16) reveals that interchanging values for the parameters κℒ and κℛ reverse the direction of τ3 and has no effect on τ4. Further, if κℒ=κℛ, then equations (3.13)–(3.16) simplify to (3.8)–(3.11) associated with the symmetric κ family of distributions.In terms of boundary conditions, the minimum value ofL-kurtosis for global pdfs is min(τ4)=1/6 where κℒ=κℛ=κ=0, which is associated with the standard logistic pdf in (2.1). As such, using (3.15) and (3.16) with κℒ=0 and κℛ∈[0,1], provided in Figure 3 is a graph of the region for feasible combinations of |τ3| and τ4 for asymmetric κℒ-κℛ global pdfs. Feasible combinations of L-skew and L-kurtosis for (3.15) and (3.16) 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 for global pdfs. An asymmetric κℒ-κℛ distributions will lie in the area above the curve graphed in the |τ3| and τ4 plane.The conventional-moment-based system forκℒ-κℛ distributions is given in the appendix. This system was used to determine the values of skew and kurtosis associated with the distributions given in Figures 1 and 2. It is worthy to point out that the conventional moment based system has a disadvantage in terms of moment existence. That is, the integral in (A.1) of the appendix reveals that for the rth moment to exist we must have 0≤κℒ<1/r and 0≤κℛ<1/r to yield appropriate polygamma functions. More specifically, if the mean, variance, skew, and kurtosis exist, then the parameters must be bounded in the range of 0≤κℒ<0.25 and 0≤κℛ<0.25. The advantage that the L-moment system has in this context is attributed to Hosking’s Theorem 1 [7] which states that if the mean (λ1) exists, then all other L-moments will have finite expectations.Figure4 provides some additional examples of κ and κℒ-κℛ distributions. These distributions will be used in the simulation portion of this study in Section 5. Note that the distributions associated with Figures 4(a) and 4(b) are local pdfs with probabilities at the point where monotonicity fails of 1-Pr{q(x)≤|-8.2084|}=0.00027 and 1-Pr{q(x)≤|12.9285|}=2.43×10-6, which do not pose any limitation for the purposes considered herein.  In the next section we first introduce the topic of the L-correlation and subsequently develop the methodology for simulating κℒ-κℛ (or κ) distributions with specified L-correlations.Two asymmetricκℒ-κℛ distributions and two symmetric κ distribution with their L-moment and conventional parameters of L-skew (τ3) and skew (α3), L-kurtosis (τ4) and kurtosis (α4), and corresponding shape parameters for (2.5) and (2.7). Note that in Tables 1–11 (a), (b), (c), (d) are denoted as 1, 2, 3, 4, respectively.
(a)
(b)
(c)
(d)


## 4.L-Correlations for κℒ-κℛ (or κ) Distributions

### 4.1. Preliminaries: TheL-Correlation

The coefficient ofL-correlation [18] is introduced by considering two random variables Yj and Yk with distribution functions F(Yj) and F(Yk), respectively. The second L-moments of Yj and Yk can alternatively be expressed as
(4.1)λ2(Yj)=2Cov(Yj,F(Yj)),(4.2)λ2(Yk)=2Cov(Yk,F(Yk)).
The second L-comoments of Yj toward Yk and Yk toward Yj are
(4.3)λ2(Yj,Yk)=2Cov(Yj,F(Yk)),(4.4)λ2(Yk,Yj)=2Cov(Yk,F(Yj)).
As such, the L-correlations of Yj toward Yk and Yk toward Yj are expressed as
(4.5)ηjk=λ2(Yj,Yk)λ2(Yj)(4.6)ηkj=λ2(Yk,Yj)λ2(Yk).
The L-correlation in (4.5) or (4.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.


### 4.2. TheL-Correlation for κℒ-κℛ (or κ) Distributions

In the context ofL-moment-based κℒ-κℛ (or κ) distributions, suppose it is desired to simulate a T-variate distribution from quantile functions of the form in (2.7) 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:
(4.7)Φ(zj)=Pr{Zj≤zj}=∫-∞zj(2π)-1/2exp{-wj22}dwj,(4.8)Φ(zk)=Pr{Zk≤zk}=∫-∞zk(2π)-1/2exp{-wk22}dwk,(4.9)fjk=(2π(1-ρjk2)1/2)-1exp{-(2(1-ρjk2))-1(zj2+zk2-2ρjkzjzk)}.
Using (4.7), it follows that the jth κℒ-κℛ distribution associated with (2.7) can be expressed as qj(g(Φ(zj))) where g(Φ(zj))=ln(Φ(zj)/(1-Φ(zj))) is standard logistic because Φ(zj)~U(0,1). As such, using (4.5), the L-correlation of qj(g(Φ(zj))) toward qk(g(Φ(zk))) can be evaluated using solved values of κℒj and κℛj for qj(g(Φ(zj))), a specified intermediate correlation (IC) ρjk in (4.9), and the following integral expressed as
(4.10)ηjk=2π∬-∞+∞sj(qj(g(Φ(zj)),κℒj,κℛj))Φ(zk)fjkdzjdzk.
We would point out that the purpose of the IC (ρjk) in (4.9) and (4.10) is to adjust for the effect of the transformation qj(g(Φ(zj))), which is induced by the κℒj and κℛj parameters, such that qj(g(Φ(zj))) has its specified L-correlation (ηjk) toward qk(g(Φ(zk))). Further, to simplify the computation, the κℒj-κℛj distribution in (4.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
(4.11)sj(qj(g(Φ(zj)),κℒj,κℛj))=δ(qj(g(Φ(zj)),κℒj,κℛj)-λ1),
where λ1 is the mean from (3.13) and δ is a constant that scales λ2 in (3.14) and in the denominator of (4.5) to 1/π  as
(4.12)δ=4(π(8-4Aκℒj+4Bκℒj+Eκℒj2-Fκℒj2-4Cκℛj+4Dκℛj+Gκℛj2+Hκℛj2)).
Analogously, the L-correlation of qk(g(Φ(zk))) toward qj(g(Φ(zj))) is expressed as
(4.13)ηkj=2π∬-∞+∞sk(qk(g(Φ(zk)),κℒk,κℛk))Φ(zj)fjkdzkdzj.
Note also for the special case that if qj(g(Φ(zj))) in (4.10) and qk(g(Φ(zk)))  in (4.13) have the same parameters, that is, κℒj=κℒk and κℛj=κℛk, then ηjk=ηkj. Provided in Algorithm 1 is Mathematica [21] source code that implements the computation of an IC (ρjk) based on (4.10). The details for simulating κℒ-κℛ (or κ) distributions with specified values of L-skew, L-kurtosis, and L-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=1 towards distribution k=2 (η12) in Figure 4 (Panels 4(a) and 4(b)). See also Table 3(a).
(*Intermediate Correlation*)
ρ
j
k
=
0.678043;
Needs[“MultivariateStatistics’”]
f
j
k
= PDF[MultinormalDistribution[{0,0},{{1,ρjk},{ρjk,1}}], {Zj,Zk}];
F
j
=
CDF[NormalDistribution[0,1],Zj];
F
k
=
CDF[NormalDistribution[0,1],Zk];
(*Parameters for Distribution 4(a) in Figure4*)
κ
ℒ
=
-
0.044817709;
κ
ℛ
=
0.1704343967;
X
=
Log[Fj/(1-Fj)];
q
ℒ
=
X
*
Exp
[
κ
ℒ
*
Abs
[
X
]
];
q
ℛ
=
X
*
Exp
[
κ
ℛ
*
Abs
[
X
]
];
(* Standardizing constantsλ1 from (3.13) and  δ from (4.12)*)
S
q
ℒ
=
δ
*
(
q
ℒ
-
λ
1
);
S
q
ℛ
=
δ
*
(
q
ℛ
-
λ
1
);
(*Compute the specifiedL-correlation*)
η
j
k
=
2
π
*NIntegrate[Piecewise[{{Sqℒ,X≤0},{Sqℛ,X>0}}]*Fk*fjk, {Zj,-10,10},
{
Z
k
,
-
10
,
10
}, Method→MultiDimensional]
0.70


## 4.1. Preliminaries: TheL-Correlation

The coefficient ofL-correlation [18] is introduced by considering two random variables Yj and Yk with distribution functions F(Yj) and F(Yk), respectively. The second L-moments of Yj and Yk can alternatively be expressed as
(4.1)λ2(Yj)=2Cov(Yj,F(Yj)),(4.2)λ2(Yk)=2Cov(Yk,F(Yk)).
The second L-comoments of Yj toward Yk and Yk toward Yj are
(4.3)λ2(Yj,Yk)=2Cov(Yj,F(Yk)),(4.4)λ2(Yk,Yj)=2Cov(Yk,F(Yj)).
As such, the L-correlations of Yj toward Yk and Yk toward Yj are expressed as
(4.5)ηjk=λ2(Yj,Yk)λ2(Yj)(4.6)ηkj=λ2(Yk,Yj)λ2(Yk).
The L-correlation in (4.5) or (4.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.


## 4.2. TheL-Correlation for κℒ-κℛ (or κ) Distributions

In the context ofL-moment-based κℒ-κℛ (or κ) distributions, suppose it is desired to simulate a T-variate distribution from quantile functions of the form in (2.7) 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:
(4.7)Φ(zj)=Pr{Zj≤zj}=∫-∞zj(2π)-1/2exp{-wj22}dwj,(4.8)Φ(zk)=Pr{Zk≤zk}=∫-∞zk(2π)-1/2exp{-wk22}dwk,(4.9)fjk=(2π(1-ρjk2)1/2)-1exp{-(2(1-ρjk2))-1(zj2+zk2-2ρjkzjzk)}.
Using (4.7), it follows that the jth κℒ-κℛ distribution associated with (2.7) can be expressed as qj(g(Φ(zj))) where g(Φ(zj))=ln(Φ(zj)/(1-Φ(zj))) is standard logistic because Φ(zj)~U(0,1). As such, using (4.5), the L-correlation of qj(g(Φ(zj))) toward qk(g(Φ(zk))) can be evaluated using solved values of κℒj and κℛj for qj(g(Φ(zj))), a specified intermediate correlation (IC) ρjk in (4.9), and the following integral expressed as
(4.10)ηjk=2π∬-∞+∞sj(qj(g(Φ(zj)),κℒj,κℛj))Φ(zk)fjkdzjdzk.
We would point out that the purpose of the IC (ρjk) in (4.9) and (4.10) is to adjust for the effect of the transformation qj(g(Φ(zj))), which is induced by the κℒj and κℛj parameters, such that qj(g(Φ(zj))) has its specified L-correlation (ηjk) toward qk(g(Φ(zk))). Further, to simplify the computation, the κℒj-κℛj distribution in (4.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
(4.11)sj(qj(g(Φ(zj)),κℒj,κℛj))=δ(qj(g(Φ(zj)),κℒj,κℛj)-λ1),
where λ1 is the mean from (3.13) and δ is a constant that scales λ2 in (3.14) and in the denominator of (4.5) to 1/π  as
(4.12)δ=4(π(8-4Aκℒj+4Bκℒj+Eκℒj2-Fκℒj2-4Cκℛj+4Dκℛj+Gκℛj2+Hκℛj2)).
Analogously, the L-correlation of qk(g(Φ(zk))) toward qj(g(Φ(zj))) is expressed as
(4.13)ηkj=2π∬-∞+∞sk(qk(g(Φ(zk)),κℒk,κℛk))Φ(zj)fjkdzkdzj.
Note also for the special case that if qj(g(Φ(zj))) in (4.10) and qk(g(Φ(zk)))  in (4.13) have the same parameters, that is, κℒj=κℒk and κℛj=κℛk, then ηjk=ηkj. Provided in Algorithm 1 is Mathematica [21] source code that implements the computation of an IC (ρjk) based on (4.10). The details for simulating κℒ-κℛ (or κ) distributions with specified values of L-skew, L-kurtosis, and L-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=1 towards distribution k=2 (η12) in Figure 4 (Panels 4(a) and 4(b)). See also Table 3(a).
(*Intermediate Correlation*)
ρ
j
k
=
0.678043;
Needs[“MultivariateStatistics’”]
f
j
k
= PDF[MultinormalDistribution[{0,0},{{1,ρjk},{ρjk,1}}], {Zj,Zk}];
F
j
=
CDF[NormalDistribution[0,1],Zj];
F
k
=
CDF[NormalDistribution[0,1],Zk];
(*Parameters for Distribution 4(a) in Figure4*)
κ
ℒ
=
-
0.044817709;
κ
ℛ
=
0.1704343967;
X
=
Log[Fj/(1-Fj)];
q
ℒ
=
X
*
Exp
[
κ
ℒ
*
Abs
[
X
]
];
q
ℛ
=
X
*
Exp
[
κ
ℛ
*
Abs
[
X
]
];
(* Standardizing constantsλ1 from (3.13) and  δ from (4.12)*)
S
q
ℒ
=
δ
*
(
q
ℒ
-
λ
1
);
S
q
ℛ
=
δ
*
(
q
ℛ
-
λ
1
);
(*Compute the specifiedL-correlation*)
η
j
k
=
2
π
*NIntegrate[Piecewise[{{Sqℒ,X≤0},{Sqℛ,X>0}}]*Fk*fjk, {Zj,-10,10},
{
Z
k
,
-
10
,
10
}, Method→MultiDimensional]
0.70


## 5. The Method and Simulation Study

To implement the method for simulatingκℒ-κℛ (or κ) distributions with specified L-moments and L-correlations we suggest the following six steps.(1)
Specify theL-moments for T transformations of the form in (2.7), that is, q1(g(Φ(z1))),…,qT(g(Φ(zT))) and obtain the solutions for the parameters of κℒj and κℛjby solving (3.15) and (3.16) 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(g(Φ(zj))) toward qk(g(Φ(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 (4.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:
(5.1)Z1=a11V1,Z2=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 [22]:
(5.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 (5.2) is less than 8×10-16.(6)
Substitute the zero-one uniform deviates,Φ(Zj), generated from Step 5 into the T equations of the form of qj(g(Φ(zj))), where g(Φ(zj))=ln(Φ(zj)/(1-Φ(zj))) is standard logistic to generate the κℒ-κℛ distributions with the specified L-moments and L-correlations.To demonstrate the steps above and evaluate the proposed method, a comparison between the proposedL-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 methods, respectively. See Algorithm 2 for an example for computing ICs for the conventional method. 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 (5.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 (2.5) and (2.7) to produce q1(g(Φ(z1))),…,q4(g(Φ(z4))) for both methods.Table 1
Specified correlation matrices for the distributions in Figure4.
(a)
1
2
3
4
1
1
2
0.70
1
3
0.70
0.70
1
4
0.85
0.70
0.70
1
(b)
1
2
3
4
1
1
2
0.40
1
3
0.50
0.40
1
4
0.60
0.50
0.40
1Table 2
Intermediate correlations for the conventional moment procedure.
(a)
1
2
3
4
1
1
2
0.850758
1
3
0.781157
0.737229
1
4
0.942741
0.721594
0.720942
1
(b)
1
2
3
4
1
1
2
0.482879
1
3
0.571693
0.428504
1
4
0.673333
0.518208
0.415569
1Table 3
Intermediate correlations for theL-moment procedure.
(a)
1
2
3
4
1
1
2
0.678043
1
3
0.678043
0.687374
1
4
0.835100
0.687374
0.678876
1
(b)
1
2
3
4
1
1
2
0.380048
1
3
0.477524
0.388439
1
4
0.576716
0.487002
0.381013
1Table 4
Cholesky decompositions for the conventional moment procedure.
(a)
a
11
=
1
a
12
=
0.850758
a
13
=
0.781157
a
14
=
0.942741
0
a
22
=
0.525558
a
23
=
0.138241
a
24
=
-
0.153076
0
0
a
33
=
0.608838
a
34
=
0.00932
0
0
0
0
a
44
=
0.296176
(b)
a
11
=
1
a
12
=
0.482879
a
13
=
0.571693
a
14
=
0.673333
0
a
22
=
0.875687
a
23
=
0.174087
a
24
=
0.220478
0
0
a
33
=
0.801786
a
34
=
-
0.009670
0
0
0
a
44
=
0.705634Table 5
Cholesky decompositions for theL-moment procedure.
(a)
a
11
=
1
a
12
=
0.678043
a
13
=
0.678043
a
14
=
0.8351
00
0
a
22
=
0.735022
a
23
=
0.309694
a
24
=
0.164812
0
0
a
33
=
0.666594
a
34
=
0.092412
0
0
0
a
44
=
0.516629
(b)
a
11
=
1
a
12
=
0.380048
a
13
=
0.477524
a
14
=
0.576716
0
a
22
=
0.924967
a
23
=
0.223745
a
24
=
0.289548
0
0
a
33
=
0.849652
a
34
=
0.048058
0
0
0
a
44
=
0.762398Algorithm 2:Mathematica source code for computing intermediate correlations for specified conventional Pearson correlations. The example is for distributions j=1 and k=4 (ρ34*) in Figure 4 (Panels 4(a) and 4(d)). See also Table 2(b).
(*Intermediate Correlation*)
ρ
j
k
=
0.673333;
Needs[“MultivariateStatistics’”]
f
j
k
=
PDF[MultinormalDistribution[{0,0},{{1,ρjk},{ρjk,1}}],{Zj,Zk}];
F
j
=
CDF[NormalDistribution[0,1],Zj];
F
k
=
CDF[NormalDistribution[0,1],Zk];
(*Parameters for Distributions 4(a) and 4(d) in Figure4*)
κ
ℒ
j
=
-0.044817709;
κ
ℛ
j
=
0.1704343967;
κ
ℒ
k
=
0.0;
κ
ℛ
k
=
0.0;
X
j
=
Log[Fj/(1-Fj)];
X
k
=
Log[Fk/(1-Fk)];
q
ℒ
j
=
X
j
*
Exp
[
κ
ℒ
j
*
Abs
[
X
j
]
];
q
ℛ
j
=
X
j
*
Exp
[
κ
ℛ
j
*
Abs
[
X
j
]
];
q
ℒ
k
=
X
k
*
Exp
[
κ
ℒ
k
*
Abs
[
X
k
]
];
q
ℛ
k
=
X
k
*
Exp
[
κ
ℛ
k
*
Abs
[
X
k
]
];
(*Standardizing constantsλ1 from (3.13) and α2 from (A.2) from the appendix*)
S
q
ℒ
j
=
(
q
ℒ
j
-
λ
1
j
)
/
α
2
j;
S
q
ℛ
j
=
(
q
ℛ
j
-
λ
1
j
)
/
α
2
j;
S
q
ℒ
k
=
(
q
ℒ
k
-
λ
1
k
)
/
α
2
k;
S
q
ℛ
k
=
(
q
ℛ
k
-
λ
1
k
)
/
α
2
k;
(*Compute the specified conventional Pearson correlation*)
ρ
j
k
*
= NIntegrate[Piecewise[{{Sqℒj,Xj≤0},{Sqℛj,Xj>0}}]* Piecewise[{{Sqℒk,Xk≤0},
{
S
q
ℛ
k
,
X
k
>
0
}
}]*
f
j
k, {Zj,-10,10}, {Zk,-10,10}, Method→MultiDimensional]
0.60In terms of the simulation, a Fortran algorithm was written for both methods to generate 25,000 independent sample estimates for the specified parameters of: (a) conventional skew (α3), kurtosis (α4), and Pearson correlation (ρjk*) and (b) L-skew (τ3), L-kurtosis (τ4), and L-correlation (ηjk). All estimates were based on sample sizes of n=25 and n=1000. The formulae used for computing estimates of α3,4 were based on Fisher’s k-statistics, that is, the formulae currently used by most commercial software packages such as SAS, SPSS, Minitab, for computing indices of skew and kurtosis (where α3,4=0 for the standard normal distribution). The formulae used for computing estimates of τ3,4 were Headrick’s Equations 2.4 and 2.6 [11]. 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 (4.5) using the empirical forms of the cdfs in (4.1) and (4.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–11 and are discussed in the next section.Table 6
Skew (α3) and Kurtosis (α4) results for the conventional moment procedure.
(a)
n = 25
Dist
Parameter
Estimate
95% Bootstrap C.I.
Standard error
Relative bias %
1
α
3
=
4.502
1.259
1.246
,
1.272
0.0065
−72.0
α
4
=
129.5
3.286
3.230
,
3.340
0.0280
−97.5
2
α
3
=
-
1.356
−0.6655
-
0.6756
,
-
0.6551
0.0053
−50.9
α
4
=
7.518
1.657
1.625
,
1.694
0.0175
−78.0
3
α
3
=
0
0.0093
-
0.0052
,
0.0228
0.0071
—
α
4
=
9.674
2.374
2.335
,
2.411
0.0194
−75.5
4
α
3
=
0
0.0068
-
0.0016
,
0.0150
0.0042
—
α
4
=
1.2
0.7011
0.6826
,
0.7202
0.0096
−41.6
(b)
n = 1000
Dist
Parameter
Estimate
95% bootstrap C.I.
Standard error
Relative bias %
1
α
3
=
4.502
3.285
3.260
,
3.311
0.0128
−27.0
α
4
=
129.5
29.16
28.56
,
29.75
0.3057
−77.5
2
α
3
=
-
1.356
−1.286
-
1.293
,
-
1.279
0.0037
−5.2
α
4
=
7.518
6.260
6.157
,
6.373
0.0555
−16.7
3
α
3
=
0
−0.0011
-
0.0110
,
0.0090
0.0051
—
α
4
=
9.674
8.016
7.897
,
8.152
0.0644
−17.1
4
α
3
=
0
−0.0003
-
0.0019
,
0.0019
0.0010
—
α
4
=
1.2
1.179
1.173
,
1.185
0.0033
−1.79Table 7
L-skew (τ3) and L-kurtosis (τ4) results.
(a)
n
=
25
Dist
Parameter
Estimate
95% bootstrap C.I.
Stand. error
Relative bias %
1
τ
3
=
0.23
0.2043
0.2026
,
0.2061
0.00091
−11.2
τ
4
=
0.25
0.2343
0.2329
,
0.2357
0.00072
−6.3
2
τ
3
=
-
0.12
−0.1075
-
0
.
1090
,
-
0
.
1059
0.00086
−10.4
τ
4
=
0.20
0.1942
0.1931
,
0.1954
0.00058
−2.9
3
τ
3
=
0
0.0015
-
0.0004
,
0.0034
0.00098
—
τ
4
=
0.25
0.2413
0.2400
,
0.2425
0.00063
−3.5
4
τ
3
=
0
0.00116
-
0.00027
,
0.00249
0.00070
—
τ
4
=
0.1667
0.1654
0.1644
,
0.1664
0.00050
−0.76
(b)
n
=
1000
Dist
Parameter
Estimate
95% bootstrap C.I.
Stand. error
Relative bias %
1
τ
3
=
0.23
0.2291
0.2287
,
0.2294
0.00017
-
0.39
τ
4
=
0.25
0
.
2494
0
.
2492
,
0.2497
0.00013
-
0.24
2
τ
3
=
-
0.12
-
0.1197
-
0.1200
,
-
0.1195
0.00013
—
τ
4
=
0.20
0.1999
0.1997
,
0.2001
0.00009
—
3
τ
3
=
0
-
0.0001
-
0.0004
,
0.0003
0.00017
—
τ
4
=
0.25
0.2498
0.2496
,
0.2500
0.00010
—
4
τ
3
=
0
-
0.00004
-
0.00026
,
0.00016
0.00011
—
τ
4
=
0.1667
0.1667
0.1665
,
0.1668
0.00007
—Table 8
Correlation (strong) results for the conventional moment procedure.
(a)
n
=
25
Parameter
Estimate
95% bootstrap C.I.
Standard error
Relative bias%
ρ
12
*
=
0.70
0.7857
0.7846, 0.7863
0.00114
12.24
ρ
13
*
=
0.70
0.7528
0.7515, 0.7541
0.00148
7.54
ρ
14
*
=
0.85
0.9083
0.9079, 0.9088
0.00125
6.86
ρ
23
*
=
0.70
0.7277
0.7264, 0.7290
0.00141
3.96
ρ
24
*
=
0.70
0.7178
0.7165, 0.7190
0.00133
2.54
ρ
34
*
=
0.70
0.7185
0.7172, 0.7198
0.00138
2.64
(b)
n
=
1000
Parameter
Estimate
95% bootstrap C.I.
Standard error
Relative bias%
ρ
12
*
=
0.70
0.7145
0.7141, 0.7149
0.00044
2.07
ρ
13
*
=
0.70
0.7073
0.7069, 0.7077
0.00037
1.04
ρ
14
*
=
0.85
0.8615
0.8610, 0.8617
0.00067
1.35
ρ
23
*
=
0.70
0.7046
0.7043, 0.7048
0.00026
0.66
ρ
24
*
=
0.70
0.7009
0.7006, 0.7011
0.00022
0.13
ρ
34
*
=
0.70
0.7011
0.7009, 0.7013
0.00022
0.16Table 9
Correlation (strong) results for theL-moment procedure.
(a)
n
=
25
Parameter
Estimate
95% bootstrap C.I.
Standard error
Relative bias%
η
12
=
0.70
0.7069
0.7053, 0.7083
0.00153
0.99
η
13
=
0.70
0.7066
0.7050, 0.7081
0.00156
0.94
η
14
=
0.85
0.8543
0.8535, 0.8552
0.00158
0.51
η
23
=
0.70
0.7078
0.7063, 0.7092
0.00147
1.11
η
24
=
0.70
0.7079
0.7065, 0.7094
0.00147
1.13
η
34
=
0.70
0.7069
0.7054, 0.7083
0.00149
0.99
(b)
n
=
1000
Parameter
Estimate
95% bootstrap C.I.
Standard error
Relative bias%
η
12
=
0.70
0.6999
0.6997, 0.7002
0.00024
—
η
13
=
0.70
0.7002
0.6999, 0.7004
0.00024
—
η
14
=
0.85
0.8500
0.8499, 0.8502
0.00024
—
η
23
=
0.70
0.6999
0.6997, 0.7001
0.00022
—
η
24
=
0.70
0.7001
0.6998, 0.7003
0.00022
—
η
34
=
0.70
0.7001
0.6999, 0.7003
0.00022
—Table 10
Correlation (moderate) results for the conventional moment procedure.
(a)
n
=
25
Parameter
Estimate
95% bootstrap C.I.
Standard error
Relative bias%
ρ
12
*
=
0.40
0.4517
0.4498, 0.4536
0.00123
12.93
ρ
13
*
=
0.50
0.5492
0.5473, 0.5512
0.00143
9.84
ρ
14
*
=
0.60
0.6505
0.6490, 0.6520
0.00134
8.42
ρ
23
*
=
0.40
0.4226
0.4203, 0.4247
0.00137
5.65
ρ
24
*
=
0.50
0.5170
0.5150, 0.5189
0.00135
3.40
ρ
34
*
=
0.40
0.4155
0.4133, 0.4178
0.00136
3.88
(b)
n
=
1000
Parameter
Estimate
95% bootstrap C.I.
Standard error
Relative bias%
ρ
12
*
=
0.40
0.4076
0.4072, 0.4079
0.00023
1.90
ρ
13
*
=
0.50
0.5061
0.5057, 0.5066
0.00030
1.22
ρ
14
*
=
0.60
0.6071
0.6067, 0.6075
0.00031
1.18
ρ
23
*
=
0.40
0.4026
0.4023, 0.4030
0.00022
0.65
ρ
24
*
=
0.50
0.5014
0.5011, 0.5017
0.00021
0.28
ρ
34
*
=
0.40
0.4007
0.4004, 0.4011
0.00021
0.18Table 11
Correlation (moderate) results for theL-moment procedure.
(a)
n
=
25
Parameter
Estimate
95% bootstrap C.I.
Standard error
Relative bias%
η
12
=
0.40
0.4075
0.4048, 0.4100
0.00157
1.88
η
13
=
0.50
0.5060
0.5036, 0.5081
0.00156
1.20
η
14
=
0.60
0.6065
0.6044, 0.6083
0.00156
1.08
η
23
=
0.40
0.4071
0.4048, 0.4096
0.00147
1.78
η
24
=
0.50
0.5083
0.5060, 0.5103
0.00148
1.66
η
34
=
0.40
0.4068
0.4042, 0.4091
0.00149
1.70
(b)
n
=
1000
Parameter
Estimate
95% bootstrap C.I.
Standard error
Relative bias%
η
12
=
0.40
0.4000
0.3996, 0.4004
0.00024
—
η
13
=
0.50
0.5001
0.4997, 0.5005
0.00024
—
η
14
=
0.60
0.6002
0.5999, 0.6004
0.00024
—
η
23
=
0.40
0.3996
0.3993, 0.4000
0.00022
—
η
24
=
0.50
0.5001
0.4998, 0.5004
0.00022
—
η
34
=
0.40
0.4000
0.3996, 0.4003
0.00022
—


## 6. Discussion and Conclusion

Headrick [11] demonstrated the advantages that L-moment ratios have over conventional-moment-based estimators in terms of relative bias and efficiency when sampling was from power method distributions with either moderate or heavy tails. Inspection of the simulation results in Table 6 and Table 7 of this study clearly indicates that this is also the case for the family of κ and κℒ-κℛ distributions. 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 28% and 2.5% of their associated population parameters, whereas the estimates of L-skew and L-kurtosis were 88.8% and 93.7% 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 relative standard errors RSE = (standard error/estimate) × 100 are substantially smaller than the conventional estimators of skew and kurtosis. For example, in terms of Distribution 1 (n=1000), inspection of Tables 6 and 7 indicates RSE measures of: RSE(α^3) =0.390% and RSE(α^4) =1.05% compared with RSE(τ^3) =0.074% and RSE(τ^4) =0.052%. This demonstrates that L-skew and L-kurtosis have more precision because they have less variance around their estimates.Presented in Tables8–11 are the results associated with the conventional Pearson and L-correlations. Overall inspection of these tables indicates that the L-correlation is substantially superior to the Pearson correlation in terms of relative bias. For example, in terms of strong correlations (n=25) the relative bias for the two heavier tailed distributions (distributions 1 and 2) was 12.24% for the Pearson correlation compared to only 0.99% 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 method is an attractive alternative to the traditional conventional-moment-based procedure. In particular, the L-moment based procedure has distinct advantages when distributions with moderate and heavy tails are used. Finally, we would note that Mathematica Version 8.0 [21] source code is available from the authors for implementing the L-moment-based method.



---

*Source: 980827-2012-09-25.xml*