Approximate Simulation Techniques and Distribution of an Extended Gamma Process

In reliability theory, many papers use a standard Gamma process to model the evolution of the cumulative deterioration of a system over time. When the variance-to-mean ratio of the system deterioration level varies over time, the standard Gamma process is not convenient any more because it provides a constant ratio. A way to overcome this restriction is to consider the extended version of a Gamma process proposed by Cinlar (J Appl Probab 17:467–480, 1980). However, based on its technicality, the use of such a process for applicative purpose requires the preliminary development of technical tools for simulating its paths and for the numerical assessment of its distribution. This paper is devoted to these two points.


Introduction
Safety and dependability is a crucial issue in many industries, which has lead to the development of a huge literature devoted to the so-called reliability theory. In the oldest literature, the lifetimes of industrial systems or components were usually directly modeled through random variables, see e.g. Barlow and Proschan (1965) for a pioneer work on the subject. Based on the development of on-line monitoring which allows the effective measurement of a system deterioration, numerous papers nowadays model the degradation in itself, which often is considered to be accumulating over time. This is done through the use of a non decreasing stochastic process (X t ) t≥0 (say), where the system lifetime τ typically corresponds to the hitting time of a given failure threshold L: τ = inf(t > 0 : X t > L).
A major quantity of interest for such a system is its reliability function R(t), which corresponds to the probability that, without any repair, the system is still in order at time t: where F X t and f X t stand for the cumulative distribution function (cdf) and probability density function (pdf) of X t , respectively. Another important indicator for prediction purpose is the system Mean Time Before Failure (MTBF), with Most of other quantities of interest (Mean Residual Lifetime, cost/profit functions, . . . ) may also be expressed with respect to the pdf or cdf of X t , e.g. see Rausand and Høyland (2004). Before being able to use a stochastic process (X t ) t≥0 for deterioration modeling purpose and make prediction over the system future behavior, there hence is a crucial need to develop tools for computing its pointwise pdf and cdf at time t and/or for simulating trajectories of (X t ) t≥0 , which provides another way to compute the aimed quantities through Monte-Carlo simulations.
Since seminal papers (Abdel-Hameed 1975) and (Ç inlar et al. 1977), the most widely used process to model the phenomena of cumulative degradation is the (standard) Gamma process, see van Noortwijk (2009) for a comprehensive presentation and overview of applications of the Gamma process to reliability theory, including more than 100 references and both scientific and empirical justification of its use. However, a notable restriction of a standard Gamma process is that its variance-to-mean ratio is constant over time, which may be restrictive within an applicative context, see Guida et al. (2012) for a real data set consisting of "sliding wear data of four metal alloy specimens" and where there is "empirical evidence that the variance-to-mean ratio is not a constant but varies with [time]". To overcome this drawback, we propose to use an Extended Gamma Process (EGP), which was introduced mostly simultaneously by Ç inlar (1980) and Dykstra and Laud (1981). Note that EGPs are also called weighted Gamma processes in the literature (Ishwaran and James 2004). An EGP is a non decreasing process with independent increments, which can be constructed as a stochastic integral with respect to a standard Gamma process. If the EGP has been used for Bayesian modeling of the hazard function (Dykstra and Laud 1981; Ishwaran and James 2004;Laud et al. 1996), up to our knowledge, it has not been much studied for cumulative degradation modeling, except in Guida et al. (2012) in a simplified setting.
A standard Gamma process is characterized by a shape function and a constant scale parameter. For an EGP, the scale parameter may vary over time. This allows for more flexibility than its standard version, for modeling purpose. However, there is a cost and the use of an EGP presents some technical difficulties. Firstly, except for specific cases, the exact simulation of such a process is generally impossible. Secondly, there is no explicit formula for the probability distribution of an EGP (which is not Gamma). These technical difficulties have lead Guida et al. (2012) to use a discrete version of an EGP. We here propose to deal with the original continuous time version.
The aim of the present paper is to develop technical tools necessary to the practical use of an EGP for cumulative degradation modeling (or for any other purpose). More precisely, we focus on the simulation of approximate sample paths and on the numerical assessment of both its pdf and cdf. Our contribution is threefold. Firstly, the series representations provided by Rosinski (2001) for standard Gamma processes are extended to EGPs. Such representations are used to generate approximate sample paths. The quality of the approximation is studied and compared to the method proposed by Ishwaran and James (2004), which is based on an alternate approximate representation of an EGP. Secondly, following the results of Veillette and Taqqu (2011) for Laplace transform inversion through Post-Widder formulas, explicit asymptotic expressions are provided for both pdf and cdf of an EGP. Thirdly, a discretization method is proposed, based on the approximation of a general EGP by a specific EGP with a piecewise constant scale function, which is easy to deal with. The discretization method allows both to simulate approximate sample paths of a general EGP and to compute approximations of its pdf and cdf. Convergence results are obtained for the approximate discretized EGP towards the initial general EGP and the quality of the approximation is studied. The method allows to compute the cdf of a general EGP at a known and adjustable precision. Up to our knowledge, a similar result was not available in the previous literature.
The paper is organized as follows. Section 2 gives an overview of the EGP. In Section 3, the approximate representation of Ishwaran and James (2004) is presented as well as four series representations of an EGP. Section 4 is devoted to the approximation of the pdf and cdf through Laplace transform inversion. Section 5 introduces and studies the discretization method. The different methods are illustrated through numerical experiments in Section 6. Conclusions are formulated in Section 7.

