The Asymmetric Alpha-Power Skew-t Distribution

In this paper, we propose a new asymmetric and heavy-tail model that generalizes both the skew-t and power-t models. Properties of the model are studied in detail. The score functions and the elements of the observed information matrix are given. The process to estimate the parameters in model is discussed by using the maximum likelihood approach. Also, the observed information matrix is shown to be non-singular at the whole parametric space. Two applications to real data sets are reported to demonstrate the usefulness of this new model.


Introduction
In recent years, there has been considerable interest in the statistical literature related to flexible families of distributions able of modeling data that present high degree of asymmetry, with kurtosis index greater or smaller than the captured by normal model. In this context, two proposals that have shown a promising behavior in this type of situations are the skew-normal (SN) distribution of Azzalini [1] and the power-normal (PN) distribution of Durrans [2]. The SN distribution has been widely studied by many authors, and its main drawback is that it presents singular Fisher information matrix, implying the inference is useless from the theory of large samples using the maximum likelihood (ML) approach. Although the PN model has a shorter asymmetry range than SN distribution, it presents non-singular information matrix and can easily be extended to censored scenarios, as it has a simple distribution function, see, for example, in Martínez-Flórez et al. [3].
The PN model is part of a wide family of distributions known as alpha-power, which has been widely studied by many authors. In addition to the normal distribution, the Birnbaum-Saunders (BS) distribution [4] has also been considered, see, for example, in Martínez-Flórez et al. [5], who propose an extension of the BS distribution based on the asymmetric alpha-power family of distributions to illustrate the applicability of the new proposal with a data set is related to the lifetimes in cycles ×10 −3 n = 101 aluminum 6061 − T6 pieces cut in parallel angle to the rotation direction of rolling at the rate of 18 cycles per second and maximum stress of 21.000 psi. More details of the PN distribution can be found in Gupta and Gupta [6] and Pewsey et al. [7].
An alternative propose for modeling asymmetric data that unifies the two previous approaches was introduced by Martínez-Flórez et al. [8]. The proposed model, which is called alpha-power skew-normal (APSN), has non-singular Fisher information matrix, and it can fit data with much more asymmetry than PN models it can handle. In addition, symmetry can be tested by using the likelihood ratio statistic, as the properties of large samples are satisfied for the ML estimator.
Another set of distributions with non-singular information matrices, useful for modeling asymmetric and heavy-tailed data, are based on generalizations of the Student-t distribution, see, for example, in [9][10][11][12][13]. Azzalini and Capitanio [9] for example, introduced a skew-t (ST) distribution as an extension of the SN model for modeling asymmetric and heavy-tailed data as follows; The random variable X is said to have the ST distribution with parameter λ and degrees of freedom ν, if X has the probability density function (PDF) given by where λ ∈ R is a parameter that controls the skewness of the distribution, and f T (·; ν) and F T (·; ν) denote the PDF and the cumulative distribution function (CDF) of a standard Student-t distribution with ν degree of freedom, respectively. The ST distribution, like an extension of the SN model, inherits the problem of the singularity of the information matrix and before this inconvenience Zhao and Kim [14] proposed the power Student-t (PT) distribution, whose information matrix is non-singular and for a given degree of freedom, the kurtosis range surpasses the kurtosis range of the skew-t model at all times.
The PT distribution is defined as follows. The random variable X is said to have the PT distribution with parameter α, and degrees of freedom ν, if X has PDF given by where α > 0 is a parameter that controls the form of the distribution, and, again, f T (·; ν) and F T (·; ν) denote the PDF and the CDF of a standard Student-t distribution, respectively. Based on the properties of the ST model, to fit data with high degree of asymmetry and the characteristic of the PN model to capture kurtosis larger than the normal model, in this paper, we introduce a new distribution for modeling asymmetric and heavy-tailed data. The proposed model possess non-singular information matrix, and it is able to fit data with far more asymmetry than ST and PT models can handle and with large sample properties satisfied for the ML estimator. The model introduced in this paper is named as alpha-power skew-t (APST) model and it extends both, ST and PT models. The APSN model by Martínez-Flórez et al. [8] is also a particular case when ν tends to infinite. Note that symmetry can be tested using the likelihood ratio statistics with its large sample chi-square distribution.
The rest of this paper is organized as follows. Section 2 introduces the APST model and some of its properties like moments are studied. In particular, skewness and kurtosis indices are computed showing that their ranges surpass those of the ST and PT models. Section 3 deals with the ML estimation for the location-scale situation and its observed information matrix is derived. The extension to censored data is also presented. Finally, two applications are shown in Section 4, revealing that the model proposed can present much improvement over competitors.

