A Logistic 𝐿-Moment-Based Analog for the Tukey 𝑔-ℎ, 𝑔, ℎ, and ℎ-ℎ System of Distributions

Abstract

This paper introduces a standard logistic L-moment-based system of distributions. The proposed system is an analog to the standard normal conventional moment-based Tukey g-h, g, h, and h-h system of distributions. The system also consists of four classes of distributions and is referred to as (i) asymmetric 𝛾-𝜅, (ii) log-logistic 𝛾, (iii) symmetric 𝜅, and (iv) asymmetric 𝜅_ℒ-𝜅_ℛ. The system can be used in a variety of settings such as simulation or modeling events—most notably when heavy-tailed distributions are of interest. A procedure is also described for simulating 𝛾-𝜅, 𝛾, 𝜅, and 𝜅_ℒ-𝜅_ℛ distributions with specified L-moments and L-correlations. The Monte Carlo results presented in this study indicate that estimates of L-skew, L-kurtosis, and L-correlation associated with the 𝛾-𝜅, 𝛾, 𝜅, and 𝜅_ℒ-𝜅_ℛ distributions are substantially superior to their corresponding conventional product-moment estimators in terms of relative bias and relative standard error.

1. Introduction

Simulating or modeling phenomena that involve heavy-tailed distributions have increasingly become of interest in many areas. Some examples include investigations associated with modeling stream flow and flood frequency [1–3], regional frequency analysis [4], aircraft landing processes [5], disaster analysis [6], finance [7], signal and image processing [8], and latent traits in the social and behavioral sciences [9]. Further, it is also common practice for methodologists to investigate the type I error and power properties associated with inferential statistics using heavy-tailed distributions [10–14]. In many cases, these investigations may only require an elementary transformation to produce heavy-tailed distributions with specified values of conventional skew and kurtosis. A system of transformations that is particularly well suited for this task is the Tukey 𝑔-ℎ, 𝑔, ℎ, and ℎ-ℎ classes, which can be used for simulating or modeling univariate and multivariate heavy-tailed distributions [14–20].

The quantile function associated with 𝑔-ℎ distributions can be succinctly described as 𝑞(𝑍)=𝑔⁻¹𝑒^𝑔𝑍−1𝑒^ℎ𝑍2, (1.1) where 𝑍∼𝑁(0,1), the parameters 𝑔≠0 and ℎ≥0 control the skew and kurtosis of 𝑞(𝑍), and the quantile functions for the one-parameter log-normal 𝑔 or symmetric ℎ distributions can be obtained by taking the lim_ℎ→0𝑞(𝑍) or lim_𝑔→0𝑞(𝑍), respectively. The class of asymmetric ℎ-ℎ distributions can be obtained by a doubling technique of ℎ distributions that is described in [18, 21]. Hence, an attractive feature of the Tukey system (𝑔-ℎ, 𝑔, ℎ, ℎ-ℎ classes) is that it is computationally efficient because each class of distributions only requires the knowledge of one or two parameters and an algorithm that generates standard normal random deviates.

One problem that arises in the context of heavy-tailed distributions is that 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 [4, 22–25]. To obviate this problem, 𝐿-moment-based estimators such as 𝐿-skew and 𝐿-kurtosis were introduced to address the limitations associated with conventional estimators [22]. Specifically, some of the advantages that 𝐿-moments have over conventional moments are that they (i) exist whenever the mean of the distribution exists, (ii) are nearly unbiased for all sample sizes and distributions, and (iii) are more robust in the presence of outliers [4, 22–25].

Thus, it would be advantageous to have a 𝐿-moment-based Tukey system of transformations for the purpose of simulating or modeling univariate and multivariate heavy-tailed distributions. However, because of the complexities associated with (1.1) and the fact that the standard normal distribution function is not available in closed form, the derivation of a 𝐿-moment-based Tukey system would be problematic. We would note that the Tukey ℎ and ℎ-ℎ classes of distributions have been characterized in the context of 𝐿-moments and the 𝐿-correlation (see [21]).

In view of the above, the present aim is to introduce an 𝐿-moment Tukey system analog that is based on the standard logistic distribution. More specifically, the system consists of four quantile functions that produce continuous symmetric and asymmetric distributions with specified values of 𝐿-skew, 𝐿-kurtosis, and 𝐿-correlation. The four classes of distributions are referred to as (i) asymmetric 𝛾-𝜅, (ii) log-logistic 𝛾, (iii) symmetric 𝜅, and (iv) asymmetric 𝜅_ℒ-𝜅_ℛ. Some of the advantages that the new system has is that the class of log-logistic (𝛾) distributions has demonstrated to be more efficient than its log-normal (𝑔) counterpart in terms of computing its hazard function or when censored data are encountered in the context of survival analysis [26]. Note that the 𝜅 and 𝜅_ℒ-𝜅_ℛ classes of distributions were recently introduced in an article by Headrick and Pant [27].

The rest of the paper is outlined as follows. In Section 2, the cumulative distribution function and probability density function as well as other properties associated with the four classes of distributions are derived. A summary of univariate 𝐿-moment theory is also provided, and the derivations of the systems of equations for specifying values of 𝐿-skew and 𝐿-kurtosis for the four classes of distributions are subsequently provided. In Section 3, the coefficient of 𝐿-correlation is introduced, and the equations are developed for determining intermediate correlations for specified 𝐿-correlations associated with the four classes of distributions. In Section 4, the steps for implementing a simulation procedure are described. Numerical examples and the results of a simulation are also provided to confirm the derivations and compare the new methodology with its conventional moment-based counterpart. In Section 5, the results of the simulation are discussed.

2. Methodology

2.1. Definitions and Properties for the System of 𝛾-𝜅, 𝛾, 𝜅, and 𝜅_ℒ-𝜅_ℛ Distributions

The derivation of the cumulative distribution function (cdf), probability density function (pdf), and other properties associated with the system of 𝛾-𝜅, 𝛾, 𝜅, and 𝜅_ℒ-𝜅_ℛ distributions begins with the following definitions.

Definition 2.1. Let 𝑋 be a random variable that has a standard logistic distribution with cdf and pdf expressed as 𝐹_𝑋(𝑥)=Ψ(𝑥)=Pr()=1𝑋≤𝑥(1+𝑒^−𝑥),𝑓_𝑋(𝑥)(𝑥)=𝑒=𝜓^−𝑥(1+𝑒^−𝑥)²,−∞<𝑥<+∞. (2.1)(2.2)

Definition 2.2. Let the quantile function 𝑞(𝑋) associated with the system of 𝛾-𝜅, 𝛾, 𝜅, and 𝜅_ℒ-𝜅_ℛ distributions be defined as 𝑞(𝑋)=𝑞_𝛾,𝜅(𝑋)=𝛾⁻¹𝑒^𝛾𝑋−1𝑒^𝜅|𝑋|,𝑞(𝑋)=𝑞_𝛾,0(𝑋)=lim_𝜅→0𝑞_𝛾,𝜅(𝑋)=𝛾⁻¹𝑒^𝛾𝑋−1,𝑞(𝑋)=𝑞_0,𝜅(𝑋)=lim_𝛾→0𝑞_𝛾,𝜅(𝑋)=𝑋𝑒^𝜅|𝑋|,𝑞(𝑋)=𝑞_{𝜅ℒ,𝜅ℛ}(𝑋)=𝑋𝑒^{𝜅ℒ|𝑋|},for𝑋≤0,𝑋𝑒^{𝜅ℛ|𝑋|},for𝑋≥0, (2.3)(2.4)(2.5)(2.6) where 𝑞(𝑋) in (2.3)–(2.6) is a strictly monotone increasing function with real-valued parameters 𝛾, 𝜅, 𝜅_ℒ, 𝜅_ℛ. Equations (2.3)–(2.6) produce distributions defined as (i) asymmetric 𝛾-𝜅 (𝛾≠0, 𝜅≥0), (ii) log-logistic 𝛾 (𝛾≠0), (iii) symmetric 𝜅 (𝜅≥0), and (iv) asymmetric 𝜅_ℒ-𝜅_ℛ (𝜅_ℒ≥0, 𝜅_ℛ≥0, 𝜅_ℒ≠𝜅_ℛ). The parameter ±𝛾 in (2.3) and (2.4) controls the degree of skew associated with a distribution. Taking the negative of 𝛾 will change the direction of the skew but not its magnitude that is, 𝑞_−𝛾,𝜅(𝑋)=−𝑞_𝛾,𝜅(−𝑋). The parameter 𝜅 in (2.3) and (2.5) controls the tail weight of a distribution where the function 𝑒^𝜅|𝑋| (i) preserves symmetry, (ii) is increasing for 𝑋≥0 and 𝜅≥0, and (iii) produces increased tail weight as the value of 𝜅 becomes larger. Analogous to 𝜅, the real-valued parameter 𝜅_ℒ (𝜅_ℛ) in (2.6) controls the left (right) tail weight of a 𝜅_ℒ-𝜅_ℛ distribution.