Definition of an EGP and First Properties
Let a : R + → R + be a measurable, increasing and right-continuous function with a(0) = 0 and let b 0 > 0. Recall that a standard (non homogeneous) Gamma process 0 (a(t), b 0 ) with a(.) as shape function and b 0 as (constant) scale parameter is a stochastic process with independent, non-negative and Gamma distributed increments. Its pdf at a specified time point t is given by e.g. see Abdel-Hameed (1975). Now, let b : R * + → R * + be a measurable positive function such that, for all t > 0: Following (Ç inlar 1980;Dykstra and Laud 1981), the process X = (X t ) t≥0 is said to be an EGP with shape function a(.) and scale function b(.) (written X ∼ (a(t), b(t)) in the sequel) if it can be represented as a stochastic integral with respect to a standard Gamma process (Y t ) t≥0 ∼ 0 (a(t), 1): and X 0 = 0. If b(.) is constant and equal to b 0 , the EGP simply reduces to a standard Gamma process 0 (a(t), b 0 ). An EGP can be proved to have independent increments and its distribution to be infinitely divisible (Ç inlar 1980). Also, an explicit formula is available for the Laplace transform of an increment, with for all t, λ ≥ 0 and h > 0. Conversely, a stochastic process with independent increments and Laplace transform of an increment provided by (3) can be proved to be an EGP ∼ (a(t), b(t)).
Finally, the mean and variance of an EGP are given by Remark 1 Let X ∼ (a(t), b(t)). So far, the shape function a was supposed to be rightcontinuous. Then, following (Ç inlar 1980), the function a can be written as a = a c + a d , where a c is the continuous part of a and a d its jump part. The EGP X is then the sum of two independent EGPs X c ∼ (a c (t), b(t)) and X d ∼ (a d (t), b(t)). This second EGP simply reduces, for each t ≥ 0, to a sum of independent Gamma variables, for which the distribution is known in full form (see Section 5). This second EGP is hence easy to deal with and can be omitted without loss of generality. Thus, the shape function a is assumed to be continuous in the sequel.
For a better understanding of the possible evolution of an EGP over time, we now look at its asymptotic behavior and we set which exists, due to the non decreasingness property of an EGP.
Proposition 1 Let X = (X t ) t≥0 be an EGP with X ∼ (a(t), b(t)). We then have: and the following results: then X ∞ is almost surely infinite.
Under assumption (6), we consequently have E(e −X ∞ ) = 0. This implies that e −X ∞ = 0 almost surely and the result.
Remark 2 In the case where the dichotomy P (X ∞ = ∞) ∈ {0, 1} is not valid any more and both P (X ∞ = ∞) > 0 and P (X ∞ < ∞) > 0 are possible, as is illustrated in the following example (third case).

Example 1
The different types of behavior are illustrated in Figs. 1, 2 and 3 (left) where a few sample paths are plotted for three different EGPs (using the discretized rate method from Section 5): The process X (1) (resp. X (2) ) corresponds to the first (resp. second) case of Proposition 1 whereas X (3) corresponds to case (7). As expected, all trajectories of Fig. 1 are stabilizing whereas all trajectories of Fig. 2 are exploding. In Fig. 3, some  In all the following, X = (X t ) t≥0 stands for an EGP ∼ (a(t), b(t)) and Y = (Y t ) t≥0 stands for a standard Gamma process 0 (a(t), 1), without any further notification. We recall that the shape function a(t) is assumed to be continuous.