The Alpha-Power Skew-t Distribution
Definition 1. The random variable X is said to have an alpha-power skew-t (APST) distribution, if X has PDF given by for x ∈ R, λ ∈ R, and α, ν ∈ R + . Functions f ST (·) and F ST (·) denote the PDF and the CDF of the standard ST distribution. A random variable having f APST (x; λ, α, ν) distribution is denoted shortly by X ∼ APST(λ, α, ν). Figure 1 displays the form of the APST distribution for some selected values of the parameters λ and α for ν = 6. Note from the figure that the asymmetry and kurtosis of the APST distribution are affected by the parameters α and λ; therefore, the APST model is more flexible to model data that can be highly skewed, as well as heavier tails than ST and PT models.
The following result provides some special cases of the model (3), which occur for different values of λ, α, and ν.
Proof. The proof of (i)-(vii) is immediate from the definition of APST distribution.

Moments
The following proposition presents an expression to compute the k-th moment of a random variable APST(λ, α, ν).

Proof. We have by definition that
which is the expected value of the function F −1 where Y follows a beta distribution with parameters α and 1.
The indices of skewness ( β 1 ) and kurtosis (β 2 ) of APST distribution can be calculated by using the moments (4) as follows, Table 1 presents the ranges of possible values for the indices of asymmetry and kurtosis for ST(λ, ν), PT(α, ν), and APST(λ, α, ν) distributions, for values of λ between −40 and 40, values of α between 0.5 and 50, and for values of ν = 2, 3, 4, 5, 6, 7. It can seen from Table 1 that the length of the admissible intervals for the skewness and the kurtosis parameters of the APST distribution are larger than the corresponding intervals of the ST and PT distributions. This is an indicator that the APST model is more flexible in terms of asymmetry and kurtosis than the ST and PT models. Table 1. Skewness and kurtosis for the models ST(λ, ν), PT(α, ν), and APST(λ, α, ν), for λ ∈ (−40, 40), α ∈ (0.5, 50) and ν = 2, . . . 7.

Distribution Function
Proof. The proof is immediate and it follows from results of Durrans [2].
The inversion method can be used to generate a random variable with APST distribution. Thus, taking λ ∈ R, α, ν ∈ R + and a random variable with uniform distribution, namely, U ∼ U(0, 1), random variable X with APST(λ, α, ν) distribution is generated by taking

Remark 1.
We consider a truncated APST(λ, α) distribution to obtain a new and useful lifetime distribution. A random variable T has a truncated alpha-power skew-t distribution (at zero), denoted by TAPST(λ, α, ν), if its PDF is given by The survival and hazard rate functions of a random variable T following a TAPST(λ, α, ν) distribution are given by and respectively.

Location and Scale Extension
We can also consider a generalization of a APST distribution by adding location and scale parameters. The following definition gives a generalization of the APST model.
The APST density of location and scale is defined as the distribution of Y = µ + σX, for µ ∈ R and σ > 0. The corresponding PDF is given by for λ ∈ R and α, ν ∈ R + . A random variable following a APST distribution of location and scale is denoted by The k-th moment of a random variable Y ∼ APST(µ, σ, λ, α, ν) can be obtained from the formula