Theorem 2.3. The cdf and pdf associated with the 𝛾-𝜅, 𝛾, 𝜅, and 𝜅_ℒ-𝜅_ℛ classes of distributions in (2.3)–(2.6) are 𝐹_𝑞(𝑋)(𝑥)=Ψ(𝑥),𝑓_𝑞(𝑋)(𝑥)=𝜓(𝑥)𝑞^′(𝑥) (2.7)(2.8) with derivatives 𝑞^′(𝑥) expressed as 𝑞^′(𝑥)=𝑞^′_𝛾,𝜅(𝑥)=𝑒^{𝛾𝑥+𝜅|𝑥|}+𝛾⁻¹𝑒^𝜅|𝑥|𝑒^𝛾𝑥𝜅𝑑|𝑥|−1,𝑞𝑑𝑥_′(𝑥)=𝑞^′_𝛾,0(𝑥)=𝑒^𝛾𝑥,𝑞^′(𝑥)=𝑞^′_0,𝜅(𝑥)=𝑒^𝜅|𝑥|+𝑥𝜅𝑒^𝜅|𝑥|𝑑|𝑥|,𝑞𝑑𝑥_′(𝑥)=𝑞^′_{𝜅ℒ,𝜅ℛ}(𝑥)=𝑒^{𝜅ℒ|𝑥|}−𝑥𝜅_ℒ𝑒^{𝜅ℒ|𝑥|},for𝑒𝑥≤0,_{𝜅ℛ|𝑥|}+𝑥𝜅_ℛ𝑒^{𝜅ℛ|𝑥|},for𝑥≥0, (2.9) where 𝑞^′(𝑥)>0 and 𝑞^′(𝑥=0)=1 for (2.9).

Proof. The requirement that 𝑞(𝑋) in (2.3)–(2.6) is a strictly monotone increasing function implies that an inverse function (𝑞⁻¹) exists and thus 𝑦=𝐹_𝑞(𝑋)(𝑥)=Ψ(𝑥). Differentiating both sides with respect to 𝑞(𝑥) yields 𝑑𝑦/𝑑𝑞(𝑥)=𝑓_𝑞(𝑥)(𝑥). Hence, 𝑑𝑦/𝑑𝑞(𝑥)=(𝑑𝑦/𝑑𝑥)/(𝑑𝑞(𝑥)/𝑑𝑥)=𝜓(𝑥)/𝑞^′(𝑥). Whence, the pdf in (2.8) integrates to one because 𝜓(𝑥) is the logistic pdf in (2.2) and will be nonnegative on its support 𝑞(𝑥) for 𝑥∈(−∞,+∞) given from Definition 2.2 that 𝛾≠0, 𝜅≥0, 𝜅_ℒ≥0, 𝜅_ℛ≥0, and where lim_𝑥→±∞𝜓(𝑥)/𝑞^′(𝑥)=0 as lim_𝑥→±∞𝜓(𝑥)=0 and lim_𝑥→±∞1/𝑞^′(𝑥)=0.

Definition 2.4. If the monotonicity assumption (𝑞^′(𝑥)>0) holds for the pdf in (2.8) for all 𝑥∈(−∞,+∞), then (2.8) is defined to be a global pdf.

Remark 2.5. Inspection of (2.2) and (2.3)–(2.9) indicates that the height of a global pdf in (2.8) for any 𝛾-𝜅, 𝛾, 𝜅, or 𝜅_ℒ-𝜅_ℛ distribution at 𝑥=0 will be 𝜓(𝑥)=1/4.

Remark 2.6. The mode associated with a global pdf in (2.8) is located at 𝑞(𝑥) where 𝑥=𝑥 is the critical number that solves 𝑑(𝜓(𝑥)/𝑞^′(𝑥))/𝑑𝑞(𝑥)=0 and globally maximizes 𝜓(𝑥)/𝑞^′(𝑥) at the mode 𝑞(𝑥). It is noted that a global pdf in (2.8) will have a global maximum because the standard logistic pdf in (2.2) has a global maximum, and the transformations in (2.3)–(2.6) are assumed to be strictly monotone increasing functions.

Remark 2.7. The median associated with a global pdf in (2.8) is located at 𝑞(𝑥=0)=0. This can be shown by letting Ψ(𝑥)=0.50 denote the 50th percentile. In general, we must have 𝑥=0 such that Ψ(0)=Pr{𝑋≤0}=0.50 holds in (2.7) for the standard logistic distribution. As such, the limit of the quantile function 𝑞(𝑥) locates the median at lim_𝑥→0𝑞(𝑥)=0 for the system of 𝛾-𝜅, 𝛾, 𝜅, and 𝜅_ℒ-𝜅_ℛ distributions.

Remark 2.8. In terms of the 𝜅_ℒ-𝜅_ℛ (or 𝜅) class of distributions, the monotonicity assumption (𝑞^′(𝑥)>0) holds locally in (2.8) for cases where 𝜅_ℒ,ℛ<0 and (1/𝜅_ℒ<𝑥<0 or 0≤𝑥<−1/𝜅_ℛ). This leads to Definition 2.9.

Definition 2.9. A 𝜅_ℒ-𝜅_ℛ (or 𝜅) distribution is defined to have a local pdf if the monotonicity assumption holds in accordance with Remark 2.8 and 1−Pr{𝑞(𝑥≤|1/𝜅_ℒ,ℛ|)}≤𝜀, where 𝜀 is a user specified threshold probability based on the cdf in (2.7) (e.g., 𝜀≤0.001).
Examples of 𝛾-𝜅, 𝛾, 𝜅, and 𝜅_ℒ-𝜅_ℛ distributions based on the cdf and pdf in (2.7) and (2.8) are presented in Figure 1 through Figure 4, respectively. The heavy-tailed distributions in Figures 1(a), 2(a), 3(b), and 4(a) are used in the simulation portion of this study in Section 4. In the next section, univariate 𝐿-moments are introduced, and the 𝐿-moments for the system of 𝛾-𝜅, 𝛾, 𝜅, and 𝜅_ℒ-𝜅_ℛ distributions are subsequently derived, and other properties are discussed.

(a)

(b)

Figure 1 

Two asymmetric 𝛾-𝜅 pdfs and their cdfs. The parameters of skew (𝛼3) and kurtosis (𝛼4) are based on (A.2) in Appendix A. The parameters of 𝐿-skew (𝜏3) and 𝐿-kurtosis (𝜏4) are based on (2.19) and (2.20).

(a)

(b)

Figure 2 

Two log-logistic 𝛾 pdfs and their cdfs. The parameters of skew (𝛼3) and kurtosis (𝛼4) are based on (A.2) in Appendix A. The parameters of 𝐿-skew (𝜏3) and 𝐿-kurtosis (𝜏4) are based on (2.23) and (2.24).

(a)

(b)

Figure 3 

Two symmetric 𝜅 pdfs and their cdfs. The parameters of skew (𝛼3) and kurtosis (𝛼4) are based on (A.2) in Appendix A. The parameters of 𝐿-skew (𝜏3) and 𝐿-kurtosis (𝜏4) are based on (2.27) and (2.28).

(a)

(b)

Figure 4 

