# Approximate simulation techniques and distribution of an extended Gamma process

Zeina Al Masry, Sophie Mercier, Ghislain Verdier

► To cite this version:

Zeina Al Masry, Sophie Mercier, Ghislain Verdier. Approximate simulation techniques and distribution of an extended Gamma process. Methodology and Computing in Applied Probability, 2017, 19 (1), pp.213-235. 10.1007/s11009-015-9474-3. hal-01247062

HAL Id: hal-01247062

https://hal.science/hal-01247062

Submitted on 27 Dec 2015

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Approximate simulation techniques and distribution of an extended Gamma process

Zeina AL MASRY · Sophie MERCIER · Ghislain VERDIER

Received: date / Accepted: date

**Abstract** 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 (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.

**Keywords** Reliability · Degradation · Process with independent increments · Series representation · Inverse Laplace transform · Post-Widder formula · Weighted Gamma process

## 1 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 \ge 0}$ (say), where the system lifetime $\tau$ typically corresponds to the hitting time of a given failure threshold $L$:

$$ \tau = \inf(t > 0 : X_t > L). $$

<footnote>Laboratoire de Mathématiques et de leurs Applications - Pau (UMR CNRS 5142)</footnote>

<footnote>Université de Pau et des Pays de l'Adour, Bâtiment IPRA, Avenue de l'Université,</footnote>

<footnote>F-64013 PAU cedex, FRANCE</footnote>

<footnote>Tel.: +33 (0)5 59 40 75 37</footnote>

<footnote>Fax: +33 (0)5 59 40 75 55</footnote>

<footnote>E-mail: zeina.almasry@univ-pau.fr, sophie.mercier@univ-pau.fr,</footnote>

<footnote>ghislain.verdier@univ-pau.fr</footnote>
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$:

$$R(t) = \mathbb{P}(\tau > t) = \mathbb{P}(X_t \le L) = F_{X_t}(L) = \int_0^L f_{X_t}(s) ds,$$

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

$$MTBF = \mathbb{E}(\tau) = \int_{0}^{\infty} \mathbb{P}(\tau > t) dt = \int_{0}^{\infty} F_{X_t}(L) dt.$$

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 \ge 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 \ge 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.

## 2 Definition of an EGP and first properties

Let $a: \mathbb{R}_+ \to \mathbb{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 $\Gamma_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

$$f_t(x) = \frac{b_0^{a(t)}}{\Gamma(a(t))} x^{a(t)-1} \exp(-b_0 x), \forall x \ge 0,$$

e.g. see Abdel-Hameed (1975).

Now, let $b: \mathbb{R}_+ \to \mathbb{R}_+$ be a measurable positive function such that, for all $t > 0$:

$$\int_{(0,t]} \frac{da(s)}{b(s)} < \infty. \qquad (1)$$

Following (Çinlar 1980; Dykstra and Laud 1981), the process $X = (X_t)_{t \ge 0}$ is said to be an EGP with shape function $a(.)$ and scale function $b(.)$ (written $X \sim \Gamma(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 \ge 0} \sim \Gamma_0(a(t), 1)$:

$$X_t = \int_{(0,t]} \frac{dY_s}{b(s)}, \forall t > 0 \qquad (2)$$

and $X_0 = 0$.

If $b(.)$ is constant and equal to $b_0$, the EGP simply reduces to a standard Gamma process $\Gamma_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

$$ \mathcal{L}_{X_{t+h}-X_t}(\lambda) := \mathbb{E} \left( e^{-\lambda(X_{t+h}-X_t)} \right) = \exp \left( - \int_{(t,t+h]} \log \left( 1 + \frac{\lambda}{b(s)} \right) da(s) \right), \quad (3) $$

for all $t, \lambda \ge 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 $\sim \Gamma(a(t), b(t))$.

Finally, the mean and variance of an EGP are given by

$$ \mathbb{E}(X_t) = \int_{(0,t]} \frac{da(s)}{b(s)} \quad \text{and} \quad \mathbb{V}(X_t) = \int_{(0,t]} \frac{da(s)}{b(s)^2}. \quad (4) $$

**Remark 1** Let $X \sim \Gamma(a(t), b(t))$. So far, the shape function $a$ was supposed to be right-continuous. 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 \sim \Gamma(a^c(t), b(t))$ and $X^d \sim \Gamma(a^d(t), b(t))$. This second EGP simply reduces, for each $t \ge 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

$$ X_{\infty} = \lim_{t \to \infty} X_t (\le \infty), $$

which exists, due to the non decreasingness property of an EGP.

**Proposition 1** Let $X = (X_t)_{t \ge 0}$ be an EGP with $X \sim \Gamma(a(t), b(t))$. We then have:

$$ \int_{(0,\infty)} \log \left( 1 + \frac{1}{b(s)} \right) da(s) \le \mathbb{E}(X_{\infty}) = \int_{(0,\infty)} \frac{da(s)}{b(s)} \le \infty \quad (5) $$

and the following results:

1. If $\mathbb{E}(X_{\infty}) < \infty$, then $X_{\infty}$ is almost surely finite.

2. If

$$ \int_{(0,\infty)} \log \left( 1 + \frac{1}{b(s)} \right) da(s) = \infty, \quad (6) $$

then $X_{\infty}$ is almost surely infinite.

**Proof** Inequation (5) and the first point are clear. For the second one, let us note that

$$ \mathbb{E}(e^{-X_{\infty}}) = \mathcal{L}_{X_{\infty}}(1) = \exp \left( - \int_{(0,\infty)} \log \left( 1 + \frac{1}{b(s)} \right) da(s) \right). $$

Under assumption (6), we consequently have $\mathbb{E}(e^{-X_{\infty}}) = 0$. This implies that $e^{-X_{\infty}} = 0$ almost surely and the result.
Remark 2 In the case where

$$ \int_{(0,\infty)} \log \left( 1 + \frac{1}{b(s)} \right) da(s) < \infty = \int_{(0,\infty)} \frac{da(s)}{b(s)}, \quad (7) $$

the dichotomy $P(X_\infty = \infty) \in \{0, 1\}$ is not valid any more and both $P(X_\infty = \infty) > 0$ and $P(X_\infty < \infty) > 0$ are possible, as is illustrated in the following example (third case).

Example 1 The different types of behavior are illustrated in Figure 1, 2 and 3 (left) where a few sample paths are plotted for three different EGPs (using the discretized rate method from Section 5): $X^{(1)} \sim \Gamma(t, (t+1)^2)$, $X^{(2)} \sim \Gamma(t, \frac{1}{(t+1)^2})$ and $X^{(3)} \sim \Gamma(\frac{1}{2}(1 - \frac{1}{(t+1)^2}), \frac{1}{\exp(1+t)-1})$. 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 Figure 1 are stabilizing whereas all trajectories of Figure 2 are exploding. In Figure 3, some trajectories are exploding and others are stabilizing. For a better insight into their behavior, the variance-to-mean ratios of the $X^{(i)}$'s are also plotted in the same figures (right). These figures illustrate the flexibility of EGPs for modeling purpose.

In all the following, $X = (X_t)_{t \ge 0}$ stands for an EGP $\sim \Gamma(a(t), b(t))$ and $Y = (Y_t)_{t \ge 0}$ stands for a standard Gamma process $\Gamma_0(a(t), 1)$, without any further notification. We recall that the shape function $a(t)$ is assumed to be continuous.

## 3 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.

### 3.1 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 \in [0, T]$, where $T > 0$. More specifically, for each $t \in [0; T]$ and each $K^* \in \mathbb{N}$, they consider:

$$ G_t^{(K^*)} = \sum_{k=1}^{K^*} 1_{[0,t]}(V_k) \frac{1}{b(V_k)} W_k^{(K^*)}, \quad (8) $$

where $(V_k)_{k=1,...,K^*}$ are i.i.d random variables (r.v.s) with distribution $\frac{da(v)}{a(T)} 1_{[0,T]}(v)$ and $(W_k^{(K^*)})_{k=1,...,K^*}$ are i.i.d r.v.s with Gamma distribution $\Gamma_0(\frac{a(T)}{K^*}, 1)$, independent from the $V_k$'s.

Ishwaran and James (2004) proved that $(G_t^{(K^*)})_{t \ge 0} \to (X_t)_{t \ge 0}$ weakly when $K^* \to \infty$ in the Skorokhod space $\mathcal{D}(0, \infty)$ of right-continuous functions with left-side limits. (To be
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Sample paths (left) and ratio $\frac{V(X_t^{(1)})}{E(X_t^{(1)})}$ (right) as a function of time</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2** Sample paths (left) and ratio $\frac{V(X_t^{(2)})}{E(X_t^{(2)})}$ (right) as a function of time</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Sample paths (left) and ratio $\frac{V(X_t^{(3)})}{E(X_t^{(3)})}$ (right) as a function of time</figcaption></figure>

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 and Vere-Jones 2007, Lemma 11.1.XI)).
One can easily check that

$$
\begin{align*}
\mathbb{E}[G_t^{(K^*)}] &= \mathbb{E}[X_t], \\
\mathbb{V}[G_t^{(K^*)}] &= \mathbb{V}[X_t] + \frac{1}{K^*} (a(T)\mathbb{V}[X_t] - (\mathbb{E}[X_t])^2)
\end{align*}
$$

so that the approximation of $X_t$ by $G_t^{(K^*)}$ is unbiased. Also, the rate of convergence of the variance is $O(1/K^*)$.

3.2 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.

**Proposition 2** Let $T > 0$ and let $(U_n)_{n \ge 1}$ be the points of a homogeneous Poisson process $M$ with parameter $a(T)$. Let also $(V_n)_{n \ge 1}$ be a sequence of i.i.d. r.v.s with distribution $H(dw) = \frac{da(v)}{a(T)}1_{[0,T]}(v)$, independent of $M$ and let $(W_n)_{n \ge 1}$ be a sequence of i.i.d r.v.s with distribution $P_W(dw)$, independent of $M$ and of the $V_n$'s. Then we have the four fol- lowing series representations of an EGP, where $\mathcal{D}$ means "is identically distributed as":

(1) Bondesson's series representation (Bond):

$$
X_t \stackrel{\mathcal{D}}{=} \sum_{n \ge 1} \frac{1}{b(V_n)} \exp(-U_n) W_n 1_{[0,t]}(V_n), \text{ for } 0 \le t \le T, \quad (9)
$$

where {$W_n$}$_{n \ge 1}$ is a sequence of i.i.d exponential r.v.s with mean 1,

(2) Rejection's series representation (Rej):

$$
X_t \stackrel{\mathcal{D}}{=} \sum_{n \ge 1} \frac{1}{b(V_n) \exp(U_n) - 1} \mathbb{1}_E (U_n, W_n) \mathbb{1}_{[0,t]} (V_n), \text{ for } 0 \le t \le T, \quad (10)
$$

where $E = \left\{ (u, w) \in \mathbb{R}_+ \times [0, 1] \mid \frac{\exp(u)}{\exp(u)-1} \exp\left(-\left(\exp(u)-1\right)^{-1}\right) \ge w \right\}$ and where
$\{W_n\}_{n \ge 1}$ is a sequence of i.i.d uniform r.v.s on $[0, 1]$,

(3) Thinning's series representation (Thin):

$$
X_t \stackrel{\mathcal{D}}{=} \sum_{n \ge 1} \frac{1}{b(V_n)} W_n 1_{(W_n U_n \le 1)} 1_{[0,t]}(V_n), \text{ for } 0 \le t \le T, \quad (11)
$$

where {$W_n$}$_{n \ge 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):

$$
X_t \stackrel{\mathcal{D}}{=} \sum_{n \ge 1} \frac{1}{b(V_n)} E_1^{-1}(U_n) 1_{[0,t]}(V_n), \text{ for } 0 \le t \le T, \quad (12)
$$

with $E_1(x) = \int_x^\infty u^{-1} \exp(-u) du$, the exponential integral function.
Proof The four series representations can be written in the same manner:

$$X_t = \sum_{n \ge 1} f_t(U_n, V_n, W_n) = \sum_{n \ge 1} \frac{1}{b(V_n)} H(U_n, W_n) \mathbb{1}_{[0,t]}(V_n), \quad (13)$$

with $f_t(u, v, w) = \frac{1}{b(v)} H(u, w) \mathbb{1}_{[0,t]}(v)$ and an appropriate choice for $H : \mathbb{R}_+^2 \to \mathbb{R}_+$. As already mentioned in Section 2, $(X_t)_{0 \le t \le 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 \ge 1}$ can be seen as the points of a Poisson random measure $N$ on $\mathbb{R}_+^3$ with intensity

$$\mu(du, dv, dw) = du \times da(v) \mathbb{1}_{[0,T]}(v) \times P_W(dw)$$

(see for example (Çınlar 2011, Corollary 3.5, p. 265)).

Considering $t_1 < \dots < t_n$, we then have $X_{t_{i+1}} - X_{t_i} = N(f_{t_{i+1}} - f_{t_i})$ for all $1 \le i \le n-1$, where $f_{t_{i+1}} - f_{t_i}$ takes range in $\mathbb{R}_+^3 \times (t_i, t_{i+1}] \times \mathbb{R}_+^3$. As the ranges of the $f_{t_{i+1}} - f_{t_i}$ are non overlapping, the increments $X_{t_{i+1}} - X_{t_i}$ are independent.

Now, based on the Laplace functional theorem for a Poisson random measure (Çınlar 2011, Theorem 2.9, p. 252) for the third line, we have for all $t \in [0, T]$ and all $h, \lambda > 0$:

$$\begin{aligned}
\mathcal{L}_{X_{t+h} - X_t}(\lambda) \\
&= \mathbb{E}[\exp(-\lambda(X_{t+h} - X_t))] \\
&= \mathbb{E}[\exp(-N(\lambda(f_{t+h} - f_t)))] \\
&= \exp\left(-\int_{\mathbb{R}^3} (1 - e^{-(\lambda(f_{t+h} - f_t))}) d\mu\right) \\
&= \exp\left(-\int_{(t, t+h]} Q(\lambda, v) da(v)\right)
\end{aligned} \quad (14)$$

where

$$Q(\lambda, v) = \int_{\mathbb{R}^2} \left( 1 - \exp \left( -\lambda \frac{1}{b(v)} H(u, w) \right) \right) du \times P_W(dw).$$

From (3) and (14), it then suffices to check that $Q(\lambda, v) = \log \left(1 + \frac{\lambda}{b(v)}\right)$ for each of the four series representations.

1) Bondesson:

In this case, $H(u, w) = \exp(-u)w$ and $P_W(dw) = \exp(-w)dw$. Then

$$\begin{aligned}
Q(\lambda, v) &= \int_0^\infty \left( \int_0^\infty \left( 1 - \exp\left(-\frac{\lambda}{b(v)} \exp(-u)w\right) \right) \exp(-w) dw \right) du \\
&= \int_0^\infty \left( 1 - \frac{b(v)}{b(v) + \lambda \exp(-u)} \right) du \\
&= \int_0^\infty \left( \frac{\lambda \exp(-u)}{b(v) + \lambda \exp(-u)} \right) du \\
&= \log \left( 1 + \frac{\lambda}{b(v)} \right).
\end{aligned}$$
2) Rejection:

Here, $H(u, w) = \frac{1}{\exp(u)-1} \mathbb{1}_{(\frac{\exp(u)}{\exp(u)-1} \exp(-\frac{1}{\exp(u)-1}) \ge w)}$ and $P_W(dw) = dw$. Accordingly,

$$ 
\begin{aligned}
Q(\lambda, v) &= \int_0^\infty \left( \int_0^\frac{\exp(u)}{\exp(u)-1} \exp\left(-\frac{1}{\exp(u)-1}\right) \left(1 - \exp\left(-\frac{\lambda}{b(v)}\frac{1}{\exp(u)-1}\right)\right) dw \right) du \\
&= \int_0^\infty \left(1 - \exp\left(-\frac{\lambda}{b(v)}\frac{1}{\exp(u)-1}\right)\right) \frac{\exp(u)}{\exp(u)-1} \exp\left(-\frac{1}{\exp(u)-1}\right) du \\
&= \int_0^\infty \left(1 - \exp\left(-\frac{\lambda}{b(v)}z\right)\right) \frac{\exp(-z)}{z} dz
\end{aligned} 
\quad (15) $$

setting $z = \frac{1}{\exp(u)-1}$ in the last line.

(15) is a Frullani's integral, which can classically be computed by differentiating with respect of $\lambda$: Setting $o(\lambda, v, z) = \left(1 - \exp\left(-\frac{\lambda}{b(v)}z\right)\right) \frac{\exp(-z)}{z}$, we get by dominated convergence that