Statistical Inference for APST Distribution
This section concerns likelihood inference about the parameter vector θ = (µ, σ, λ, α, ν) of the location-scale family defined in Equation (9). Let Y = (Y 1 , . . . , Y n ) be a random sample of the distribution APST(µ, σ, λ, α, ν). The log-likelihood function for θ = (µ, σ, λ, α, ν) can be written as follows, where z i = (y i − µ)/σ. Thus, by differentiating the log-likelihood function, we obtain the following score equations, i +ν for i = 1, . . . , n, and g(x; ν) is the function defined by Equations (11)-(15) include nonlinear functions; therefore, it is not possible to obtain explicit forms of the maximum likelihood estimators (MLEs), and they must be calculated by using numerical methods. In this work, we used the maxLik function of R Development Core Team [15] which uses the Newton-Raphson optimization method. The elements of the observed information matrix are easily obtained after calculating the second derivative of the log-likelihood function and multiplying by −1, that is, where θ = (µ, σ, λ, α, ν) . This elements are given in the Appendix A. To find the standard errors (EE) of the MLEs and calculate confidence intervals, the information matrix I (or Fisher information) must be calculated, which is defined as the expected value of the second derived from the log-likelihood function or less the expected value of the Hessian matrix; from this matrix, we calculate the EE as the diagonal elements of the inverse of this matrix. The elements of the I matrix are obtained as The role of the Fisher information in the asymptotic theory of maximum-likelihood estimation was emphasized by Ronald Fisher following some initial results by Francis Edgeworth, see Lehman and Casella [16] and Frieden [17] for more details. The Fisher-information matrix is used to calculate the covariance matrices associated with maximum-likelihood estimates, and it can also be used in the formulation of test statistics, such as the Wald test.
As the expected value under the APST distribution and the second-order derivatives are not direct, numerical methods must be used to obtain the explicit form of the information matrix I. Therefore, we use the observed information matrix to calculate the standard errors in the rest of the document.
When ν tends to infinite the ST distribution converges to the SN distribution and we recall that the information matrix of a random variable X ∼ SN(µ, σ, λ) which is denoted by I λ (ϕ), where ϕ = (µ, σ, λ) , is singular for λ = 0. Therefore, it is convenient to use a centered parameterization of the ST distribution proposed by Arellano-Valle and Azzalini [18].
The centered parameterization of the SN distribution was proposed as an alternative to the problem of singularity of the information matrix of the SN when λ = 0. Arellano-Valle and Azzalini [19] proposed a second representation of the SN by defining a new random variable X as where µ ∈ R and σ > 0 are parameters of the random variable X and Z ∼ SN(λ). This representation is called centered parameterization, as E[X] = µ and Var[X] = σ 2 and it is denoted by CSN(µ, σ, γ 1 ), Under the centered parameterization model, µ, σ, and γ 1 = β 1 represent the mean, the standard deviation and the skewness index of X, respectively. If Z ∼ SN(λ) then E[Z] = bδ and Var[Z] = 1 − (bδ) 2 , where b = √ 2/π and δ = λ/ √ 1 + λ 2 ; it has that the random variable X can be written as X = µ + σZ which has SN(λ 1 , λ 2 , λ) distribution, where with c = {2/(4 − π)} 1/3 . Under this denomination, the information matrix can be written as where D is a matrix that represents the derivative of the parameters of the standard representation (λ 1 , λ 2 and λ) regarding to the new parameters (µ, σ and γ 1 ). It also follows that the information matrix converges to a diagonal matrix Σ −1 c = diag(σ 2 , σ 2 /2, 6) when λ → 0. This guarantees the existence and uniqueness of the MLEs of λ 1 and λ 2 for each fixed value of λ.

Extension to Censored Data
Based on the goodness of the APST distribution to fit asymmetric and heavy-tailed data, in this section we introduce the censored APST model which we will be denote by CAPST.

Definition 3.
Suppose that the random variable Y follows APST distribution, and consider a random sample Y = (Y 1 , Y 2 , . . . , Y n ) where only the Y i values greater than a constant k are recorded. In addition, for values Y i ≤ k only the value of k is recorded. Therefore, for i = 1, 2, . . . , n, the observed values Y o i can be written as The resulting sample is said to be a censored APST, and we say that Y is a censored random variable APST. We will use the notation Y ∼ CAPST(θ), where θ = (µ, σ, λ, α, ν) .

From Definition 3 it follows that
. For convenience, we choose to work with the case of left-censored data; however, the followings results can be extended to other types of censorship.
The estimates of the parameters of the model can be obtained via maximum likelihood method, where the log-likelihood function is given by where x i = (y i − µ)/σ; ∑ 1 and ∑ 0 are the sum over censored individuals and uncensored individuals, respectively; and n 1 is the number of uncensored individuals.

Real Data Applications
In this section, we illustrate the applicability of the proposed model in Section 2 by analyzing two data sets. We use the statistical software R [15], version 3.5.3 with the package maxLike for maximizing the corresponding likelihood functions. For comparing purposes of various models, the AIC Akaike [20], BIC Schwarz [21], and corrected AIC (CAIC) Bozdogan [22] information criteria were used.