Two asymmetric 𝜅ℒ-𝜅ℛ distributions pdfs and their cdfs. The parameters of skew (𝛼3) and kurtosis (𝛼4) are based on (A.2) in Appendix A. The parameters of 𝐿-skew (𝜏3) and 𝐿-kurtosis (𝜏4) are based on (2.32) and (2.33). The distribution in (b) is a local pdf in the range of 𝑞(𝑥)∈(−8.085,+∞) where (2.7) yields 1−Pr{𝑞(𝑥)≤|−8.085|}=0.0003.

2.2. Preliminaries on Univariate 𝐿-Moments

Let 𝑌₁,…,𝑌_𝑗,…,𝑌_𝑛 be independent and identically distributed random variables each with continuous pdf 𝑓_𝑌(𝑦), cdf 𝐹_𝑌(𝑦), order statistics denoted as 𝑌_1∶𝑛≤⋯≤𝑌_𝑗∶𝑛≤⋯≤𝑌_𝑛∶𝑛, and 𝐿-moments defined in terms of either linear combinations of (i) expectations of order statistics or (ii) probability weighted moments (𝛽_𝑖). For the purposes considered herein, the first four 𝐿-moments associated with 𝑌_𝑗∶𝑛 are expressed as [4, pages 20–22] 𝜆₁=𝐸𝑌_1∶1=𝛽₀,𝜆₂=12𝐸𝑌_2∶2−𝑌_1∶2=2𝛽₁−𝛽₀,𝜆₃=13𝐸𝑌_3∶3−2𝑌_2∶3+𝑌_1∶3=6𝛽₂−6𝛽₁+𝛽₀,𝜆₄=14𝐸𝑌_4∶4−3𝑌_3∶4+3𝑌_2∶4−𝑌_1∶4=20𝛽₃−30𝛽₂+12𝛽₁−𝛽₀, (2.10)(2.11)(2.12)(2.13) where the 𝛽_𝑖 are determined from 𝛽_𝑖=𝑦𝐹_𝑌𝑦^𝑖𝑓_𝑌𝑦𝑑𝑦, (2.14) where 𝑖=0,…,3. The coefficients associated with 𝛽_𝑖 in (2.14) are obtained from shifted orthogonal Legendre polynomials and are computed as shown in [4, page 20] or in [22].

The 𝐿-moments 𝜆₁ and 𝜆₂ in (2.10) and (2.11) are measures of location and scale and are the arithmetic mean and one-half the coefficient of mean difference, respectively. Higher order 𝐿-moments are transformed to dimensionless quantities referred to as 𝐿-moment ratios defined as 𝜏_𝑟=𝜆_𝑟/𝜆₂ for 𝑟≥3, and where 𝜏₃ and 𝜏₄ are the analogs to the conventional measures of skew and kurtosis. In general, 𝐿-moment ratios are bounded in the interval −1<𝜏_𝑟<1 as is the index of 𝐿-skew (𝜏₃) where a symmetric distribution implies that all 𝐿-moment ratios with odd subscripts are zero. Other smaller boundaries can be found for more specific cases. For example, the index of 𝐿-kurtosis (𝜏₄) has the boundary condition for continuous distributions of [28] 5𝜏²₃−14<𝜏₄<1. (2.15)

2.3. 𝐿-Moments for the System of 𝛾-𝜅, 𝛾, 𝜅, and 𝜅_ℒ-𝜅_ℛ Distributions

The derivation of the first four 𝐿-moments associated with the 𝛾-𝜅 class of distributions begins by defining the probability weighted moments based on (2.14) in terms of (2.1)–(2.3) as 𝛽_𝑖=^+∞_−∞𝑞𝑥,𝛾,𝜅{Ψ(𝑥)}^𝑖𝜓(𝑥)𝑑𝑥. (2.16)

In general, to obtain finite values of 𝜆₁, 𝜆₂, 𝜏₃, 𝜏₄ based on (2.16), we must have (i) 𝛾+𝜅<1, (ii) 𝜅<1, and (iii) 1+𝛾>𝜅. As such, integrating (2.16) for 𝑖=0,1,2,3 and using (2.10)–(2.13) gives 𝜆₁, 𝜆₂, 𝜏₃, and 𝜏₄ as 𝜆₁=−𝛾−𝜅ℎ₁+𝛾−𝜅ℎ₂+−𝛾+𝜅ℎ₃+2𝜅ℎ₄+𝛾+𝜅ℎ₅−2𝜅ℎ₆,2𝛾𝜆₂=2𝛾−𝛾+𝜅²ℎ₁+𝛾−𝜅²ℎ₂−ℎ₃+𝛾+𝜅²ℎ₅,2𝛾𝜏₃=2𝛾²−𝛾+𝜅³ℎ₁+𝛾−𝜅³ℎ₂−𝛾−𝜅³ℎ₃+2𝜅³ℎ₄+𝛾+𝜅³ℎ₅−2𝜅³ℎ₆2𝛾𝜆₂,𝜏₄=2𝛾1+5𝛾²+5𝜅(1−3𝜅)−𝛾+𝜅²1+5𝛾+𝜅²ℎ₁+1+5𝛾−𝜅²𝛾−𝜅²ℎ₂−ℎ₃+𝛾+𝜅²1+5𝛾+𝜅²ℎ₅/12𝛾𝜆₂ (2.17)(2.18)(2.19)(2.20) with harmonic number functions 𝐻[⋅] of ℎ₁=𝐻[(1/2)(−1−𝛾−𝜅)], ℎ₂=𝐻[(1/2)(−1+𝛾−𝜅)], ℎ₃=𝐻[(1/2)(𝛾−𝜅)], ℎ₄=𝐻[(1/2)(−1−𝜅)], ℎ₅=𝐻[(−1/2)(𝛾+𝜅)], and ℎ₆=𝐻[−(1/2)𝜅]. As such, given user specied values of 𝜏₃ and 𝜏₄, (2.19) and (2.20) can be numerically solved to obtain the parameters for 𝛾 and 𝜅.

Taking the limit of (2.17)–(2.20) as 𝜅→0 gives the class associated with the log-logistic 𝛾 distributions in (2.4) as 𝜆₁=𝜋csc−1𝛾𝜋𝛾,𝜆₂=𝛾𝜋csc𝛾𝜋,𝜏₃=𝛾,𝜏₄=161+5𝛾². (2.21)(2.22)(2.23)(2.24)

Analogously, taking the limit of (2.17)–(2.20) as 𝛾→0 gives the system associated with the symmetric 𝜅 family of distributions in (2.5) as in [27] 𝜆₁=0,𝜆₂=122+𝜅−4𝑝₁+4𝑝₂+𝜅𝑝₃−𝑝₄,𝜏₃=0,𝜏₄=161+5𝜅²+10𝜅1+2𝜅+2𝜅²ℎ₆−ℎ₄2𝜆₂ (2.25)(2.26)(2.27)(2.28) with polygamma functions 𝑃[⋅] of 𝑝₁=𝑃[0,1/2−𝜅/2], 𝑝₂=𝑃[0,1−𝜅/2], 𝑝₃=𝑃[1,1/2−𝜅/2], and 𝑝₄=𝑃[1,1−𝜅/2].

The 𝐿-moments for the class of asymmetric 𝜅_ℒ-𝜅_ℛ distributions in (2.6) can be determined by separately evaluating and summing two integrals of the form in (2.16) as 𝛽_𝑖=⁰_−∞𝑞𝑥,𝜅_ℒ{Ψ(𝑥)}^𝑖𝜓(𝑥)𝑑𝑥+₀^+∞𝑞𝑥,𝜅_ℛ{Ψ(𝑥)}^𝑖𝜓(𝑥)𝑑𝑥. (2.29)