$$ \frac{\partial Q(\lambda, v)}{\partial \lambda} = \int_0^\infty \frac{\partial o(\lambda, v, z)}{\partial \lambda} dz = \frac{1}{b(v) + \lambda}, $$

which provides the result.

3) Thinning:

$H(u, w) = w \mathbb{1}_{(uw \le 1)}$ and $P_W(dw) = \exp(-w) dw$.

$$ 
\begin{aligned}
Q(\lambda, v) &= \int_0^\infty \left( \int_0^\frac{1}{w} \left(1 - \exp\left(-\frac{\lambda}{b(v)}w\right)\right) \exp(-w) du \right) dw \\
&= \int_0^\infty \left(1 - \exp\left(-\frac{\lambda}{b(v)}w\right)\right) \frac{\exp(-w)}{w} dw.
\end{aligned} 
$$

We recognize (15), which allows to conclude.

4) Inversion of Lévy measure:

$H(u, w) = E_1^{-1}(u)$ is independent on $w$. Then

$$ 
\begin{aligned}
Q(\lambda, v) &= \int_0^\infty \left(1 - \exp\left(-\frac{\lambda}{b(v)}E_1^{-1}(u)\right)\right) du \\
&= \int_0^\infty \left(1 - \exp\left(-\frac{\lambda}{b(v)}z\right)\right) \frac{\exp(-z)}{z} dz
\end{aligned} 
$$

setting $z = E_1^{-1}(u)$. We recognize (15) again, which allows to conclude. ■