Application 1: Volcano Heights Data
Consider the data set related to heights of 1520 volcanoes in the world which is available in website dx.doi.org/10.5479/si.GVP.VOTW4-2013 [23]. Table 2 presents the summary statistics for the data set. It can be noted that the asymmetry and kurtosis indices seem to indicate that the use of an asymmetric and heavy-tailed model is appropriate to analyze this data set. We analyzed these data by fitting the Student-t, ST, PT, and APST distributions.  Table 3 shows the parameter estimates, together with their corresponding standard errors (SE). Note that the values of the standard errors of the µ and σ estimates for the APST model are smaller than the corresponding standard errors of the respective parameters for the Student-t, ST, and PT models. Table 3 also presents some model selection criteria, together with the values of the log-likelihood. The AIC, BIC, and CAIC criteria indicate that the APST model seems to provide better fit to the volcanoes heights data than the T, ST, and PT models, supporting the asymmetry assertion of the volcano's heights variable. Figure 2 shows the graphs QQplot of the fitted models. It can be clearly seen from the figure that the APST model fits the data better than the Student-t, ST, and PT models. In addition, we can use the likelihood ratio (LR) test statistic to conform our claim. To do this, we consider the following hypotheses, H 0 : (λ, α) = (0, 1) (T(µ, σ, ν)) v.s H 1 : (λ, α) = (0, 1) (APST(µ, σ, λ, α, ν)), The value of the LR test statistic is −2 log(Λ) = −2 T (θ) − APST (θ) =134.823 and comparing this quantity with χ 2 2 =5.9914, the null hypotheses is rejected. The APST model is also compared with the ST and PT models by considering the hypotheses H 01 : α = 1 (ST(µ, σ, λ, ν)) v.s H 11 : α = 1 (APST(µ, σ, λ, α, ν)),

Application 2: Stellar Abundances Data
The second data set is related to measurements for 68 solar-type stars, which are available in the package astrodatR of the software R [24] under the name Stellar abundances. These data were previously analyzed Mattos et al. [25] by using the Scale Mixture of Skew Normal Censored Regression (SMSNCR) models. We take only the response variable: log N(Be), which represents the log of the abundance of beryllium scaled to Sun's abundance (i.e., the Sun has log N(Be) = 0.0) In astronomical research, a previously identified sample of objects (stars, galaxies, quasars, X-ray sources, etc.) is observed at some new wavebands. According to Feigelson [24], due to limited sensitivities, some objects may be undetected, leading to upper limits in their derived luminosities. For this dataset we have 12 left-censored data points, i.e., 12 undetected beryllium measurement, that represents 19.35% of observations. Table 4 presents the ML estimates for the parameters of the censored Studen-t (CT), censored skew-t (CST), censored power-t (CPT), and censored alpha-power skew-t (CAPST) models, together with their corresponding standard errors. Table 4 also compares the fit of the four models using the model selection criteria (AIC, CAIC and BIC). Note that, again, the CAPST model with heavy tails have better fit than the CT, CST, and CPT models.
To identify atypical observations and/or model mispecification, we analyzed the transformation of the martingale residual, r MT i , proposed in Barros et al. [26]. These residuals are defined by where r Mi = δ i + log S(y i ;θ) is the martingal residual proposed by Ortega et al. [27], where δ i = 0, 1 indicates whether the i-th observation is censored or not, respectively; sign(r Mi ) denotes the sign of r Mi ; and S(y i ;θ) = Pθ(Y i > y i ) represents the survival function evaluated at y i , whereθ are the MLE for θ. The plots of r MT i with generated confidence envelopes are presented in Figure 4. From this figure, we can see clearly that the CST, CPT, and CAPST models fit better the data than the CT model, since, in that cases, there are not observations which lie outside the envelopes. The Figure 5 shows the graph of the densities of the different models fitted to the stellar abundances data. From the figure, the CAPST model seems to fit better the stellar abundances data than CT, CST and CPT models.

Conclusions
In this work, a new asymmetric model has been introduced. It is based on the combination of skew-t [1] and power-t [2] models. The new model presents greater ranges of asymmetry and kurtosis, which is very useful for modeling skewed and heavy-tailed data. The problem of estimating the parameters in the model is dealt by using the maximum likelihood approach which is also used for developing large sample properties for the estimators. The elements of the observed information matrix are analytically obtained. The likelihood ratio statistics can be used for testing the APST null hypothesis since the Student-t, ST, and PT models are special cases of the model entertained. Two applications to volcano heights data and stellar abundances data indicate that the proposed model can be a useful alternative to the ST and PT models.