For an asymmetric 𝜅_ℒ-𝜅_ℛ distribution with a global pdf to have finite values of 𝜆₁, 𝜆₂, 𝜏₃, 𝜏₄ based on (2.29) will require 𝜅_ℒ<1 and 𝜅_ℛ<1. As such, integrating (2.29) to obtain 𝛽₀,…,𝛽₃ and subsequently substituting these terms into (2.10)–(2.13) yield as in [27] 𝜆₁=142𝑝₅−2𝑝₆−2𝑝₇+2𝑝₈−𝜅_ℒ𝑝₉+𝜅_ℒ𝑝₁₀+𝜅_ℛ𝑝₁₁−𝜅_ℛ𝑝₁₂,𝜆₂=144+𝜅_ℒ−4𝑝₅+4𝑝₆+𝜅_ℒ𝑝₉−𝑝₁₀+4+𝜅_ℛ−4𝑝₇+4𝑝₈+𝜅_ℛ𝑝₁₁−𝑝₁₂,𝜏₃=4𝜅_ℛ−4𝜅_ℒ+𝜅²_ℒ6𝑝₅−6𝑝₆+𝜅_ℒ𝑝₁₀−𝑝₉+𝜅²_ℛ−6𝑝₇+6𝑝₈+𝜅_ℛ𝑝₁₁−𝑝₁₂4𝜆₂,𝜏₄=22+5𝜅_ℒ1+3𝜅_ℒ+5𝜅_ℛ1+3𝜅_ℛ−4𝜅_ℒ+10𝜅³_ℒ𝑝₅−4𝜅_ℛ+10𝜅³_ℛ𝑝₇+5𝜅_ℒ4+40𝜅²_ℒ𝑝₆+𝜅_ℒ1+5𝜅²_ℒ𝑝₉−𝑝₁₀+𝜅_ℛ4+40𝜅²_ℛ𝑝₈+𝜅_ℛ1+5𝜅²_ℛ𝑝₁₁−𝑝₁₂/24𝜆₂ (2.30)(2.31)(2.32)(2.33) with polygamma functions 𝑃[⋅] of 𝑝₅=𝑃[0,1/2−𝜅_ℒ/2], 𝑝₆=𝑃[0,1−𝜅_ℒ/2], 𝑝₇=𝑃[0,1/2−𝜅_ℛ/2], 𝑝₈=𝑃[0,1−𝜅_ℛ/2], 𝑝₉=𝑃[1,1/2−𝜅_ℒ/2], 𝑝₁₀=𝑃[1,1−𝜅_ℒ/2], 𝑝₁₁=𝑃[1,1/2−𝜅_ℛ/2], 𝑝₁₂=𝑃[1,1−𝜅_ℛ/2].

As with the 𝛾-𝜅 class of distributions, given specified values of 𝜏₃ and 𝜏₄, (2.32) and (2.33) can be numerically solved to obtain the corresponding values of 𝜅_ℒ and 𝜅_ℛ. Inspection of (2.32) and (2.33) reveals that interchanging values for the parameters of 𝜅_ℒ and 𝜅_ℛ reverses the direction of 𝜏₃ and has no effect on 𝜏₄. As such, a graph of the region for feasible combinations of |𝜏₃| and 𝜏₄ for (2.32) and (2.33) is provided in Figure 5. Feasible combinations of 𝐿-skew and 𝐿-kurtosis will lie in the region above the curve graphed in the |𝜏₃|, 𝜏₄ plane of Figure 5. Note that the curve in Figure 5 was graphed by setting 𝜅_ℒ=0 with 𝜅_ℛ∈[0,1] in (2.32) and (2.33).

Figure 5 

Graph of the region for feasible combinations of (absolute value) 𝐿-skew |𝜏3| and 𝐿-kurtosis 𝜏4 for global 𝜅ℒ-𝜅ℛ pdfs. Asymmetric 𝜅ℒ-𝜅ℛ distributions will lie in the area above the curve graphed in the |𝜏3| and 𝜏4 plane.

The conventional moment-based systems for the 𝛾-𝜅 and 𝜅_ℒ-𝜅_ℛ classes of distributions are given in Appendices A and B, respectively. These systems were used to determine the values of skew and kurtosis associated with the distributions given in Figures 1–4. It is worthy to point out that the conventional moment-based systems have a disadvantage in terms of moment existence. That is, the integral in (A.1) of Appendix A reveals that for the 𝑟th moment to exist we must have in general (i) 𝛾+𝜅<1/𝑟, (ii) 𝜅<1/𝑟 and (iii) 1/𝑟+𝛾<𝜅. For example, if the mean, variance, skew, and kurtosis exist, then we must have 1/𝑟=0.25. And, analogously for the 𝜅_ℒ-𝜅_ℛ class of distributions in Appendix B, the parameters are bounded in the range of 0≤𝜅_ℒ<0.25 and 0≤𝜅_ℛ<0.25. The advantage that the 𝐿-moment system has in this context is attributed to Hosking’s Theorem 1 [22] which states that if the mean (𝜆₁) exists, then all other 𝐿-moments will have finite expectations. We note that the conventional moment-based systems for the log-logistic 𝛾 and symmetric 𝜅 classes of distributions can be obtained by simplifying the system as described in Appendix A. In the next section we first introduce the topic of the 𝐿-correlation and subsequently develop the methodology for simulating 𝛾-𝜅, 𝛾, 𝜅, and 𝜅_ℒ-𝜅_ℛ distributions with specified 𝐿-correlations.

3. 𝐿-Correlations for the System of 𝛾-𝜅, 𝛾, 𝜅, and 𝜅_ℒ-𝜅_ℛ Distributions

The coefficient of 𝐿-correlation [29] is introduced by considering two random variables 𝑌_𝑗 and 𝑌_𝑘 with distribution functions 𝐹(𝑌_𝑗) and 𝐹(𝑌_𝑘), respectively. The second 𝐿-moments of 𝑌_𝑗 and 𝑌_𝑘 can alternatively be expressed as 𝜆₂𝑌_𝑗=2Cov𝑌_𝑗,𝐹𝑌_𝑗,𝜆₂𝑌_𝑘=2Cov𝑌_𝑘,𝐹𝑌_𝑘. (3.1)(3.2)

The second 𝐿-comoments of 𝑌_𝑗 toward 𝑌_𝑘 and 𝑌_𝑘 toward 𝑌_𝑗 are 𝜆₂𝑌_𝑗,𝑌_𝑘=2Cov𝑌_𝑗,𝐹𝑌_𝑘,𝜆₂𝑌_𝑘,𝑌_𝑗=2Cov𝑌_𝑘,𝐹𝑌_𝑗. (3.3)(3.4)

As such, the 𝐿-correlations of 𝑌_𝑗 toward 𝑌_𝑘 and 𝑌_𝑘 toward 𝑌_𝑗 are expressed as 𝜂_𝑗𝑘=𝜆₂𝑌_𝑗,𝑌_𝑘𝜆₂𝑌_𝑗,𝜂_𝑘𝑗=𝜆₂𝑌_𝑘,𝑌_𝑗𝜆₂𝑌_𝑘. (3.5)(3.6)

The 𝐿-correlation in (3.5) (or (3.6)) is bounded such that −1≤𝜂_𝑗𝑘≤1 where a value of 𝜂_𝑗𝑘=1 (𝜂_𝑗𝑘=−1) indicates a strictly increasing (decreasing) monotone relationship between the two variables. In general, we would also note that 𝜂_𝑗𝑘≠𝜂_𝑘𝑗.

In the context of the 𝐿-moment-based 𝛾-𝜅, 𝛾, 𝜅, and 𝜅_ℒ-𝜅_ℛ classes of distributions, suppose that it is desired to simulate a 𝑇-variate distribution based on quantile functions of the forms in (2.3)–(2.6) with a specified 𝐿-correlation matrix and where each distribution has its own specified values of 𝜏₃ and 𝜏₄. Let 𝑍₁,…,𝑍_𝑇 denote standard normal variables where the distribution functions and bivariate density function associated with 𝑍_𝑗 and 𝑍_𝑘 are expressed as Φ𝑧_𝑗=Pr𝑍_𝑗≤𝑧_𝑗=^𝑧𝑗_−∞(2𝜋)^−1/2exp−𝑤²_𝑗2𝑑𝑤_𝑗,Φ𝑧_𝑘=Pr𝑍_𝑘≤𝑧_𝑘=^𝑧𝑘_−∞(2𝜋)^−1/2exp−𝑤²_𝑘2𝑑𝑤_𝑘,𝑓_𝑗𝑘=2𝜋1−𝜌²_𝑗𝑘^1/2−1exp−21−𝜌²_𝑗𝑘⁻¹𝑧²_𝑗+𝑧²_𝑘−2𝜌_𝑗𝑘𝑧_𝑗𝑧_𝑘. (3.7)(3.8)(3.9)

