|
Article Excerpt
1. INTRODUCTION
Real data often show significant departures from normality assumptions. It is widely acknowledged that heavy-tailed distributions are frequently encountered in empirical studies (Hampel, Ronchetti, Rousseeuw, and Stahel 1986, sec. 1.2), as are asymmetric distributions (Hill and Dixon 1982; Azzalini 1986; Hampel et al. 1986, p. 23; Azzalini and Della Valle 1996; Azzalini and Capitanio 1999). For cases where the assumption of normality is not tenable, more flexible models can be adopted to accommodate skewness and heavy tails. An attractive feature of using flexible models is that they permit likelihood inference while coping with asymmetry and kurtosis. Flexible models that include the normal distribution as a special case are especially important, because they allow continuous variation from normality to nonnormality. Thus such models can adapt to distributions that are in a neighborhood of the normal model, while offering the possibility of checking the assumption of normality through a formal test on a few model parameters.
Using a flexible model to handle nonnormality has certain benefits. When the deviation from normality involves skewness, a model is needed to establish meaningful location and scale parameters. Even when skewness is absent, a model is necessary to deduce an interpretable scale parameter in the presence of leptokurtosis. Furthermore, a model permits efficient testing of the normality assumption, and can be used for constructing prediction intervals for future observations.
Subbotin (1923) introduced the exponential power (EP) distribution with density
[f.sub.EP](x; [mu], [sigma], [alpha]) = ([sigma]c)[.sup.-1] exp(-|z|[.sup.[alpha]]/[alpha]), -[infinity] < x < [infinity],
where -[infinity] 0, [alpha] > 1, z = (x - [mu])/[sigma], and c = 2[[alpha].sup.1/[alpha]-1][GAMMA](1/[alpha]). The normal model is obtained when [alpha] = 2, and symmetric, heavy-tailed distributions are produced as the parameter [alpha] varies in the interval (1, 2). This model has been studied extensively by Box (1953), Turner (1960), Vianelli (1963), and others.
Another model, introduced by Azzalini (1985) to accommodate asymmetry, is the skew normal (SN) distribution, whose density is given by
[f.sub.SN](x; [mu], [sigma], [lambda]) = (2/[sigma])[PHI]([lambda]z)[phi](z), -[infinity] < x < [infinity],
where -[infinity] 0, z = (x - [mu])/[sigma], and [phi](z) and [PHI](z) are the normal density and distribution functions. The skewness is regulated by the parameter [lambda]; [lambda] = yields the normal model. A multivariate extension of this model has been developed by Azzalini and Della Valle (1996) and Azzalini and Capitanio (1999).
To handle both skewness and heavy tails simultaneously, Azzalini (1986) proposed the skew exponential power (SEP) distribution, which has density
[f.sub.SEP](x; [mu], [sigma], [lambda], [alpha]) = (2/[sigma])[PHI](w)[f.sub.EP](x; [mu], [sigma], [alpha]), -[infinity] < x < [infinity],
where -[infinity] 0, [alpha] > 1, w = sign(z)|z|[.sup.[alpha]/2] X [lambda](2/[alpha])[.sup.1/2], and z = (x - [mu])/[sigma]. The SEP distribution reduces to the EP when [lambda] = 0, to the SN when [alpha] = 2, and to the normal when ([lambda], [alpha]) = (0, 2). The SEP distribution is flexible, incorporating a wide range of models in a neighborhood of the normal distribution. In the SEP distribution, [mu] and [sigma] are location and scale parameters, and they do not in general denote the mean and standard deviation.
Azzalini (1986) investigated the properties of the SEP distribution. In particular, he gave formulas for the moments. The even moments of Z ~ SEP(0, 1, [lambda], [alpha]) can be expressed conveniently in terms of the gamma function,
E([Z.sup.2m]) = [[alpha].sup.2m/[alpha]][GAMMA]((2m + 1)/[alpha])/[GAMMA](1/[alpha]).
[FIGURE 1 OMITTED]
Expressions for the odd moments involve infinite series expansions (Azzalini, personal communication),
E([Z.sup.2m+1]) = [[2[[alpha].sup.(2m+1)/[alpha]][lambda]]/[[square root of [pi]][GAMMA](1/[alpha])(1 + [[lambda].sup.2])[.sup.s+1/2]]] X [[infinity].summation over (n=0)][[[GAMMA](s + n + 1/2)]/[(2n + 1)!!]]([2[[lambda].sup.2]]/[1 + [[lambda].sup.2]])[.sup.n], (1)
where s = 2(m + 1)/[alpha] and (2n + 1)!! = 1 * 3 ... (2n - 1) * (2n + 1). In addition, for an extensive grid of ([lambda], [alpha]) values, Azzalini (1986) plotted the corresponding values of ([[gamma].sub.1], [[gamma].sub.2]), where [[gamma].sub.1] is the standardized skewness and [[gamma].sub.2] is the excess in kurtosis. He found the ([[gamma].sub.1], [[gamma].sub.2]) values to be wide ranging; the [[gamma].sub.1] values varied between -2 and 2, whereas the [[gamma].sub.2] values extended between -.4 and 6. Figure 1 shows the SEP densities for [lambda] = .5, 1.5, 3, 5 and [alpha] = 1.6; Figure 2 shows the densities for [lambda] = 1.5 and [alpha] = 1.3, 1.6, 1.9. Only positive values of [lambda] are considered in these plots; when the sign of [lambda] is reversed, the density is reflected about the origin.
Inferential aspects of the EP distribution have been studied by many authors (see Chiodi 2000 for a review). In particular, Agro (1995) proved consistency, asymptotic normality, and efficiency of the maximum likelihood estimators (MLEs). The information matrix of the MLEs in the SN distribution has been derived by Azzalini (1985); other estimation issues of the SN distribution have been considered by Chiogna (1997) and by Azzalini and Capitanio (1999). In contrast, properties of likelihood inference for the SEP distribution still need investigation.
[FIGURE 2 OMITTED]
This article concerns likelihood inference for the SEP distribution; in particular, maximum likelihood estimation of [theta] = ([mu], [sigma], [lambda], [alpha]) is considered. The information matrix is derived in Section 2. Special problems that arise in the normal case ([lambda], [alpha]) = (0, 2) are discussed in Section 3. The finite-sample performance of MLEs and likelihood ratio test statistics are studied numerically in Section 4. Examples in which the SEP distribution is used for robust estimation are given in Section 5. Technical details are given in the Appendixes.
2. MAXIMUM LIKELIHOOD ESTIMATION
Let [X.sub.1],..., [X.sub.n] be a random sample drawn from the SEP distribution. The log-likelihood function for [theta] is [summation]l([theta], [X.sub.i]), where l([theta], X) is the log-likelihood for [theta] based on a single observation X, that is,
l([theta], X) = -(1/[alpha] - 1)ln[alpha] - ln[GAMMA](1/[alpha]) - ln[sigma] + ln[PHI](W) - |Z|[.sup.[alpha]]/[alpha],
Z = (X - [mu])/[sigma], and W = sign(Z)|Z|[.sup.[alpha]/2][lambda](2/[alpha])[.sup.1/2]. The score function is [summation][S.sub.[theta]]([theta], [X.sub.i]), where [S.sub.[theta]]([theta], X) = [partial derivative]l([theta], X)/[partial derivative][theta] = ([S.sub.[mu]], [S.sub.[sigma]], [S.sub.[lambda]], [S.sub.[alpha]]) is given by
[S.sub.[mu]] = -[1/[sigma]][[[phi](W)]/[[PHI](W)]][[[partial derivative]W]/[[partial derivative]Z]] + [1/[sigma]]sign(Z)|Z|[.sup.[alpha]-1],
[S.sub.[sigma]] = -[1/[sigma]][[[phi](W)]/[[PHI](W)]][[[partial derivative]W]/[[partial derivative]Z]]Z + [1/[sigma]](|Z|[.sup.[alpha]] - 1).
[S.sub.[lambda]] = [[[phi](W)]/[[PHI](W)]][[[partial derivative]W]/[[partial derivative][lambda]]],
and
[S.sub.[alpha]] = [[[phi](W)]/[[PHI](W)]][[[partial derivative]W]/[[partial derivative][alpha]]] - [[ln[alpha] + ([alpha] - 1) + [PSI](1/[alpha]) - [alpha]|Z|[.sup.[alpha]]ln|Z| + |Z|[.sup.[alpha]]]/[[alpha].sup.2]],
with [partial derivative]W/[partial derivative]Z = |Z|[.sup.[alpha]/2-1][lambda](2/[alpha])[.sup.-1/2], [partial derivative]W/[partial derivative][lambda] = W/[lambda], [partial derivative]W/[partial derivative][alpha] = W(ln|Z| - 1/[alpha])/2, and [PSI](r) = [partial derivative]ln[GAMMA](r)/[partial derivative]r.
Let [^.[theta]] = ([^.[mu]], [^.[sigma]], [^.[lambda]], [^.[alpha]]) be the MLE. It is easily verified that [^.[sigma]] = ([[summation].sub.i=1.sup.n] |[X.sub.i] - [^.[mu]]|[.sup.[^.[alpha]]]/n)[.sup.1/[^.[alpha]]].
The information matrix involves the quantities
[[eta].sub.k] = E[{[[phi](W)]/[[PHI](W)]}[.sup.2]|Z|[.sup.k]],
[[xi].sub.k] = E[{[[phi](W)]/[[PHI](W)]}[.sup.2]|Z|[.sup.k]ln|Z|],
[[tau].sub.k] = E[{[[phi](W)]/[[PHI](W)]}[.sup.2]sign(Z)|Z|[.sup.k]ln|Z|],
[[upsilon].sub.k] = E[{[[phi](W)]/[[PHI](W)]}[.sup.2]|Z|[.sup.k](ln|Z|)[.sup.2]],
[[zeta].sub.k] = E[{[[phi](W)]/[[PHI](W)]}[.sup.2]sign(Z)|Z|[.sup.k]],
[[rho].sub.k] = E{sign(Z)|Z|[.sup.k]ln|Z|},
and [delta] = [[lambda].sup.2]/(1 + [[lambda].sup.2]). The elements of the expected information matrix are shown in Appendix A to be
[I.sub.[mu][mu]] = [[[alpha][[lambda].sup.2][[eta].sub.[alpha]-2]]/[2[[sigma].sup.2]]] + [[[[alpha].sup.1-2/[alpha]]([alpha] - 1)]/[[sigma].sup.2]][[[GAMMA](1 - 1/[alpha])]/[[GAMMA](1/[alpha])]],
[I.sub.[mu][sigma]] = [[[[alpha].sup.2-1/[alpha]][[delta].sup.1/2](1 + [delta])]/[2[[sigma].sup.2][GAMMA](1/[alpha])]] + [[[alpha][[lambda].sup.2][[zeta].sub.[alpha]-1]]/[2[[sigma].sup.2]]],
[I.sub.[mu][lambda]] = [[[[alpha].sup.1-1/[alpha]][[delta].sup.1/2](1...
|
