The Gamma Generalized Pareto Distribution with Applications in Survival Analysis

We study a three-parameter model named the gamma generalized Pareto distribution. This distribution extends the generalized Pareto model, which has many applications in areas such as insurance, reliability, finance and many others. We derive some of its characterizations and mathematical properties including explicit expressions for the density and quantile functions, ordinary and incomplete moments, mean deviations, Bonferroni and Lorenz curves, generating function, Rényi entropy and order statistics. We discuss the estimation of the model parameters by maximum likelihood. A small Monte Carlo simulation study and two applications to real data are presented. We hope that this distribution may be useful for modeling survival and reliability data.


Introduction
The Pareto distribution is a very known statistical model, widely used for accommodate heavy-tailed distributions and many of its applications in areas such as economics, biology and physics can be found in the literature (Alzaatreh, A., et al., 2012). An important extension of the Pareto model is the generalized Pareto (GP for short) distribution. It is believed that this distribution was pioneered by Pickands, J.(1975) in the extreme values context as the distribution of the sample of surpluses above a certain high level. Some of its applications are addressed to areas such as insurance, reliability, finance, meteorology and environment, among others (De Zea Bermudez, P. & Kotz, S., 2010). In a recent paper, (Nadarajah, S. & Gupta, A. K., 2007) list many practical situations in which the GP distribution has been applied. Here, we mention some of these: lifetime data analysis, coupon collector's problem, analysis of radio audience data, analysis of rainfall time series, regional flood frequency analysis, drought modeling, value at risk and measuring liquidity risk of open-end funds, among others.
The cumulative distribution function (cdf) of the GP distribution is given by The corresponding probability density function (pdf) is given by where ξ ∈ R and σ > 0 are the shape and scale parameters, respectively. The support of g(x; ξ, σ) is x > 0 for ξ ≥ 0 and 0 < x < −σ/ξ for ξ < 0.
Here, we adopt the same parametrization considered by (Song, J. & Song, S., 2012) for the GP distribution. The case ξ = 0 reduces to the exponential distribution. In general, it is considered that the scale parameter has a direct effect on the tails of this distribution, so that it has heavy-tailed when ξ > 0, medium-tailed when ξ = 0, and short-tailed when ξ < 0 (Song, J. & Song, S., 2012).
In recent years, several studies on generalizations of continuous distributions by introducing new parameters have been considered. Particularly, great attention has been given to generalizations by using special generators (Tahir, M. & Nadarajah, S., 2013). In this context, some classes well-established in the literature are: the Marshall-Olkin class proposed by (Marshall, A. W. & Olkin, I., 1997), the beta class pioneered by (Eugene, N., et al., 2002), the gamma class proposed by (Zografos, K. & Balakrishnan, N., 2009) and (Ristić, M. M. & Balakrishnan, N., 2012) and the Kumaraswamy class defined by (Cordeiro, G. M. & De Castro, M., 2011).
For a given parent distribution with cdf G(x) and pdf g(x), x ∈ R, (Zografos, K. & Balakrishnan, N., 2009) proposed the gamma-G family with an extra shape parameter a > 0 and cdf F(x) given by where denotes the lower incomplete gamma function and Γ(·) is the gamma function. The pdf of this family is given by Its hazard rate function (hrf) can be expressed as where denotes the upper incomplete gamma function. Note that Γ(a) = γ(a, z) + Γ(a, z). Several mathematical properties about this class were studied by .
Due to its wide applicability, some recent works extended the GP distribution using special generators to obtain more flexibility. Two important extensions are the beta GP and Kumaraswamy GP distributions investigated by ( Mahmoudi, E., 2011) and (Nadarajah, S. & Eljabri, S., 2013), respectively. Based on a similar motivation, we study here another promising extension of the GP distribution using the gamma generator called the gamma generalized Pareto (GGP for short) distribution. We present some of its mathematical properties and emphasize its application in the context of survival analysis.
The paper is organized as follows. In Section , we provide the cdf, pdf and hrf of the GGP distribution. Some important results about transformations involving this distribution are addressed in Section . Some of its mathematical properties are derived in Section including a mixture representation for its pdf and explicit expressions for the quantile function (qf), ordinary and incomplete moments, mean deviations, Bonferroni and Lorenz curves, generating function, Rényi entropy and order statistics. In Section , we use the maximum likelihood method to estimate the parameters of the GGP distribution. A simulation study and two applications are addressed in Section . Section provides some concluding remarks.
The corresponding pdf has the form International Journal of Statistics and Probability Vol. 6, No. 3;2017 and its hrf is given by Setting ξ = 0 in equation (6), we have the gamma exponential distribution, which was studied by (Risti'c, M. M. & Balakrishnan, N., 2012) and (Nascimento, D. C., 2014), whereas setting ξ < 0 leads to a support for X that depends on unknown parameters. For this reason, we consider only the GGP model with ξ > 0 and positive support. We emphasize that (Nadarajah, S.& Eljabri, S., 2013) gave the cdf (Dahiya, R. C., and Gurland, J., 1972) outside the context of the gamma generator but they did not study any of its properties.
Henceforth, a random variable X having cdf (6) is denoted by X ∼ GGP ( a, ξ, σ ) and we write F(x) = F(x; a, ξ, σ) in order to simplify the notation. Plots of the pdf and hrf of X for selected parameter values are displayed in Figures and , respectively.

Some Important Results
, where a, ξ and σ are positive real numbers. Then, the random variable Y = log [ has the log-gamma distribution with parameters a and 1/ξ, namely Y ∼ log-gamma ( a, 1/ξ ) .

Proof. Consider the function h(x)
and its inverse function h −1 (x) = σ [exp(e x ) − 1]/ξ. The density function of Y, say f Y (y), is obtained from the transformation method and equation (7) gives Further, if W ∼ gamma(a, 1/ξ), then log(W) ∼ log-gamma(a, 1/ξ) and its density function agrees with the last equation.
In the context of survival analysis, the Weibull distribution is commonly used to model data with monotone hrf. An important generalization of the Weibull model is the gamma-exponentiated Weibull (GEW) distribution proposed by (Pinho, L. G., et al., 2012). The GEW model can be characterized as follows.

Proof. Consider the function
From equation (7) and using the transformation method, the density function of Y, f Y (y), can be expressed as The last equation agrees with equation (5) of (Pinho, L. G., et al., 2012) by taking δ = a, α = 1/ξ, k = 1 and λ = ξ.
International Journal of Statistics and Probability Vol. 6, No. 3;2017 Theorem 2.3 [Transformation] Let Y be a gamma distributed random variable with parameters a and 1/ξ. Then, the random variable X = σ(e Y − 1)/ξ has the GGP distribution with parameters a, ξ and σ.

Mathematical Properties
Many structural properties of the GGP distribution can be derived using the concept of exponentiated distributions. The class of exponentiated distributions was initially studied by (Lehmann, E. L., 1953), who defined the so-called Lehmann alternatives. However, this class has received special attention in the literature over the last twenty years and several papers have been published. An excellent review of these publications is provided by (Tahir, M. & Nadarajah, S., 2013).

A Useful Representation
For any baseline cdf G(x), a random variable Y is said to have the exponentiated-G ("exp-G" for short) distribution with power parameter α > 0, say Y ∼ G α , if its cdf and pdf are given by Using a result pioneered by , we can express the density of the GGP distribution in terms of exp-GP densities. Thus, for a real non-integer a > 0, equation (7) can be expressed as where, the coefficients b k are given by and the constants p j,k can be determined recursively by for k = 1, 2, . . . and p j,0 = 1. Here, h a+k (x) denotes the exp-GP density with power parameter a + k given by Equation (9) reveals that the GGP density function is a linear combination of exp-GP densities. This result is important to derive some mathematical properties of X such as the ordinary and incomplete moments, generating function and mean deviations from those of the exp-GP distribution.

Quantile Function
The qf is obtained as the inverse of the cumulative function. It follows from equation (6) that the qf of the GGP distribution can be expressed as for 0 < u < 1, where W −1 (a, u) is the qf of the gamma distribution with shape parameter a and scale parameter one evaluated at u. The function W −1 (a, u) can be expressed as a power series expansion (Cordeiro, G. M., 2013) where t 0 = 0, t 1 = 1 and any coefficient t i+1 for i ≥ 1 is determined by the cubic recurrence equation It follows directly from (12) that the median of X is simply respectively.

Moments
The nth moment of X can be determined based on (9) as Setting u = (1 + ξx/σ) −1/ξ , E(X n ) can be expressed as Using the binomial expansion for (u −ξ − 1) n and the generalized binomial expansion for (1 − u) a+k−1 , the nth moment of X reduces to

Generating Function
The moment generating function (mgf) of X can be obtained using the fact that the GGP density function is a linear combination of exp-GP densities. Thus, Setting u = (1 + ξx/σ) −1/ξ , M(t) can be expressed as Using the generalized binomial expansion for (1 − u) a+k−1 , M(t) reduces to Thus, it follows that (for t < 0) ] .

Incomplete Moments
Frequently, in applied works it is of interest to know the nth incomplete moment of a random variable. The nth incomplete moment of X is given by T n (z) = ∫ z −∞ x n f (x) dx, and using equation (9), we can write Setting u = (1 + ξx/σ) −1/ξ , T n (z) reduces to Using the binomial expansion for (u −ξ − 1) n and the generalized binomial expansion for (1 − u) a+k−1 , T n (z) reduces to A particularly important case is the first incomplete moment. From equation (13), we have An important application of equation (14) refers to the mean deviations of X about the mean and the median. These quantities are defined by respectively, where µ ′ 1 = E(X), x 1/2 is the median of X that follows from (12) as x 1/2 = Q(1/2), F(µ ′ 1 ) is easily obtained from (6) and the quantities T 1 (µ ′ 1 ) and T 1 (x 1/2 ) are given by (14). Equation (14) can also be used to determine the Bonferroni and Lorenz curves, which have several applications in economics, reliability, demography and medicine. These curves are defined by B(π) = T 1 (q)/(πµ ′ 1 ) and L(π) = T 1 (q)/µ ′ 1 , where q = Q(π) is calculated by (12) for a given probability π.

Rényi Entropy
The Rényi entropy of X is given by where γ > 0, γ 1 and f (x) is the pdf (7).
For a general cdf G(x) satisfying equation (3),  proved that where and p j,k is defined by (10).
For the GP cdf, the integral given in equation (17) becomes Setting u = ( 1 + ξx/σ ) −1/ξ , the integral I k can be expressed as Then, equation (16) can be reduced to .
Then, from equations (15) and (18), the entropy of X can be given as

Order Statistics
Here, we derive an explicit expression for the density of the ith order statistic X i:n , say f i:n (x), in a random sample of size n from the GGP distribution. Using results by , we can write where Here, h a( j+i)+r+k (x) denotes the pdf of the exp-GP distribution with power parameter a( j + i) + r + k and it can be easily obtained by a reparameterization of (11).
Equation (19) reveals that the density of the ith order statistic of the GGP distribution is a mixture of exp-GP densities. This equation can be used, for example, to obtain the moments and mgf of the GGP order statistics.
The sth moment of X i:n is given by Using similar results to those presented in Section , we have The mgf of X i:n can be expressed as which, using similar results to those presented in Section , reduces (for t < 0) to ] .