Using (3.7), it follows that the 𝑗th distribution associated with (2.3)–(2.6) can be expressed as 𝑞_𝑗(𝑔(Φ(𝑍_𝑗))), where 𝑔(Φ(𝑍_𝑗))=ln(Φ(𝑍_𝑗)/(1−Φ(𝑍_𝑗))) is standard logistic because Φ(𝑍_𝑗)∼𝑈(0,1). As such, using (3.5), the 𝐿-correlation of 𝑞_𝑗(𝑔(Φ(𝑍_𝑗))) toward 𝑞_𝑘(𝑔(Φ(𝑍_𝑘))) can be evaluated using the solved value(s) of the parameter(s) (i.e., 𝛾_𝑗-𝜅_𝑗; 𝛾_𝑗; 𝜅_𝑗; 𝜅_ℒ𝑗-𝜅_ℛ𝑗) for 𝑞_𝑗(𝑔(Φ(𝑍_𝑗))), a specified intermediate correlation (IC) 𝜌_𝑗𝑘 in (3.9), and the following integral generally expressed as 𝜂_𝑗𝑘=2√𝜋^+∞_−∞𝑠_𝑗𝑞_𝑗𝑔Φ𝑧_𝑗,𝛾_𝑗,𝜅_𝑗,𝜅_ℒ𝑗,𝜅_ℛ𝑗Φ𝑧_𝑘𝑓_𝑗𝑘𝑑𝑧_𝑗𝑑𝑧_𝑘. (3.10)

We would point out that the purpose of the IC (𝜌_𝑗𝑘) in (3.9) and (3.10) is to adjust for the effect of the transformation 𝑞_𝑗(𝑔(Φ(𝑍_𝑗))), which is induced by the parameters, such that 𝑞_𝑗(𝑔(Φ(𝑍_𝑗))) has its specified 𝐿-correlation (𝜂_𝑗𝑘) toward 𝑞_𝑘(𝑔(Φ(𝑍_𝑘))). Further, to simplify the computation, the quantile function 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) where 𝜆₁ is a mean from (2.17), (2.21), (2.25), or (2.30) and 𝛿 is a constant that scales 𝜆₂ in (2.18), (2.22), (2.26), or (2.31) and in the denominator of (3.5) to 1/√𝜋 as 𝛿=𝛿_{𝛾𝑗,𝜅𝑗}=2𝛾_𝑗√𝜋2𝛾_𝑗−ℎ₁𝛾_𝑗+𝜅_𝑗²+𝛾_𝑗−𝜅_𝑗²ℎ₂−ℎ₃+ℎ₅𝛾_𝑗+𝜅_𝑗²,𝛿=𝛿_𝛾𝑗=sin𝛾_𝑗𝜋𝛾_𝑗𝜋^3/2,𝛿=𝛿_𝜅𝑗=2√𝜋2+𝜅_𝑗−4𝑝₁+4𝑝₂+𝜅_𝑗𝑝₃−𝑝₄,𝛿=𝛿_{𝜅ℒ𝑗,𝜅ℛ𝑗}=4/√𝜋4−4𝑝₅𝜅_ℒ𝑗+4𝑝₆𝜅_ℒ𝑗+𝑝₉𝜅²_ℒ𝑗−𝑝₁₀𝜅²_ℒ𝑗−4𝑝₇𝜅_ℛ𝑗+4𝑝₈𝜅_ℛ𝑗+𝑝₁₁𝜅²_ℛ𝑗−𝑝₁₂𝜅²_ℛ𝑗. (3.12)

Analogously, the 𝐿-correlation of 𝑞_𝑘(𝑔(Φ(𝑍_𝑘))) toward 𝑞_𝑗(𝑔(Φ(𝑍_𝑗))) is expressed as 𝜂_𝑘𝑗=2√𝜋^+∞_−∞𝑠_𝑘𝑞_𝑘𝑔Φ𝑧_𝑘,𝛾_𝑘,𝜅_𝑘,𝜅_ℒ𝑘,𝜅_ℛ𝑘Φ𝑧_𝑗𝑓_𝑗𝑘𝑑𝑧_𝑘𝑑𝑧_𝑗. (3.13)

Note also for the special case that if 𝑞_𝑗(𝑔(Φ(𝑍_𝑗))) in (3.10) and 𝑞_𝑘(𝑔(Φ(𝑍_𝑘))) in (3.13) have the same parameters, that is, 𝛾_𝑗=𝛾_𝑘; 𝜅_𝑗=𝜅_𝑘; 𝜅_ℒ𝑗=𝜅_ℒ𝑘; 𝜅_ℛ𝑗=𝜅_ℛ𝑘, then 𝜂_𝑗𝑘=𝜂_𝑘𝑗. Source code written in Mathematica [30] that implements the computation of an IC (𝜌_𝑗𝑘) based on (3.10) is provided in Algorithm 1. The details for simulating 𝛾-𝜅, 𝛾, 𝜅, and 𝜅_ℒ-𝜅_ℛ distributions with specified values of 𝐿-skew, 𝐿-kurtosis, and 𝐿-correlations are described in the next section.

(∗Intermediate Correlation for Distributions 1 and 2. See Table 3.∗) 𝜌𝑗𝑘= 0.380206;   Needs[“MultivariateStatistics’”] 𝑓𝑗𝑘 = PDF[MultinormalDistribution[{0,0},{{1,𝜌𝑗𝑘},{𝜌𝑗𝑘,1}}],{𝑍𝑗,𝑍𝑘}]; 𝐹𝑗 = CDF[NormalDistribution[0, 1],𝑍𝑗]; 𝐹𝑘 = CDF[NormalDistribution[0, 1],𝑍𝑘];   (∗Parameters for Figure 1(a)∗) 𝛾𝑗=0.15; 𝜅𝑗=0.075;   𝑋𝑗 = Log[𝐹𝑗/(1 −𝐹𝑗)];   (∗Quantile function from (2.3)*) 𝑞𝛾𝑗,𝜅𝑗=𝛾𝑗−1∗(Exp[𝛾𝑗𝑋𝑗]−1)∗Exp[𝜅𝑗∗Abs[𝑋𝑗]];   (∗Standardizing constants 𝜆1 from (2.17) and 𝛿𝛾𝑗,𝜅𝑗 from (3.12)*) 𝑆𝑞𝛾𝑗,𝜅𝑗=𝛿𝛾𝑗,𝜅𝑗∗(𝑞𝛾𝑗,𝜅𝑗−𝜆1);   (∗Compute the specified 𝐿-correlation in Table 1 for Distributions 1 and 2∗) 𝜂𝑗𝑘=2√𝜋∗NIntegrate[𝑆𝑞𝛾𝑗,𝜅𝑗∗𝐹𝑘∗𝑓𝑗𝑘,{𝑍𝑗,−10,10},{𝑍𝑘,−10,10} Method→MultiDimensional]   0.40

Algorithm 1 

Mathematica source code for computing intermediate correlations for specified 𝐿-correlations. The example is for distribution 𝑗=1 towards distribution 𝑘=2 (𝜂₁₂) in Figures 1(a) and 2(a). See Tables 1 and 3.

4. The Procedure for Simulation and Monte Carlo Study