Approximate simulation of $(X_t)_{0 \le t \le 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$:

$$ X_t^{(B)} = \sum_{n \ge 1} \frac{1}{b(V_n)} H(U_n, W_n) \mathbb{1}_{[0,t]}(V_n) \mathbb{1}_{[0,B]}(U_n) $$
for all $t \in [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

$$ \tilde{X}_t^{(B)} = X_t - X_t^{(B)} = \sum_{n \ge 1} \frac{1}{b(V_n)} H(U_n, W_n) \mathbb{1}_{[0,t]}(V_n) \mathbb{1}_{(B,\infty)}(U_n) \quad (16) $$

for all $t \in [0, T]$, to be the remainder of the truncated series. By construction, both $(\tilde{X}_t^{(B)})_{t \in [0,T]}$ and $(X_t^{(B)})_{t \in [0,T]}$ have independent increments. The moments of $\tilde{X}_t^{(B)}$ are given in the following proposition.

**Proposition 3** Recall that $\tilde{X}_t^{(B)}$ is the remainder of the truncated series (16). We have

(1) Bondesson's series representation:

$$ \tilde{X}_t^{(B)} \sim \Gamma \left( a(t), \frac{b(t)}{\exp(-B)} \right); $$

$$ \mathbb{E}[\tilde{X}_t^{(B)}] = \exp(-B) \mathbb{E}[X_t]; $$

$$ \mathbb{V}[\tilde{X}_t^{(B)}] = \exp(-2B) \mathbb{V}[X_t], $$

(2) Rejection's series representation:

$$ \mathbb{E}[\tilde{X}_t^{(B)}] = E[X_t] \left[ 1 - \exp \left( - \frac{1}{\exp(B) - 1} \right) \right]; $$

$$ \mathbb{V}[\tilde{X}_t^{(B)}] = V[X_t] \left[ 1 - \left( 1 + \frac{1}{\exp(B) - 1} \right) \exp \left( - \frac{1}{\exp(B) - 1} \right) \right], $$

(3) Thinning's series representation:

$$ \mathbb{E}[\tilde{X}_t^{(B)}] = \mathbb{E}[X_t] \left[ 1 - B + B \exp \left( - \frac{1}{B} \right) \right]; $$

$$ \mathbb{V}[\tilde{X}_t^{(B)}] = \mathbb{V}[X_t] \left[ 1 + \exp \left( - \frac{1}{B} \right) + 2B \exp \left( - \frac{1}{B} \right) - 2B \right], $$

(4) Inversion of Lévy measure's series representation:

$$ \mathbb{E}[\tilde{X}_t^{(B)}] = \mathbb{E}[X_t](1 - E_1(B)); $$

$$ \mathbb{V}[\tilde{X}_t^{(B)}] = \mathbb{V}[X_t](1 - \exp(-E_1(B)) - E_1(B) \exp(-E_1(B))). $$

**Proof** Starting from (16), the Laplace transform of $\tilde{X}_t^{(B)}$ is computed using the Laplace functional theorem for Poisson random measures. In the Bondesson case, we recognize the Laplace transform of an EGP $\sim \Gamma(a(t), \frac{b(t)}{\exp(-B)})$. For the three other representations, the mean and variance are obtained through the first and second order derivatives of the Laplace transform at $\lambda = 0$. ■

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 \in [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.
<table><caption>**Table 1** Asymptotic equivalents of $E(\tilde{X}_t^{(B)})/E(X_t)$ and $V(\tilde{X}_t^{(B)})/V(X_t)$ when $B \to \infty$</caption><thead><tr><th>Method</th><th>E(Χ̃<sub>t</sub><sup>(B)</sup>)/E(X<sub>t</sub>)</th><th>V(Χ̃<sub>t</sub><sup>(B)</sup>)/V(X<sub>t</sub>)</th></tr></thead><tbody><tr><td>Bondesson</td><td>exp(-B)</td><td>exp(-2B)</td></tr><tr><td>Rejection</td><td>exp(-B)</td><td>exp(-2B)/2</td></tr><tr><td>Thinning</td><td>1/2B</td><td>1/6B<sup>2</sup></td></tr><tr><td>Inversion of Lévy measure</td><td>E<sub>1</sub>(B)</td><td>E<sub>1</sub>(B)<sup>2</sup>/2</td></tr></tbody></table>

From Table 1, we can note that the Inversion of Lévy measure's method converges very fast ($E_1(B) = o(\exp(-B))$ when $B \to \infty$). 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.

## 4 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 \ge 0$,

$$F_{X_t}(x) = \lim_{H \to \infty} \sum_{k=0}^{H} \left(-\frac{H}{x}\right)^k \frac{1}{k!} \mathcal{L}_{X_t}^{(k)}\left(\frac{H}{x}\right), \quad (17)$$

$$f_{X_t}(x) = \lim_{H \to \infty} -\frac{1}{H!} \left(-\frac{H}{x}\right)^{H+1} \mathcal{L}_{X_t}^{(H)}\left(\frac{H}{x}\right), \quad (18)$$

where $\mathcal{L}_{X_t}^{(k)}$ denotes the k-th derivative of the Laplace transform $\mathcal{L}_{X_t}$ for $k \ge 1$ and $\mathcal{L}_{X_t}^{(0)} = \mathcal{L}_{X_t}$.

Let $\phi_t$ be the Laplace exponent of $X_t$, with

$$\phi_t(\lambda) = -\ln(\mathcal{L}_{X_t}(\lambda)), \forall \lambda \ge 0.$$

Following Veillette and Taqqu (2011), Leibniz's formula provides:

$$\mathcal{L}_{X_t}^{(k+1)}(\lambda) = -\sum_{j=0}^{k} \binom{k}{j} \mathcal{L}_{X_t}^{(j)}(\lambda) \phi_t^{(k+1-j)}(\lambda), \forall \lambda, t, k \ge 0, \quad (19)$$

with $\phi_t^{(j)}$ the j-th derivative of $\phi_t$. An expression of $\phi_t^{(j)}$ is given in the following Lemma in the case of an EGP.
**Lemma 1** Let $\lambda \ge 0, j \ge 1$ and let $X \sim \Gamma(a(t), b(t))$.

Then

$$ \phi_t^{(j)}(\lambda) = (-1)^{j-1}(j-1)! m_t^{(j)}(\lambda) $$

where

$$ m_t^{(j)}(\lambda) = \int_0^t \frac{da(s)}{(b(s) + \lambda)^j} $$

*Proof* Based on (3) and (15) and setting $u = z/b(s)$ in the third line we have:

$$ \begin{aligned} \phi_t(\lambda) &= \int_{(0,t]} \log \left(1 + \frac{\lambda}{b(s)}\right) da(s) \\ &= \int_{(0,t]} \left( \int_0^\infty \left(1 - \exp\left(-\frac{\lambda}{b(s)}z\right)\right) \frac{\exp(-z)}{z} dz \right) da(s) \\ &= \int_0^\infty \left(1 - e^{-\lambda u}\right) \mu_t(du) \end{aligned} $$

where

$$ \mu_t(du) = \frac{1}{u} \left( \int_{(0,t]} e^{-b(s)u} da(s) \right) du $$

is a measure on $\mathbb{R}_+$ such that $\int_{\mathbb{R}_+} (1 \wedge x) \mu_t(dx) < \infty$. From Veillette and Taqqu (2011), $\phi_t^{(j)}(\lambda)$ is directly calculated by

$$ \phi_t^{(j)}(\lambda) = (-1)^{j+1} \int_0^\infty x^j e^{-\lambda x} \mu_t(dx) \quad (20) $$

for all $j \ge 1$. This provides:

$$ \begin{aligned} \phi_t^{(j)}(\lambda) &= (-1)^{j+1} \int_0^t \left( \int_0^\infty x^{j-1} e^{-(b(s)+\lambda)x} dx \right) da(s) \\ &= (-1)^{j+1} \int_0^t \frac{\Gamma(j)}{(b(s)+\lambda)^j} da(s), \end{aligned} $$

renormalizing the pdf of $\Gamma_0(j, b(s) + \lambda)$ to get the last line. ■

Based on (19) and Lemma 1, $\mathcal{L}_{X_t}^{(k)}(\lambda)$ can be recursively computed through

$$ \mathcal{L}_{X_t}^{(k+1)}(\lambda) = -k! \sum_{j=0}^{k} \mathcal{L}_{X_t}^{(j)}(\lambda) \frac{(-1)^{k-j}}{j!} m_t^{(k+1-j)}(\lambda). $$

This allows to compute an approximation of the pdf and cdf of an EGP through (17 – 18), taking $H$ large enough.

## 5 Discretized rate function method

In case of a piecewise constant scale function $b(\cdot)$, the process $(X_t)_{t \ge 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.
## 5.1 Construction of the approximate process $X^{(\varepsilon)}$

Let $T > 0$. The following assumption is considered:

$$b(.) \text{ is continuous on } (0, T] \text{ and } \exists m > 0 \text{ such that } \forall t \in (0, T], b(t) \ge m. \quad (21)$$

Let $\varepsilon > 0$. The following piecewise constant approximation of $\frac{1}{b(.)}$ on $(0, T]$, denoted by $\frac{1}{b^{(\varepsilon)}(.)}$, is considered:

$$\forall t \in (0, T], \quad \frac{1}{b^{(\varepsilon)}(t)} = \sum_{i=0}^{n(\varepsilon)} \frac{1}{b_i} \mathbb{1}_{[l_i, l_{i+1})}(t) \quad (22)$$

where $n(\varepsilon)$ is such that $l_{n(\varepsilon)} \le T < l_{n(\varepsilon)+1}$, $l_0 = 0$, and the $l_i$'s, for $i = 1, \dots, n(\varepsilon) + 1$, are recursively defined by:

$$l_{i+1} = \sup \left\{ l \in (l_i, T] : \forall l' \in [l_i, l], \left| \frac{1}{b(l_i)} - \frac{1}{b(l')} \right| < \varepsilon \right\}. \quad (23)$$

Note that assumption (21) insures that $n(\varepsilon)$ is finite. The constants $b_i$, $i = 0, \dots, n(\varepsilon)$ are next defined by:

$$\forall i = 0, \dots, n(\varepsilon) - 1, \quad \frac{1}{b_i} = \frac{1}{l_{i+1} - l_i} \int_{l_i}^{l_{i+1}} \frac{1}{b(s)} ds, \quad (24)$$

and

$$\frac{1}{b_{n(\varepsilon)}} = \frac{1}{T - l_{n(\varepsilon)}} \int_{l_{n(\varepsilon)}}^{T} \frac{1}{b(s)} ds. \quad (25)$$

Now that $b^{(\varepsilon)}(.)$ is constructed, we set

$$X_t^{(\varepsilon)} = \int_{(0,t]} \frac{dY_s}{b^{(\varepsilon)}(s)}, \quad \forall t \in [0, T], \quad (26)$$

where $(Y_t)_{t \ge 0}$ is the same standard Gamma process $\Gamma_0(a(t), 1)$ as in Equation (2) defining $X$. Then $(X_t^{(\varepsilon)})_{t \in [0, T]}$ is the restriction to $[0, T]$ of an EGP $\sim \Gamma(a(t), b^{(\varepsilon)}(t))$ with the same shape function $a(.)$ as $X$ and with $b^{(\varepsilon)}(.)$ as scale function.

**Remark 3** Based of the definition (2) of an EGP, we prefer to construct an approximation of $\frac{1}{b(.)}$ instead of $b(.)$. Also, by construction, for all $t$ in $(0, T]$,

$$\left| \frac{1}{b(t)} - \frac{1}{b^{(\varepsilon)}(t)} \right| \le \varepsilon. \quad (27)$$
## 5.2 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^{(\epsilon)}$, for $t \in [0, T]$.

**Proposition 4** Let

$$ \tilde{X}_t^{(\epsilon)} = X_t - X_t^{(\epsilon)} = \int_{(0,t]} \left( \frac{1}{b(s)} - \frac{1}{b^{(\epsilon)}(s)} \right) dY_s, \quad \forall t \in [0, T]. \quad (28) $$

Then:

$$ \mathbb{E}[|\tilde{X}_t^{(\epsilon)}|] \le \epsilon a(t), \quad (29) $$

$$ \mathbb{V}[\tilde{X}_t^{(\epsilon)}] \le \epsilon^2 a(t) \quad (30) $$

for all $t \in [0, T]$.

**Proof** Setting $\tilde{b}^{(\epsilon)}(s) = (1/b(s) - 1/b^{(\epsilon)}(s))^{-1}$ if $b(s) \ne b^{(\epsilon)}(s)$ and $b^{(\epsilon)}(s) = \infty$ otherwise, we have

$$ \tilde{X}_t^{(\epsilon)} = \int_0^t \frac{1}{\tilde{b}^{(\epsilon)}(s)} dY_s $$

where $\frac{1}{|\tilde{b}^{(\epsilon)}(s)|} \le \epsilon$, based on (27).

We easily derive that

$$ |\tilde{X}_t^{(\epsilon)}| \le \int_0^t \frac{1}{|\tilde{b}^{(\epsilon)}(s)|} dY_s \le \epsilon Y_t \quad (31) $$

from where we derive (29).

Also, one can check that (4) is still valid for $\tilde{X}_t^{(\epsilon)}$ with $b(t)$ substituted by $\tilde{b}^{(\epsilon)}$, so that

$$ \mathbb{V}[\tilde{X}_t^{(\epsilon)}] = \int_{(0,t]} \frac{1}{(\tilde{b}^{(\epsilon)}(s))^2} da(s) \le \epsilon^2 a(t). $$

<figure><img src="image_5.png" /></figure>

We now provide some convergence results.

**Proposition 5** Let $X \sim \Gamma(a(t), b(t))$, $(\epsilon_n)_{n \ge 1}$ be a sequence of positive real numbers and let $(X^{(\epsilon_n)})_{n \ge 1}$ be a sequence of EGPs with shape function $a(\cdot)$ and scale function $b^{(\epsilon_n)}(\cdot)$. Suppose that assumption (21) is satisfied.

1) If $\sum_{n \ge 1} \sqrt{\epsilon_n} < \infty$, then $\sup_{t \in [0, T]} |X_t^{(\epsilon_n)} - X_t| \xrightarrow{n \to \infty} 0$ almost surely.

2) If $\epsilon_n \xrightarrow{n \to \infty} 0$, then $\sup_{t \in [0, T]} |X_t^{(\epsilon_n)} - X_t| \xrightarrow{n \to \infty} 0$ in quadratic mean.
Proof 1) According to (31), we have:

$$
\begin{align*}
\mathbb{P} \left( \sup_{t \in [0,T]} |X_t^{(\epsilon_n)} - X_t| > \sqrt{\epsilon_n} \right) &\le \mathbb{P} \left( \epsilon_n \sup_{t \in [0,T]} Y_t > \sqrt{\epsilon_n} \right) \\
&= \mathbb{P} \left( Y_T > \frac{1}{\sqrt{\epsilon_n}} \right) \\
&\le \mathbb{E} (Y_T) \sqrt{\epsilon_n},
\end{align*}
$$

based on Markov inequality for the last line. We derive that

$$
\sum_{n \ge 1} \mathbb{P} \left( \sup_{t \in [0,T]} |X_t^{(\epsilon_n)} - X_t| > \sqrt{\epsilon_n} \right) < \infty,
$$

which concludes the first point, based on e.g. (Çınlar 2011, Proposition 2.7, p. 98).

2) In the same way, we have:

$$
\mathbb{E} \left[ \left( \sup_{t \in [0,T]} |X_t^{(\epsilon_n)} - X_t| \right)^2 \right] \le \mathbb{E} \left[ \left( \epsilon_n \sup_{t \in [0,T]} Y_t \right)^2 \right] = \epsilon_n^2 \mathbb{E} [Y_T^2],
$$

which allows to conclude. ■

Remark 4 Note that a consequence of point 2 is that, under the same assumptions,
$(X_t^{(\epsilon)})_{t \in [0,T]} \longrightarrow (X_t)_{t \in [0,T]}$ weakly in the Skorokhod space $\mathcal{D}([0,T])$.

5.3 Approximation of the pdf and cdf of $X_t^{(\epsilon)}$ (and of $X_t$)

As previously mentioned, $X^{(\epsilon)}$ is an EGP with a piecewise constant scale function. For each $t$, $X_t^{(\epsilon)}$ can be written as a sum of independent Gamma distributed random variables with different scale parameters. Several methods are available in the literature to evaluate the pdf of such a sum, see Nadarajah (2008). In Guida et al (2012), the authors use the method from Moschopoulos (1985) to evaluate the pdf of their discrete time EGP. After performing some numerical comparisons (not provided here), we finally chose to use the method from (Peppas 2011, Sub. 4.2), which provides an expression of both pdf and cdf of $X_t^{(\epsilon)}$ in terms of an infinite integral. For all $t$ in $[0, T]$ and all $x \ge 0$, this writes:

$$
f_{X_t^{(\epsilon)}}(x) = \lim_{K \to \infty} f_{X_t^{(\epsilon)}}^{(K)}(x) \text{ and } F_{X_t^{(\epsilon)}}(x) = \lim_{K \to \infty} F_{X_t^{(\epsilon)}}^{(K)}(x),
$$