Estimation
In this section, we emphasize how to obtain the MLEs of the parameters a, ξ and σ of the GGP distribution. Let x 1 , . . . , x n be a sample of size n from this distribution. The log-likelihood function for the vector of parameters θ ⊤ = (a, ξ, σ) ⊤ can be expressed as The elements of the score vector are given by where ψ(a) = d log[Γ(a)]/da and z i = 1 + ξx i /σ.
The MLE θ of θ is obtained by solving simultaneously the nonlinear equations U ξ (θ) = 0, U σ (θ) = 0 and U a (θ) = 0. These estimators can be obtained numerically by maximizing the log-likelihood function by means of algorithms for non-linear optimization. For confidence intervals and hypothesis tests, it is important to determine the observed information matrix. Particularly, for the GGP distribution, the 3 × 3 observed information matrix is J(θ) = {−U rs }, where U rs = ∂ 2 ℓ(θ)/(∂θ r ∂θ s ) for r, s ∈ {a, ξ, σ}. The elements of J(θ) are listed in the Appendix.

Applied Results
In this section, a simulation study and two applications to real data are considered. The parameter estimation was made using the R software (www.r-project.org) by the maxLik package with simulated annealing method, which is useful for finding a global maximum when the objective function has some local maximums.

Simulation Study
Here, we present the results of a Monte Carlo simulation study, which was carried out to evaluate the performance of the MLEs for the parameters a, ξ and σ of the GGP distribution in finite samples. The results are obtained from 10, 000 Monte Carlo replications. We take samples of sizes n = 100, 200, 300, 500 and 800. In the data generation process, we consider equation (12) and the true parameters values a = 2, ξ = 0.5 and σ = 1.5.
In Table 1, we provide the means, standard error estimates, root mean square errors (RMSEs) and mean absolute errors (MAEs) of a, ξ and σ for the GGP model. Based on the results, we can verify that the estimators are biased even for moderate samples. Thus, the use of bias-correction techniques can improve the performance of these estimators in finite samples.