Simulation Techniques of an EGP
The goal of this section is to present methods for generating approximate sample paths of an EGP on a given compact set [0, T ]. The section is divided into two parts. The first part quickly recalls the method proposed by Ishwaran and James (2004). In the second part, the series representations of a standard Gamma process provided in Rosinski (2001) are extended to the case of an EGP. Algorithms based on these series representations are further provided in Section 6 for approximate sample paths generation.

Ishwaran and James's Method
The method proposed by Ishwaran and James (2004) roughly boils down to some discretization of the stochastic integral (2) for t ∈ [0, T ], where T > 0. More specifically, for each t ∈ [0; T ] and each K * ∈ N, they consider: where K * , 1), independent from the V k 's. Ishwaran and James (2004) −→ (X t ) t≥0 weakly when K * → ∞ in the Skorokhod space D (0, ∞) of right-continuous functions with left-side limits. (To be more specific, they show the weak convergence of the random measure with corresponding cumulative process defined by (8) towards a weighted gamma measure with shape measure a(dt) and scale function b(t). The weak convergence of the cumulative process (8) towards an EGP is then a consequence of (Daley DJ and Vere-Jones 2007, Lemma 11.1.XI)). One can easily check that

Series Representations
The following results extend those from (Rosinski 2001, Section 6) devoted to standard Gamma processes, see Rosinski (2001) for more references and historical remarks. (1) Bondesson's series representation (Bond):

Proposition 2 Let T > 0 and let (U n ) n≥1 be the points of a homogeneous Poisson process
where {W n } n≥1 is a sequence of i.i.d exponential r.v.s with mean 1, (2) Rejection's series representation (Rej): Thinning's series representation (Thin): where {W n } n≥1 is a sequence of i.i.d exponential r.v.s with mean 1, (4) Inversion of Lévy measure's series representation (ILM): Proof The four series representations can be written in the same manner: and an appropriate choice for H : R 2 + → R + . As already mentioned in Section 2, (X t ) 0≤t≤T is an EGP if its increments are independent and the Laplace transform of an increment can be written as in (3). Let us first note that (U n , V n , W n ) n≥1 can be seen as the points of a Poisson random measure N on R 3 + with intensity (see for example (Ç ınlar E 2011, Corollary 3.5, p. 265)). Considering Now, based on the Laplace functional theorem for a Poisson random measure (Ç ınlar E 2011, Theorem 2.9, p. 252) for the third line, we have for all t ∈ [0, T ] and all h, λ > 0: where From (3) and (14), it then suffices to check that Q(λ, v) = log 1 + λ b(v) for each of the four series representations.
2) Rejection: setting z = 1 exp(u)−1 in the last line. Equation (15) is a Frullani's integral, which can classically be computed by differentiating with respect of λ: which provides the result. 3) Thinning: We recognize (15), which allows to conclude. 4) Inversion of Lévy measure: We recognize (15) again, which allows to conclude.
Approximate simulation of (X t ) 0≤t≤T is done by truncating the infinite series. As is classically done in the Lévy case, we retain only the points of the Poisson process M which belong to a compact set [0, B], where B > 0: for all t ∈ [0, T ], where we use the notations of (13). This allows to better control the truncation error than simply retaining a fixed number of terms in the series, see e.g. Imai and Kawai (2013) for further details in the Lévy case. We also set for all t ∈ [0, T ], to be the remainder of the truncated series. By construction, both have independent increments. The moments ofX (B) t are given in the following proposition.

Proposition 3 Recall thatX (B)
t is the remainder of the truncated series (16). We have (1) Bondesson's series representation: (2) Rejection's series representation: Thinning's series representation: Proof Starting from (16) Proposition 3 allows to determine the minimal value of B which ensures a given relative precision on both mean and variance of X t for all t ∈ [0, T ] (analytically for Bondesson's series representation and numerically for the three other procedures). Asymptotic equivalents of the series remainders are provided in Table 1 when B tends to infinity, which are obtained through Taylor expansions.   (−B)) when B → ∞). However, the computing time of the Inversion of Lévy measure's method may be long because of the numerical computation of the inverse of the exponential integral function, for which no explicit formula is available. As far as the Thinning's method is concerned, it clearly converges more slowly than the others. Thus, Rejection's and Bondesson's methods seem to be the most suitable ones among the four series representations to approximate the paths of an EGP.