where

$$
f_{X_t^{(\epsilon)}}^{(K)}(x) = \frac{1}{\pi} \int_0^K \frac{\cos\left(\sum_{p=0}^P \alpha_p \arctan(u/b_p) - xu\right)}{\prod_{p=0}^P (1 + (u/b_p)^2)^{\alpha_p/2}} du, \quad (32)
$$

$$
F_{X_t^{(\epsilon)}}^{(K)}(x) = \frac{1}{2} - \frac{1}{\pi} \int_0^K \frac{\sin\left(\sum_{p=0}^P \alpha_p \arctan(u/b_p) - xu\right)}{u \prod_{p=0}^P (1 + (u/b_p)^2)^{\alpha_p/2}} du \quad (33)
$$

with $P = P(\epsilon, t)$ the single integer such that $l_P < t \le l_{P+1}$, $\alpha_p = a(l_{p+1}) - a(l_p)$ for $p = 0, \dots, P-1$ and $\alpha_P = a(t) - a(l_P)$.
Taking $K$ large enough, $(32 - 33)$ provide approximations for the pdf and cdf of $X_t^{(\epsilon)}$, and consequently of $X_t$.

Using that $|\sin(\cdot)| \le 1$ and that $(1 + (u/b_p)^2)^{\alpha_p/2} \ge (u/b_p)^{\alpha_p}$, we easily get an upper bound for the approximation error, with:

$$ 
\begin{aligned}
E_{X_t^{(\epsilon)}}^{(K)}(x) &= |F_{X_t^{(\epsilon)}}(x) - F_{X_t^{(\epsilon)}}(x)| \\
&\le \frac{1}{\pi} \int_K \int_0^\infty \frac{1}{u \prod_{p=0}^P (u/b_p)^{\alpha_p}} du \\
&\le \frac{\prod_{p=0}^P b_p^{\alpha_p}}{\pi} \int_K \int_0^\infty \frac{1}{u^{a(t)+1}} du \\
&\le \frac{\prod_{p=0}^P b_p^{\alpha_p}}{\pi a(t) K^{a(t)}}
\end{aligned} 
\quad (34) $$

for all $K \in \mathbb{N}$, due to $\sum_{p=0}^P \alpha_p = a(t)$ for the second-to-last line.

## 5.4 Bounds for the cdf of $X_t$

Consider two other piecewise constant approximations of $\frac{1}{b(t)}$ defined by

$$ \frac{1}{b^{(\epsilon,-)}(t)} = \sum_{i=0}^{n(\epsilon)} \frac{1}{b_i^{-}} 1_{[l_i, l_{i+1})}(t) \quad \text{and} \quad \frac{1}{b^{(\epsilon,+)}(t)} = \sum_{i=0}^{n(\epsilon)} \frac{1}{b_i^{+}} 1_{[l_i, l_{i+1})}(t) $$

for all $t \in (0, T]$, where

$$ \frac{1}{b_i^{-}} = \inf \left\{ \frac{1}{b(l)}, l \in [l_i, \min\{l_{i+1}, T\}) \right\} \quad \text{and} \quad \frac{1}{b_i^{+}} = \sup \left\{ \frac{1}{b(l)}, l \in [l_i, \min\{l_{i+1}, T\}) \right\}. $$

The following bounds are obtained:

$$ \frac{1}{b^{(\epsilon,-)}(t)} \le \frac{1}{b(t)} \le \frac{1}{b^{(\epsilon,+)}(t)}, \quad \text{for all } t \le T. \quad (35) $$

The induced EGPs are denoted by $X^{(\epsilon,-)}$ and $X^{(\epsilon,+)}$, respectively. They satisfy

$$ X_t^{(\epsilon,-)} \le X_t \le X_t^{(\epsilon,+)} \quad \text{for all } t \le T $$

so that

$$ F_{X_t^{(\epsilon,+)}(x)} \le F_{X_t(x)} \le F_{X_t^{(\epsilon,-)}(x)}, \forall t \le T, \forall x \ge 0. \quad (36) $$

**Theorem 1** Using the notations of Subsection 5.3 substituting $X^{(\epsilon)}$ by $X^{(\epsilon,-)}$ and by $X^{(\epsilon,+)}$, we set:

$$ m_t(x, \epsilon, K) = F_{X_t^{(\epsilon,+)}(x)}^{(K)} - E_{X_t^{(\epsilon,+)}(x)}^{(K)}, $$

$$ M_t(x, \epsilon, K) = F_{X_t^{(\epsilon,-)}(x)}^{(K)} + E_{X_t^{(\epsilon,-)}(x)}^{(K)} $$
for all $K \in \mathbb{N}^*$ and all $\varepsilon > 0$. Then,

$$m_t(x, \varepsilon, K) \le F_{X_t}(x) \le M_t(x, \varepsilon, K) \quad (37)$$

and

$$\left| F_{X_t}(x) - \frac{m_t(x, \varepsilon, K) + M_t(x, \varepsilon, K)}{2} \right| \le \frac{M_t(x, \varepsilon, K) - m_t(x, \varepsilon, K)}{2}. \quad (38)$$

*Proof* The bounds (37) are easily obtained starting from Equation (36) and using that

$$F_{X_t^{(\varepsilon,-)}}(x) \le F_{X_t^{(\varepsilon,-)}}^{(K)}(x) + E_{X_t^{(\varepsilon,-)}}^{(K)}(x) = M_t(x, \varepsilon, K)$$

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 $\varepsilon$ 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^{(\varepsilon,-)}(\cdot)$ and $b^{(\varepsilon,+)}(\cdot)$ satisfy (27). The theoretical results of Subsection 5.2 are then still valid for $X^{(\varepsilon,-)}$ and $X^{(\varepsilon,+)}$.

Finally, remark that for any fixed $s > 0$, the increment process $(X_{t+s} - X_s)_{t \ge 0}$ still is an EGP $\sim \Gamma'(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 \in (s, s + t]$ in the series representations from Proposition 2).

## 5.5 When assumption (21) is not satisfied

Assumption (21) allows to construct the piecewise constant function $b^{(\varepsilon)}(\cdot)$ specified in Subsection 5.1. When this assumption is not satisfied, alternatives for the construction of $b^{(\varepsilon)}(\cdot)$ have to be considered.

* An interesting case is when function $b(\cdot)$ satisfies $\lim_{t \to 0^+} b(t) = 0$. If the following assumption is checked,

$$\exists \eta > 0 \text{ such as } \int_0^\eta \frac{ds}{b(s)} < \infty,$$

then the piecewise constant approximation can be constructed from $l_0 = 0, l_1$ such that:

$$\int_0^{l_1} \frac{ds}{b(s)} \le \varepsilon \quad \text{and} \quad \int_0^{l_1} \frac{da(s)}{b(s)} \le \varepsilon,$$

and the $l_i$'s, for $i = 2, \dots, n(\varepsilon) + 1$, recursively defined according to (23). Note that in this case, inequality (27) is not satisfied any more.