To implement the method for simulating 𝛾-𝜅, 𝛾, 𝜅, and 𝜅_ℒ-𝜅_ℛ distributions with specified 𝐿-moments and 𝐿-correlations we suggest the following six steps.(1)Specify the 𝐿-moments for 𝑇 transformations of the forms in (2.3)–(2.6) that is, 𝑞₁(𝑔(Φ(𝑍₁))),…,𝑞_𝑇(𝑔(Φ(𝑍_𝑇))) and obtain the solutions for the parameters of 𝛾-𝜅, 𝛾, 𝜅, or 𝜅_ℒ-𝜅_ℛ distributions by solving (2.19), (2.20); (2.23), (2.24); (2.27), (2.28); (2.32), (2.33) using the specified values of 𝐿-skew (𝜏₃) and 𝐿-kurtosis (𝜏₄) for each distribution. Specify a 𝑇×𝑇 matrix of 𝐿-correlations (𝜂_𝑗𝑘) for 𝑞_𝑗(𝑔(Φ(𝑍_𝑗))) toward 𝑞_𝑘(𝑔(Φ(𝑍_𝑘))), where 𝑗<𝑘∈{1,2,…,𝑇}.(2)Compute the (Pearson) intermediate correlations (ICs) 𝜌_𝑗𝑘 by substituting the solutions of the parameters from Step (1) into (3.10) and then numerically integrate to solve for 𝜌_𝑗𝑘 (see Algorithm 1 for an example). Repeat this step separately for all 𝑇(𝑇−1)/2 pairwise combinations of correlations.(3)Assemble the ICs into a 𝑇×𝑇 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 generate 𝑇 standard normal variables (𝑍₁,…,𝑍_𝑇) correlated at the intermediate levels as follows: 𝑍₁=𝑎₁₁𝑉₁,𝑍₂=𝑎₁₂𝑉₁+𝑎₂₂𝑉₂,⋮𝑍_𝑗=𝑎_1𝑗𝑉₁+𝑎_2𝑗𝑉₂+⋯+𝑎_𝑖𝑗𝑉_𝑖+⋯+𝑎_𝑗𝑗𝑉_𝑗,⋮𝑍_𝑇=𝑎_1𝑇𝑉₁+𝑎_2𝑇𝑉₂+⋯+𝑎_𝑖𝑇𝑉_𝑖+⋯+𝑎_𝑗𝑇𝑉_𝑗+⋯+𝑎_𝑇𝑇𝑉_𝑇, (4.1) where 𝑉₁,…,𝑉_𝑇 are independent standard normal random variables and where 𝑎_𝑖𝑗 represents the element in the 𝑖th row and the 𝑗th column of the matrix associated with the Cholesky factorization performed in Step (3).(5)Substitute 𝑍₁,…,𝑍_𝑇 from Step (4) into the following Taylor series-based expansion for the standard normal cdf [30]: Φ𝑍_𝑗=12+𝜙𝑍_𝑗𝑍_𝑗+𝑍³_𝑗3+𝑍⁵_𝑗()+𝑍3⋅5(^7𝑗),3⋅5⋅7+⋯ (4.2) where 𝜙(𝑍_𝑗) denotes the standard normal pdf and where the absolute error associated with (4.2) is less than 8×10⁻¹⁶.(6)Substitute the zero-one uniform deviates, Φ(𝑍_𝑗), generated from Step (5) into the 𝑇 equations of the form of 𝑞_𝑗(𝑔(Φ(𝑍_𝑗))), where 𝑔(Φ(𝑍_𝑗))=ln(Φ(𝑍_𝑗)/(1−Φ(𝑍_𝑗))) is standard logistic to generate the 𝛾-𝜅, 𝛾, 𝜅, and 𝜅_ℒ-𝜅_ℛ distributions with the specified 𝐿-moments and 𝐿-correlations.

To demonstrate the steps above and evaluate the proposed method, a comparison between the proposed 𝐿-moment and conventional product-moment-based procedures is subsequently described. Specifically, the heavy-tailed distributions in Figures 1(a), 2(a), 3(b), and 4(a) are used as a basis for a comparison using the specified correlation matrix in Table 1. Tables 2 and 3 give the solved IC matrices for the conventional moment and 𝐿-moment-based methods, respectively. See Algorithm 2 for an example of 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 𝑍₁,…,𝑍₄ with the specified ICs by making use of the formulae given in (4.1) of Step (4) with 𝑇=4. The values of 𝑍₁,…,𝑍₄ are subsequently transformed to Φ(𝑍₁),…,Φ(𝑍₄) using (4.2) and then substituted into equations of the forms in (2.3)–(2.6) to produce 𝑞₁(𝑔(Φ(𝑍₁))),…,𝑞₄(𝑔(Φ(𝑍₄))) for both methods.

Table 1 

Specified correlation matrix for the distributions in Figures 1(a), 2(a), 3(b), and 4(a).

Table 2 

Intermediate correlations for the conventional moment procedure.

Table 3 

Intermediate correlations for the L-moment procedure.

Table 4 

Cholesky decompositions for the conventional moment procedure.

Table 5 

Cholesky decompositions for the L-moment procedure.

(∗Intermediate Correlation for Distributions 1 and 4. See Table 2.∗) 𝜌𝑗𝑘= 0.835142;   Needs[“MultivariateStatistics’”] 𝑓𝑗𝑘=PDF[MultinormalDistribution[{0,0},{{1,𝜌𝑗𝑘},{𝜌𝑗𝑘,1}}],{𝑍𝑗,𝑍𝑘}]; 𝐹𝑗=CDF[NormalDistribution[0,1],𝑍𝑗]; 𝐹𝑘=CDF[NormalDistribution[0,1],𝑍𝑘];   (∗Parameters for Figure 1(a) and Figure 4(a) in Figure 4∗) 𝛾𝑗=0.15; 𝜅𝑗=0.075; 𝜅ℒ𝑘=0.0475; 𝜅ℛ𝑘=0.1358;   𝑋𝑗=Log[𝐹𝑗/(1 −𝐹𝑗)]; 𝑋𝑘=Log[𝐹𝑘/(1−𝐹𝑘)];   (∗Quantile functions from (2.3) and (2.6)*) 𝑞𝛾𝑗,𝜅𝑗=𝛾𝑗−1∗(Exp[𝛾𝑗𝑋𝑗]−1)∗Exp[𝜅𝑗∗Abs[𝑋𝑗]]; 𝑞ℒ𝑘=𝑋𝑘∗Exp[𝜅ℒ𝑘∗Abs[𝑋𝑘]]; 𝑞ℛ𝑘=𝑋𝑘∗Exp[𝜅ℛ𝑘∗Abs[𝑋𝑘]];   (∗Standardizing constants 𝜆1𝑗, 𝜆1𝑘 ((2.17), (2.30)) and 𝛼2𝑗, 𝛼2𝑘 ((A.2), (B.2)) from Appendices A and B∗) 𝑆𝑞𝛾𝑗,𝜅𝑗=(𝑞𝛾𝑗,𝜅𝑗−𝜆1𝑗)/𝛼2𝑗; 𝑆𝑞ℒ𝑘=(𝑞ℒ𝑘−𝜆1𝑘)/𝛼2𝑘; 𝑆𝑞ℛ𝑘=(𝑞ℛ𝑘−𝜆1𝑘)/𝛼2𝑘;   (∗Compute the specified Pearson correlation in Table 1 for Distributions 1 and 4.∗) 𝜌∗𝑗𝑘 = NIntegrate[𝑆𝑞𝛾𝑗,𝜅𝑗∗Piecewise[{{𝑆𝑞ℒ𝑘, 𝑋𝑘≤0}, {𝑆𝑞ℛ𝑘, 𝑋𝑘>0}}]∗𝑓𝑗𝑘, {𝑍𝑗, −10, 10}, {𝑍𝑘, −10, 10}, Method→MultiDimensional]   0.80

Algorithm 2 

Mathematica source code for computing intermediate correlations for specified conventional Pearson correlations. The example is for distributions 𝑗=1 and 𝑘=4 (𝜌^∗₁₄) in Figures 1(a) and 4(a). See Tables 1 and 2.