Post-Widder Formulas
As already noted in Section 2, the distribution of an EGP is infinitely divisible and its Laplace transform is available in full form (see (3)). This allows to use the method from Veillette and Taqqu (2011) to compute the pdf and cdf of an EGP, by inverting its Laplace transform through the Post-Widder method. We use the following formulas from (Masol and Teugels 2010, Section 4) for the pdf and cdf (see also Veillette and Taqqu (2011) for the cdf): for all x ≥ 0, where L (k) X t denotes the k-th derivative of the Laplace transform L X t for k ≥ 1 and L (0) Let φ t be the Laplace exponent of X t , with Following Veillette and Taqqu (2011), Leibniz's formula provides: with φ (j ) t the j -th derivative of φ t . An expression of φ (j ) t is given in the following Lemma in the case of an EGP.
Lemma 1 Let λ ≥ 0, j ≥ 1 and let X ∼ (a(t), b(t)) . Then φ (j ) Proof Based on (3) and (15) and setting u = z/b(s) in the third line we have: for all j ≥ 1. This provides: renormalizing the pdf of 0 (j, b(s) + λ) to get the last line.
Based on (19) and Lemma 1, L (k) X t (λ) can be recursively computed through This allows to compute an approximation of the pdf and cdf of an EGP through (17-18), taking H large enough.

Discretized Rate Function Method
In case of a piecewise constant scale function b(.), the process (X t ) t≥0 can easily be constructed from standard Gamma processes. The simulation of its paths is hence immediate (see e.g. van Noortwijk (2009)). Also, the random variable X t simply is the sum of standard Gamma variables, and different tools are available in the literature to compute both its pdf and cdf (see Nadarajah (2008) for a review). Based on this, the aim of this section is to propose an approximation of an EGP with a general scale function by another EGP with a piecewise constant scale function.
Remark 3 Based of the definition (2) of an EGP, we prefer to construct an approximation

Quality of the Approximation of X by X (ε)
We first look at the moments of the residual part in the approximation of X t by X (ε) t , for t ∈ [0, T ].
We easily derive that from where we derive (29). Also, one can check that (4) is still valid forX We now provide some convergence results.
Proposition 5 Let X ∼ (a(t), b(t)), (ε n ) n≥1 be a sequence of positive real numbers and let (X (ε n ) ) n≥1 be a sequence of EGPs with shape function a(.) and scale function b (ε n ) (.). Suppose that assumption (21) is satisfied.
Proof 1) According to (31), we have: based on Markov inequality for the last line. We derive that n≥1 P sup which concludes the first point, based on e.g. (Ç ınlar E 2011, Proposition 2.7, p. 98). 2) In the same way, we have: which allows to conclude.
Remark 4 Note that a consequence of point 2 is that, under the same assumptions, X

Approximation of the pdf and cdf of X (ε) t (and of X t )
with P = P ( , t) the single integer such that l P < t ≤ l P +1 , α p = a(l p+1 ) − a(l p ) for p = 0, . . . , P − 1 and α P = a(t) − a(l P ).
Taking K large enough, (32-33) provide approximations for the pdf and cdf of X (ε) t , and consequently of X t .
Using that | sin(.)| ≤ 1 and that (1 + (u/b p ) 2 ) α p /2 ≥ (u/b p ) α p , we easily get an upper bound for the approximation error, with: for all K ∈ N, due to P p=0 α p = a(t) for the second-to-last line.

Bounds for the cdf of X t
Consider two other piecewise constant approximations of 1 b(t) defined by and Proof The bounds (37) are easily obtained starting from Equation (36) and using that for the upper bound, and similar arguments for the lower one. Inequality (38) is a direct consequence of (37).
The previous theorem provides computable bounds for the cdf F X t of an EGP. These bounds can be made as tight as necessary, taking ε small enough and K large enough. This method hence provides a way to numerically assess the cdf F X t at a known precision. This is used in Section 6 to get numerical reference results. Remark also that the two piecewise constant approximations b (ε,−) (.) and b (ε,+) (.) satisfy (27). The theoretical results of Subsection 5.2 are then still valid for X (ε,−) and X (ε,+) .
Finally, remark that for any fixed s > 0, the increment process (X t+s − X s ) t≥0 still is an EGP ∼ (a(t + s) − a(s), b(t + s)) (based on its independent increments and on its pointwise Laplace transform provided by (3)). Thus, all the results of the paper can be used to compute both pdf and cdf of an EGP increment, or to simulate it. (For simulation purpose, one can also only retain the terms for which V n ∈ (s, s + t] in the series representations from Proposition 2).