* If function $b(\cdot)$ 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.
# 6 Numerical Experiments

## 6.1 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_{X_t(\varepsilon)}^{(K)}$ as provided by (33), where we recall that $X_t^{(\varepsilon)}$ is defined by (26). There hence are two parameters to define $(\varepsilon$ and $K)$.

All these approaches are investigated in this section. We set $X \sim \Gamma(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 \in [0,T]}$. Note that for Method 2, only Bondesson's series representation algorithm is given.

**Algorithm 1** Method 1: Ishwaran and James's approach

1. Fix $K^*$.

For $i$ in 1: $N$, repeat the following steps:

2. Simulate $K^*$ realizations $Z_1, \dots, Z_{K^*}$ according to the uniform distribution on $[0, 1]$.

3. Compute $V_j = a^{-1}(Z_j a(T))$ for $j = 1, \dots, K^*$ (inverse transform sampling).

4. Simulate $K^*$ realizations $W_1^{(K^*)}, \dots, W_{K^*}^{(K^*)}$ according to the Gamma distribution $\Gamma_0 \left( \frac{a(T)}{K^*}, 1 \right)$.

5. Compute $G_{t,i}^{(K^*)} = \sum_{j=1}^{K^*} 1_{[0,t]}(V_j) \frac{1}{b(V_j)} W_j^{(K^*)}$.

**Algorithm 2** Method 2: Bondesson's series representation

1. Fix $B$.

For $i$ in 1: $N$, repeat the following steps:

2. Simulate the total number of jumps $N_B$ according to the Poisson distribution $P(a(T)B)$.

3. Simulate $N_B$ realizations $U_1, \dots, U_{N_B}$ according to the uniform distribution on $[0, B]$.

4. Simulate $N_B$ realizations $W_1, \dots, W_{N_B}$ according to the exponential distribution with parameter 1.

5. Simulate $N_B$ realizations $Z_1, \dots, Z_{N_B}$ according to the uniform distribution on $[0, 1]$.

6. Compute $V_j = a^{-1}(Z_j a(T))$ for $j = 1, \dots, N_B$.

7. Compute $X_{t,i}^{(B)} = \sum_{j=1}^{N_B} \frac{1}{b(V_j)} \exp(-U_j) W_j 1_{[0,t]}(V_j)$.
**Algorithm 3 Method 3: Discretized rate method**

1. Fix $\epsilon$.

2. Construct the sequence $(l_j)_{j=1,...,n(\epsilon)+1}$ using (23).

3. Approximate the rate function using (22).

For $i$ in 1: $N$, repeat the following steps:

4. Simulate $Y_j$ according to the Gamma distribution $\Gamma_0 (a(l_{j+1}) - a(l_j), b^{(\epsilon)}(l_j))$ for $j = 0, ..., n(\epsilon) - 1$ and $Y_{n(\epsilon)}$ according to the distribution $\Gamma_0 (a(T) - a(n(\epsilon)), b^{(\epsilon)}(n(\epsilon)))$.

5. Compute $X_{t,i}^{(\epsilon)} = \sum_{j=1}^{n(\epsilon)} Y_j 1_{\{t \le l_{j+1}\}}$.

## 6.2 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 $\lambda_j = \mathcal{L}_{X_t}^{-1}(0.01 * j)$, $j = 1, ..., 100$ and we define the relative error on the Laplace transform by:

$$E_L(t) = \frac{1}{100} \sum_{j=1}^{100} \frac{|\mathcal{L}_{X_t}(\lambda_j) - \bar{\mathcal{L}}_{X_t}^{(\epsilon)}(\lambda_j)|}{\mathcal{L}_{X_t}(\lambda_j)}, \quad (39)$$

where $\bar{\mathcal{L}}_{X_t}^{(\epsilon)}(\lambda_j)$ is the empirical Laplace transform of $X_t^{(\epsilon)}$ at point $\lambda_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.

<table><caption>**Table 2** Mean (standard deviation) of the relative errors for the four series representations of Method 2 for $a(t) = t$, $b(t) = (t+1)^2$ and $t = 1$</caption><thead><tr><th></th><th>ILM</th><th>Bond</th><th>Thin</th><th>Rej</th></tr></thead><tbody><tr><td>simulation parameter</td><td><i>B</i> = 10</td><td><i>B</i> = 10</td><td><i>B</i> = 10</td><td><i>B</i> = 10</td></tr><tr><td>cpu time (s)</td><td>3</td><td>1</td><td>1</td><td>1</td></tr><tr><td>mean</td><td>0.0027 (0.0022)</td><td>0.0028 (0.0020)</td><td>0.0246 (0.0033)</td><td>0.0026 (0.0021)</td></tr><tr><td>variance</td><td>0.0087 (0.0066)</td><td>0.0089 (0.0064)</td><td>0.0084 (0.0064)</td><td>0.0086 (0.0062)</td></tr><tr><td>Laplace transform</td><td>0.0017 (0.0011)</td><td>0.0021 (0.0008)</td><td>0.0511 (0.0022)</td><td>0.0021 (0.0010)</td></tr></tbody></table>

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.

<table><caption>**Table 3** Mean (standard deviation) of the relative errors for the different methods of simulation of $X_t$ for $a(t) = t$, $b(t) = (t+1)^2$ and $t = 1$</caption><thead><tr><th></th><th>DSI</th><th>Bond</th><th>Rej</th><th>DRF</th></tr></thead><tbody><tr><td>simulation parameters</td><td>$K^* = 20$</td><td>$B = 10$</td><td>$B = 10$</td><td>$\epsilon = 0.035$</td></tr><tr><td>cpu time (s)</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>mean</td><td>0.0026 (0.0020)</td><td>0.0028 (0.0020)</td><td>0.0026 (0.0021)</td><td>0.0028 (0.0022)</td></tr><tr><td>variance</td><td>0.0108 (0.0077)</td><td>0.0089 (0.0064)</td><td>0.0086 (0.0062)</td><td>0.0088 (0.0070)</td></tr><tr><td>Laplace transform</td><td>0.0020 (0.0013)</td><td>0.0021 (0.0008)</td><td>0.0021 (0.0010)</td><td>0.0017 (0.0011)</td></tr></tbody></table>

<table><caption>**Table 4** Mean (standard deviation) of the relative errors for the different methods of simulation of $X_t$ for $a(t) = 2t$, $b(t) = \frac{1}{(t+1)^2}$ and $t = 2$</caption><thead><tr><th></th><th>DSI</th><th>Bond</th><th>Rej</th><th>DRF</th></tr></thead><tbody><tr><td>simulation parameters</td><td>$K^* = 30$</td><td>$B = 10$</td><td>$B = 10$</td><td>$\epsilon = 0.25$</td></tr><tr><td>cpu time (s)</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>mean</td><td>0.0014 (0.0011)</td><td>0.0014 (0.0011)</td><td>0.0014 (0.0011)</td><td>0.0015 (0.0011)</td></tr><tr><td>variance</td><td>0.0301 (0.0070)</td><td>0.0056 (0.0040)</td><td>0.0056 (0.0041)</td><td>0.0055 (0.0041)</td></tr><tr><td>Laplace transform</td><td>0.0038 (0.0012)</td><td>0.0010 (0.0007)</td><td>0.0011 (0.0007)</td><td>0.0011 (0.0007)</td></tr></tbody></table>