In terms of the simulation, a Fortran algorithm was written for both methods to generate 25,000 independent sample estimates for the specified parameters of (i) conventional skew (𝛼₃), kurtosis (𝛼₄), and Pearson correlation (𝜌^∗_𝑗𝑘) and (ii) 𝐿-skew (𝜏₃), 𝐿-kurtosis (𝜏₄), and 𝐿-correlation (𝜂_𝑗𝑘). All estimates were based on sample sizes of 𝑛=25 and 𝑛=1000. The formulae used for computing estimates of 𝛼_3,4 were based on Fisher’s 𝑘-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 of 𝜏_3,4 were Headrick’s Equations (2.4) and (2.6) [25]. The estimate for 𝜌^∗_𝑗𝑘 was based on the usual formula for the Pearson product moment of correlation statistic, and the estimate for 𝜂_𝑗𝑘 was computed based on (3.5) using the empirical forms of the cdfs in (3.1) and (3.3). The estimates for 𝜌^∗_𝑗𝑘 and 𝜂_𝑗𝑘 were both transformed using Fisher’s 𝑧^′ transformation. Bias-corrected accelerated bootstrapped average (mean) estimates, confidence intervals (CIs), and standard errors were subsequently obtained for the estimates associated with the parameters (𝛼_3,4, 𝜏_3,4, 𝑧^′_𝜌^∗_𝑗𝑘, 𝑧^′_𝜂𝑗𝑘) using 10,000 resamples via the commercial software package Spotfire S+ [31]. The bootstrap results for the estimates of the means and CIs associated with 𝑧^′_𝜌^∗_𝑗𝑘 and 𝑧^′_𝜂𝑗𝑘 were transformed back to their original metrics (i.e., estimates for 𝜌^∗_𝑗𝑘 and 𝜂_𝑗𝑘). Further, if a parameter (𝑃) was outside its associated bootstrap CI, then an index of relative bias (RB) was computed for the estimate (𝐸) as RB=((𝐸−𝑃)/𝑃)×100. Note that the small amount of bias associated with any bootstrap CI containing a parameter was considered negligible and thus not reported. The results of the simulation are reported in Tables 6–13 and are discussed in the next section.

Table 6 

Skew (𝛼₃) and Kurtosis (𝛼₄) results for the conventional moment procedure. Sample size of 𝑛=25.

Table 7 

Skew (𝛼₃) and Kurtosis (𝛼₄) results for the conventional moment procedure. Sample size of 𝑛=1000.

Table 8 

𝐿-skew (𝜏₃) and 𝐿-kurtosis (𝜏₄) results. Sample size of 𝑛=25.

Table 9 

𝐿-skew (𝜏₃) and 𝐿-kurtosis (𝜏₄) results. Sample size of 𝑛=1000.

Table 10 

Correlation results for the conventional moment procedure. Sample size of 𝑛=25.

Table 11 

Correlation results for the conventional moment procedure. Sample size of 𝑛=1000.

Table 12 

Correlation results for the L-moment procedure. Sample size of 𝑛=25.

Table 13 

Correlation results for the L-moment procedure. Sample size of 𝑛=1000.

5. Discussion and Conclusion

One of the primary advantages that 𝐿-moments have over conventional moment-based estimators is that they can be far less biased when sampling is from distributions with more severe departures from normality (e.g., [21, 25]). Inspection of the simulation results in Tables 6 and 7 of this study clearly indicates that this is also the case for the system of 𝛾-𝜅, 𝛾, 𝜅, and 𝜅_ℒ-𝜅_ℛ distributions. Specifically, the superiority that estimates of 𝐿-moment ratios (𝜏₃,𝜏₄) have over their corresponding conventional moment-based counterparts (𝛼₃,𝛼₄) is obvious. For example, with samples of size 𝑛=25 the estimates of skew and kurtosis for Distribution 1 (Tables 6 and 8) were, on average, only 26.80% and 2.38% of their associated population parameters whereas the estimates of 𝐿-skew and 𝐿-kurtosis were 87.80% and 94.81% of their respective parameters. It is also evident from Tables 6 and 8 that 𝐿-skew and 𝐿-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, inspection of Tables 7 and 9 indicates RSE measures of RSE(𝛼₃)=0.436 and RSE(𝛼₄)=1.12 compared with RSE(𝜏₃)=0.109 and RSE(𝜏₄)=0.048. This demonstrates that 𝐿-skew and 𝐿-kurtosis have more precision because they have less variance around their estimates.

The results associated with the conventional Pearson and 𝐿-correlations are presented in Tables 10–13. Overall inspection of these tables indicates that the 𝐿-correlation is substantially superior to the Pearson correlation in terms of relative bias. For example, in terms of a moderate correlation (Table 10, 𝑛=25, 𝜌^∗₁₂=0.40) the relative bias for Distributions 1 and 2 was 9.73% for the Pearson correlation compared to only 1.85% for the 𝐿-correlation (Table 12, 𝑛=25, 𝜂₁₂=0.40). Further, for large sample sizes (Table 13, 𝑛=1000), the 𝐿-correlation bootstrap CIs contained all of the population parameters whereas the Pearson correlation CIs contained none of the parameters (Table 11).

In summary, the proposed 𝐿-moment-based system of 𝛾-𝜅, 𝛾, 𝜅, and 𝜅_ℒ-𝜅_ℛ distributions is an attractive alternative to the traditional conventional-moment-based system. In particular, the 𝐿-moment-based system has distinct advantages when heavy-tailed distributions are of concern. Finally, we would note that Mathematica Version 8.0 [30] source code is available from the authors for implementing the 𝐿-moment-based method.

Appendices

A. System of Conventional Moment-Based Equations for 𝛾-𝜅Distributions

The moments (𝜇_𝑟=1,…,4) associated with the class of 𝛾-𝜅 distributions in (2.3) can be determined from 𝜇_𝑟=^+∞_−∞𝑞𝑥,𝛾,𝜅^𝑟𝑓_𝑋(𝑥)𝑑𝑥. (A.1)

In general, to obtain defined values of 𝜇_𝑟 based on (A.1) we must have (i) 𝛾+𝜅<1/𝑟, (ii) 𝜅<1/𝑟, and (iii) 1/𝑟+𝛾<𝜅.

The mean, variance, skew, and kurtosis are defined in general as in [17] 𝛼₁=𝜇₁,𝛼²₂=𝜇₂−𝜇²₁,𝛼₃=𝜇₃−3𝜇₂𝜇₁+2𝜇³₁𝛼³₂,𝛼₄=𝜇₄−4𝜇₃𝜇₁−3𝜇²₂+12𝜇₂𝜇²₁−6𝜇⁴₁𝛼⁴₂. (A.2)

The moments associated with the location and scale parameters in (A.2) are 𝛼₁=𝜇₁=−𝛾−𝜅ℎ₁₁+𝛾−𝜅ℎ₁₂+−𝛾+𝜅ℎ₁₃+2𝜅h₁₄+𝛾+𝜅ℎ₁₅−2𝜅ℎ₁₆2𝛾 (A.3) with harmonic number functions 𝐻[⋅] of ℎ₁₁=𝐻[(1/2)(−1−𝛾−𝜅)], ℎ₁₂=𝐻[(1/2)(−1+𝛾−𝜅)], ℎ₁₃=𝐻[(1/2)(𝛾−𝜅)], ℎ₁₄=𝐻[(1/2)(−1−𝜅)], ℎ₁₅=𝐻[(−1/2)(𝛾+𝜅)], ℎ₁₆=𝐻[(−1/2)𝜅], and 𝜇₂=−𝛾−2𝜅ℎ₂₁−2𝜅ℎ₂₂−𝛾+𝜅ℎ₂₃+𝛾+𝜅ℎ₂₄+𝛾+2𝜅ℎ₂₅−𝛾+2𝜅ℎ₂₆+𝛾−2𝜅ℎ₂₇+𝛾−𝜅ℎ₂₈+−𝛾+𝜅ℎ₂₉+2𝜅ℎ₂₁₀/2𝛾² (A.4) with harmonic number functions 𝐻[⋅] of ℎ₂₁=𝐻[(1/2)(−1+𝛾−2𝜅)], ℎ₂₂=𝐻[−1/2−𝜅], ℎ₂₃=𝐻[−1/2−𝛾−𝜅], ℎ₂₄=𝐻[−𝛾−𝜅], ℎ₂₅=𝐻[−1/2−𝛾/2−𝜅], ℎ₂₆=𝐻[−𝛾/2−𝜅], ℎ₂₇=𝐻[𝛾/2−𝜅], ℎ₂₈=𝐻[−1/2+𝛾−𝜅], ℎ₂₉=𝐻[𝛾−𝜅], and ℎ₂₁₀=𝐻[−𝜅].