When Assumption (21) is not Satisfied
Assumption (21) allows to construct the piecewise constant function b (ε) (.) specified in Subsection 5.1. When this assumption is not satisfied, alternatives for the construction of b (ε) (.) have to be considered.
• An interesting case is when function b(.) satisfies lim t→0 + b(t) = 0. If the following assumption is checked, then the piecewise constant approximation can be constructed from l 0 = 0, l 1 such that: and the l i 's, for i = 2, . . . , n(ε) + 1, recursively defined according to (23). Note that in this case, inequality (27) is not satisfied any more. • If function b(.) is not continuous on (0, T ] but only piecewise continuous, the procedure described in Subsection 5.1 can be applied on each continuous section and inequality (27) is still valid.

Summary and Algorithms
To sum up the above sections, three main methods are considered to simulate an EGP: • Method 1: Discretization of stochastic integral (DSI) from Subsection 3.1, • Method 2: Series representations from Subsection 3.2 : ILM, Bond, Thin, Rej, • Method 3 : Discretized rate function (DRF) from Subsection 5.1.
Furthermore, two approaches to approximate the pdf and cdf of an EGP are proposed: • Approach 1: Post-Widder formula from Section 4, • Approach 2: DRF method + approximation of the cdf and pdf through the infinite integral formulas from Subsection 5.3. This method is referred to as DRFI in the sequel.
For the DRFI method, the cdf of X t is approximated by F as provided by (33), where we recall that X (ε) t is defined by (26). There hence are two parameters to define (ε and K). All these approaches are investigated in this section. We set X ∼ (a(t), b(t)) to be an EGP, where a(t) is assumed to be one-to-one. We first state the algorithms resulting from the three methods previously mentioned to simulate N approximate sample paths of (X t ) t∈ [0,T ] . Note that for Method 2, only Bondesson's series representation algorithm is given.

Approximate the rate function using (22).
For i in 1 : N , repeat the following steps : 4. Simulate Y j according to the Gamma distribution 0 a(l j +1 ) − a(l j ), b (ε) (l j ) for j = 0, . . . , n(ε) − 1 and Y n(ε) according to the distribution

Comparison of the Simulation Methods
The different approaches for the sample paths generation of an EGP are compared. For each method, N = 10 5 approximate sample paths are generated and the relative errors between the theoretical and empirical mean, variance and Laplace transform are computed. Let us recall that the Laplace transform of X t fully characterizes its distribution. In order to scan a good part of this distribution, we introduce λ j = L −1 X t (0.01 * j), j = 1, . . . , 100 and we define the relative error on the Laplace transform by: whereL X (ε) t (λ j ) is the empirical Laplace transform of X (ε) t at point λ j . The previous simulation is repeated 500 times and for each method, we compute the mean and the standard deviations of the tree relative errors (mean, variance and Laplace transform) based on these 500 sets of N = 10 5 approximate sample paths.
As a first step, the four series representations of Method 2 are compared and the results provided in Table 2. As expected from Table 1, the thinning method seems less accurate than the others. Moreover, the computing time for the inverse Lévy measure method is thrice the computing time of Bondesson and Rejection methods for a similar precision. Thus, only Bondesson and Rejection approaches are maintained to continue the comparison.
Tables 3, 4 and 5 provide the results for the DSI method from Subsection 3.1, the two selected series representations (Bondesson and Rejection) and for the proposed approach of the discretized rate function for different shape and rate functions (see the legends of the tables). The parameters of the four methods were adjusted to have similar computing  times. The results for the discretized rate function are quite similar to the series representations methods, while the DSI method is sometimes less accurate for the variance or for the Laplace transform (see Table 4). The discretized rate method hence seems to behave as well as Bondesson's and Rejection methods for simulation purpose, whereas the DSI method seems a little below.