2 Application to Real Survival and Reliability Data
Here, we present two applications to real data to illustrate the potentiality of the GGP model. To compare the performance of the GGP model, we use two generalizations of the GP distribution well-established in the literature. We consider the four-parameter beta generalized Pareto (BGP) model proposed by (Mahmoudi, E., 2011) and the Kumaraswamy generalized Pareto (KGP) model pioneered by (Nadarajah, S. & Eljabri, S., 2013). The density functions of the KGP and BGP models are given by and respectively, where z = 1 + ξx/σ and B(·, ·) is the beta function.

Conclusion
In this paper, we study a three-parameter model named gamma generalized Pareto (GGP) distribution, which consists of a major extension of the generalized Pareto (GP) distribution pioneered by (Pickands, J., 1975). We give some results International Journal of Statistics and Probability Vol. 6, No. 3;2017 Time between failures connecting the GGP model with other well-established distributions in the literature, such as the gamma and log-gamma distributions. We also provide several mathematical properties of the GGP model including explicit expressions for the density and quantile functions, ordinary and incomplete moments, mean deviations, Bonferroni and Lorenz curves, generating function, Rényi entropy and order statistics. We discuss the maximum likelihood method to estimate the model parameters and provide the elements of the score vector and the observed information matrix. We also present a Monte Carlo simulation study to evaluate the performance of the maximum likelihood estimators for the GGP model. Finally, two applications illustrate the potentiality of the GGP distribution to fit survival data.