The moments associated with the shape parameters of skew and kurtosis in (A.2) are 𝜇₃=3−𝛾+𝜅ℎ₃₁−𝛾+3𝜅ℎ₃₂+𝛾−3𝜅ℎ₃₃−𝛾−3𝜅ℎ₃₄+𝛾−𝜅ℎ₃₅+2𝜅ℎ₃₆+2𝛾+3𝜅ℎ₃₇+−2𝛾−3𝜅ℎ₃₈+𝛾+3𝜅ℎ₃₉+−2𝛾+3𝜅ℎ₃₁₀+2𝛾−3𝜅ℎ₃₁₁+−𝛾+𝜅ℎ₃₁₂−2𝜅ℎ₃₁₃+𝛾+𝜅ℎ₃₁₄/2𝛾³ (A.5) with harmonic number functions 𝐻[⋅] of ℎ₃₁=𝐻[(1/2)(−1−3𝛾−3𝜅)], ℎ₃₂=𝐻[(1/2)(−1−𝛾−3𝜅)], ℎ₃₃=𝐻[(1/2)(−1+𝛾−3𝜅)], ℎ₃₄=𝐻[(1/2)(𝛾−3𝜅)], ℎ₃₅=𝐻[(1/2)(−1+3𝛾−3𝜅)], ℎ₃₆=𝐻[−1/2−3𝜅/2], ℎ₃₇=𝐻[−1/2−𝛾−3𝜅/2], ℎ₃₈=𝐻[−𝛾−3𝜅/2], ℎ₃₉=𝐻[−𝛾/2−3𝜅/2], ℎ₃₁₀=𝐻[−1/2+𝛾−3𝜅/2], ℎ₃₁₁=𝐻[𝛾−3𝜅/2], ℎ₃₁₂=𝐻[3(𝛾−𝜅)/2], ℎ₃₁₃=𝐻[−3𝜅/2], ℎ₃₁₄=𝐻[−3(𝛾+𝜅)/2], and 𝜇₄=2−𝛾−4𝜅−2𝜅ℎ₄₁−𝛾+𝜅ℎ₄₂+3𝛾+4𝜅ℎ₄₃+−3𝛾−4𝜅ℎ₄₄−3𝛾+2𝜅ℎ₄₅+3𝛾+2𝜅ℎ₄₆+𝛾+4𝜅ℎ₄₇−𝛾+4𝜅ℎ₄₈+𝛾−4𝜅ℎ₄₉+3𝛾−2𝜅ℎ₄₁₀+−3𝛾−2𝜅ℎ₄₁₁+−3𝛾+4𝜅ℎ₄₁₂+3𝛾−4𝜅ℎ₄₁₃+𝛾−𝜅ℎ₄₁₄+−𝛾+𝜅ℎ₄₁₅+2𝜅ℎ₄₁₆+𝛾+𝜅ℎ₄₁₇/𝛾⁴ (A.6) with harmonic number functions 𝐻[⋅] of ℎ₄₁=𝐻[−1/2−2𝜅], ℎ₄₂=𝐻[−1/2−2𝛾−2𝜅], ℎ₄₃=𝐻[−1/2−3𝛾/2−2𝜅], ℎ₄₄=𝐻[−3𝛾/2−2𝜅], ℎ₄₅=𝐻[−1/2−𝛾−2𝜅], ℎ₄₆=𝐻[−𝛾−2𝜅], ℎ₄₇=𝐻[−1/2−𝛾/2−2𝜅], ℎ₄₈=𝐻[−𝛾/2−2𝜅], ℎ₄₉=𝐻[𝛾/2−2𝜅], ℎ₄₁₀=𝐻[−1/2+𝛾−2𝜅], ℎ₄₁₁=𝐻[𝛾−2𝜅], ℎ₄₁₂=𝐻[−1/2+3𝛾/2−2𝜅], ℎ₄₁₃=𝐻[3𝛾/2−2𝜅], ℎ₄₁₄=𝐻[−1/2+2𝛾−2𝜅], ℎ₄₁₅=𝐻[2𝛾−2𝜅], ℎ₄₁₆=𝐻[−2𝜅], and ℎ₄₁₇=𝐻[−2𝛾−2𝜅].

The moments for the asymmetric 𝛾 class of distributions in (2.4) can be obtained by taking the limit of (A.3)–(A.6) as 𝜅→0. Analogously, the moments for the symmetric 𝜅 class of distributions in (2.5) can be obtained by taking the limit of (A.3)–(A.6) as 𝛾→0.

B. System of Conventional Moment-Based Equations for 𝜅_ℒ-𝜅_ℛ Distributions

The moments (𝜇_𝑟=1,…,4) associated with the asymmetric 𝜅_ℒ-𝜅_ℛ class of distributions in (2.6) can be determined from 𝜇_𝑟=⁰_−∞𝑞𝑥,𝜅_ℒ^𝑟𝑓_𝑋(𝑥)𝑑𝑥+₀^+∞𝑞𝑥,𝜅_ℛ^𝑟𝑓_𝑋(𝑥)𝑑𝑥. (B.1)

Using the expressions in (A.2), the moments for the location and scale parameters are 𝛼₁=𝜇₁=142𝑝₁₁−2𝑝₁₂−2𝑝₁₃+2𝑝₁₄−𝜅_ℒ𝑝₁₅+𝜅_ℒ𝑝₁₆+𝜅_ℛ𝑝₁₇−𝜅_ℛ𝑝₁₈,𝜇₂=12𝜅_ℒ1𝐻𝑍3,2−𝜅_ℒ−𝜅_ℒ𝐻𝑍3,1−𝜅_ℒ+𝑝₁₅+𝑝₁₆+12𝜅_ℛ1𝐻𝑍3,2−𝜅_ℛ−𝜅_ℛ𝐻𝑍3,1−𝜅_ℛ+𝑝₁₇−𝑝₁₉, (B.2)(B.3) where the polygamma functions 𝑃[⋅] are: 𝑝₁₁=𝑃[0,1/2−𝜅_ℒ/2], 𝑝₁₂=𝑃[0,1−𝜅_ℒ/2], 𝑝₁₃=𝑃[0,1/2−𝜅_ℛ/2], 𝑝₁₄=𝑃[0,1−𝜅_ℛ/2], 𝑝₁₅=𝑃[1,1/2−𝜅_ℒ/2], 𝑝₁₆=𝑃[1,1−𝜅_ℒ/2], 𝑝₁₇=𝑃[1,1/2−𝜅_ℛ/2], 𝑝₁₈=𝑃[1,1−𝜅_ℛ/2], and 𝑝₁₉=𝑃[1,1−𝜅_ℛ]. The notation HZ[⋅] in (B.3) is the Hurwitz zeta function. The moments related to skew and kurtosis are as follows: 𝜇₃=3162𝑝₃₁−2𝑝₃₂−𝜅_ℒ𝑝₃₃+𝜅_ℒ𝑝₃₄+38𝑝₃₅1+2𝑍3,21−3𝜅_ℛ+3𝜅_ℛ𝑍14,21−3𝜅_ℛ−3𝜅_ℛ𝑍4,1−3𝜅_ℛ2,𝜇₄=18𝑝₄₁−𝑝₄₂−𝜅_ℒ𝑝₄₃+𝜅_ℒ𝑝₄₄+182𝑝₄₅−2𝑝₄₆−𝜅_ℛ𝑝₄₇+𝜅_ℛ𝑝₄₈, (B.4)(B.5) where the polygamma functions 𝑃[⋅] are 𝑝₃₁=𝑃[2,1/2−3𝜅_ℒ/2], 𝑝₃₂=𝑃[2,1−3𝜅_ℒ/2], 𝑝₃₃=𝑃[3,1/2−3𝜅_ℒ/2], 𝑝₃₄=𝑃[3,1−3𝜅_ℒ/2], 𝑝₃₅=𝑃[2,1−3𝜅_ℛ/2], 𝑝₄₁=𝑃[3,1/2−2𝜅_ℒ], 𝑝₄₂=𝑃[3,1−2𝜅_ℒ], 𝑝₄₃=𝑃[4,1/2−2𝜅_ℒ], 𝑝₄₄=𝑃[4,1−2𝜅_ℒ], 𝑝₄₅=𝑃[3,1/2−2𝜅_ℛ], 𝑝₄₆=𝑃[3,1−2𝜅_ℛ], 𝑝₄₇=𝑃[4,1/2−2𝜅_ℛ], and 𝑝₄₈=𝑃[4,1−2𝜅_ℛ]. The notation 𝑍[⋅] in (B.4) is the zeta function.