Comparison of the Approximation Methods for the cdf and pdf
We here compare the quality of the numerical assessment of the cdf/pdf of an EGP through both Post-Widder's formula and the DRFI method. We begin with the cdf, for which we recall that we are able to get reference results, with an exact control of the error (see Subsection 5.4). An EGP X ∼ (a(t), b(t)) is considered with a(t) = t and b(t) = (t + 1) 2 . Reference results are first computed for F X t (x) with x ∈ [0.01 : 0.1 : 4] and t = 2, up to a precision lower than 3.10 −6 , obtained for K = 800 and ε = 10 −6 . The results are provided in Table 6 for the mean absolute error on F X t (x) with respect to the reference results for x ∈ [0.01 : 0.1 : 4] and both Post-Widder's and DRFI methods. For H 100 in Post-Widder's method and ε = 10 −2 , K = 12 for the DRFI approximation, both methods provide similar results with similar computation times. However, the DRFI method easily yields more accurate results by decreasing ε and increasing K whereas the highest precision by Post-Widder's method is obtained for H 100 in the present example, with a clearly lower precision for H = 200. This problem of poor convergence when H → ∞ for Post-Widder's method has already been observed in (Masol and Teugels 2010).  DRFI (ε = 10 −3 , K = 800) < precision 30 Fig. 4 The cdf of X t as a function of x for a(t) = t, b(t) = (t + 1) 2 and t = 10 Keeping the same a(t) and b(t), three approximations of the cdf of X 10 are next plotted in Fig. 4: a non-parametric estimation obtained from a sample of size 10 5 generated by Rejection method, Post-Widder's estimation (H = 100) and the DRFI approximation (ε = 0.001 and K = 100). We observe the good superposition of the DRFI approximation and of the non-parametric estimation. These results are in concordance with the convergence results obtained in Subsection 5.2, which imply the convergence of the cdf F (ε) X 10 at each continuity point of F X 10 . In the zoomed part of Fig. 4 (right), one can observe that Post-Widder's method seems a little less accurate than the proposed DRFI approximation. Figure 5 represents the corresponding approximations for the pdf of X 10 . Even if there is no theoretical convergence result for the DRFI approximation, the closeness of the new approach with the non-parametric one gives us good reason to believe that our approximation is efficient, keeping in mind, however, that the non-parametric estimation is obtained from approximate realizations of X t . Note that, here again, the zoomed part of Fig. 5 (right) shows that Post-Widder's method seems a little less accurate than the DRFI approximation.

Conclusion
Different tools based on some discretization of the rate function of an EGP have been presented here, firstly, for the approximate simulation of an EGP, and secondly, for the numerical assessment of the cdf/pdf of an EGP. Three discretization schemes have been proposed: one provides the best approximation among the three (b (ε) (t)) and the other two (b (ε,+) (t), b (ε,−) (t)) provide bounds for X t and for its cdf F X t .
As far as simulation procedures are concerned, it seems that our approximate simulation scheme behaves as well as two of the most usual ones (Bondesson's and rejection methods), previously developed in the context of subordinators. Also, our simulation scheme seems to behave better than the one proposed by Ishwaran and James (2004), specifically developed for the simulation of an EGP.
As for the numerical assessment of the cdf of an EGP, we could not find in the literature any available procedure with a similar control on the precision (except from Veillette and Taqqu (2011), which provides a refinement of the Post-Widder's method as well as some asymptotic bounds for the error when H → ∞). Beyond this control on the precision, the results provided by the proposed method have been compared to those obtained by Post-Widder's formula, at the advantage of our method. Even if we do not have any similar control on the precision for the pdf, a good behavior of the proposed approximation has been numerically observed, here again at its advantage when compared to Post-Widder.
To sum up, the discretized rate function method seems to behave well, both for simulating approximate paths and for the numerical assessment of the cdf/pdf of an EGP.
Before being able to use the proposed tools in an applied situation (for e.g. deterioration modelling in an industrial reliability context), another important issue that requires further study concerns the development of statistical estimation procedures for an EGP. In his seminal paper, Ç inlar (1980) proposes an iterative procedure that seems difficult to use in practice (since it needs notably a test for deciding whether a sample path comes from an ordinary Gamma process), whereas Guida et al. (2012), in a parametric context, apply an approximate maximum likelihood method. The study of estimation procedures for an EGP is in progress and will be the subject of a future paper.