<table><caption>**Table 5** Mean (standard deviation) of the relative errors for the different methods of simulation of $X_t$ for $a(t) = 1 - \exp(-t)$, $b(t) = 1 - \exp(-t - 1)$ and $t = 1$</caption><thead><tr><th></th><th>DSI</th><th>Bond</th><th>Rej</th><th>DRF</th></tr></thead><tbody><tr><td>simulation parameters</td><td>$K^* = 20$</td><td>$B = 30$</td><td>$B = 30$</td><td>$\epsilon = 0.021$</td></tr><tr><td>cpu time (s)</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>mean</td><td>0.0030 (0.0024)</td><td>0.0032 (0.0024)</td><td>0.0031 (0.0025)</td><td>0.0030 (0.0023)</td></tr><tr><td>variance</td><td>0.0085 (0.0066)</td><td>0.0092 (0.0070)</td><td>0.0086 (0.0065)</td><td>0.0088 (0.0066)</td></tr><tr><td>Laplace transform</td><td>0.0019 (0.0012)</td><td>0.0019 (0.0011)</td><td>0.0019 (0.0013)</td><td>0.0019 (0.0011)</td></tr></tbody></table>

## 6.3 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 Sub-section 5.4).

An EGP $X \sim \Gamma(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 \in [0.01: 0.1: 4]$ and $t = 2$, up to
a precision lower than $3 \cdot 10^{-6}$, obtained for $K = 800$ and $\varepsilon = 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 \in [0.01 : 0.1 : 4]$ and both Post-Widder's and DRFI methods. For $H \simeq 100$ in Post-Widder's method and $(\varepsilon = 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 $\varepsilon$ and increasing $K$ whereas the highest precision by Post-Widder's method is obtained for $H \simeq 100$ in the present example, with a clearly lower precision for $H = 200$. This problem of poor convergence when $H \to \infty$ for Post-Widder's method has already been observed in Masol and Teugels (2010).

<table><caption>**Table 6** Mean absolute error on $F_{X_t}(x)$ for $a(t) = t$, $b(t) = (t+1)^2$, $t = 2$ and $x \in [0.01 : 0.1 : 4]$</caption><thead><tr><th>Method</th><th>Mean absolute error</th><th>cpu time</th></tr></thead><tbody><tr><td>Post-Widder (H = 50)</td><td>0.0033</td><td>2</td></tr><tr><td>Post-Widder (H = 100)</td><td>0.0015</td><td>4</td></tr><tr><td>Post-Widder (H = 110)</td><td>0.0016</td><td>4</td></tr><tr><td>Post-Widder (H = 200)</td><td>0.0645</td><td>7</td></tr><tr><td>DRFI (ε = 10<sup>-2</sup>, K = 12)</td><td>0.0014</td><td>3</td></tr><tr><td>DRFI (ε = 10<sup>-2</sup>, K = 200)</td><td>2.10<sup>-5</sup></td><td>15</td></tr><tr><td>DRFI (ε = 10<sup>-2</sup>, K = 800)</td><td>2.10<sup>-5</sup></td><td>15</td></tr><tr><td>DRFI (ε = 10<sup>-3</sup>, K = 800)</td><td>&lt; precision</td><td>30</td></tr></tbody></table>

Keeping the same $a(t)$ and $b(t)$, three approximations of the cdf of $X_{10}$ are next plotted in Figure 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 ($\varepsilon = 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}}^{(\varepsilon)}$ at each continuity point of $F_{X_{10}}$. In the zoomed part of Figure 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 Figure 5 (right) shows that Post-Widder's method seems a little less accurate than the DRFI approximation.

## 7 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^{(\varepsilon)}(t)$) and the other two ($b^{(\varepsilon,+)}(t)$, $b^{(\varepsilon,-)}(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
<figure><img src="image_6.png" /><figcaption>**Fig. 4** The cdf of $X_t$ as a function of $x$ for $a(t) = t$, $b(t) = (t+1)^2$ and $t = 10$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 5** The pdf of $X_t$ as a function of $x$ for $a(t) = t$, $b(t) = (t+1)^2$ and $t = 10$</figcaption></figure>

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 \to \infty$). 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.

References

Abdel-Hameed M (1975) A gamma wear process. IEEE Trans Rel 24(2):152-153

Barlow RE, Proschan F (1965) Mathematical theory of reliability, Classics in Applied Mathematics, vol 17.
Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, with contributions by Larry
C. Hunter, 1996

Çinlar E (1980) On a generalization of gamma processes. J Appl Probab 17:467-480

Çinlar E (2011) Probability and stochastics, Graduate Texts in Mathematics, vol 261. Springer, New York

Çinlar E, Bažant ZP, Osman E (1977) Stochastic process for extrapolating concrete creep. J of the Engrg Mech Div, ASCE 103(EM6):1069-1088

Daley DJ, Vere-Jones D (2007) An introduction to the theory of point processes - Volume II: general theory and structure, vol 2. Springer Science & Business Media

Dykstra RL, Laud P (1981) A bayesian nonparametric approach to reliability. Ann Stat 9(2):356-367

Guida M, Postiglione F, Pulcini G (2012) A time-discrete extended gamma process for time-dependent degradation phenomena. Reliab Eng Syst Safe 105:73-79

Imai J, Kawai R (2013) Numerical inverse Lévy measure for infinite shot noise series representation. J Comput Appl Math 253:264-283

Ishwaran H, James LF (2004) Computational methods for multiplicative intensity models using weighted gamma processes: proportional hazards, marked point processes, and panel count data. J Amer Statist Assoc 99(465):175-190

Laud PW, Smith AFM, Damien P (1996) Monte Carlo methods for approximating a posterior hazard rate process. Stat Comput 6(1):77-83

Masol V, Teugels JL (2010) Numerical accuracy of real inversion formulas for the Laplace transform. J Comput Appl Math 233:2521-2533

Moschopoulos PG (1985) The distribution of the sum of independent gamma random variables. Ann Inst Statist Math 37(1):541-544

Nadarajah S (2008) A review of results on sums of random variables. Acta Appl Math 103(2):131-140

van Noortwijk JM (2009) A survey of the application of gamma processes in maintenance. Reliab Eng Syst Safe 94:2-21

Peppas K (2011) Advanced Trends in Wireless Communications, InTech, chap Performance Analysis of Maximal Ratio Diversity Receivers over Generalized Fading Channels

Rausand M, Høyland A (2004) System Reliability Theory: Models, Statistical Methods, and Applications, 2nd edn. Wiley Series in Probability and Statistics, Wiley-Interscience [John Wiley & Sons], Hoboken, NJ

Rosinski J (2001) Series representations of Lévy processes from the perspective of point processes. In: Lévy processes - Theory and Applications, Eds Barndorff-Nielsen, OE, Mikosch, T, Resnick, SI, Birkhäuser pp 401-415

Veillette MS, Taqqu MS (2011) A technique for computing the pdfs and cdfs of nonnegative infinitely divisible random variables. J Appl Probab 48:217-237