1 Introduction

The Lévy process is a fundamental tool for the study of continuous time stochastic phenomena Bertoin (1997), Barndorff-Nielsen et al. (2001), Applebaum (2009). The most familiar example is of course Brownian motion, but it is well known that the Gaussian assumption is inadequate for modelling of real-world phenomena, which are often more heavy-tailed than the Gaussian. Applications are found in areas such as financial modelling Mandelbrot (1963), Fama (1965), Cont and Tankov (2003), communications Azzaoui and Clavier (2010), Fahs and Abou-Faycal (2012), de Freitas et al. (2017), Liebeherr et al. (2012), Shevlyakov and Kim (2006), Warren and Thomas (1991), signal processing Nikias and Shao (1995), image analysis Achim et al. (2001), Achim et al. (2006) and audio processing Godsill and Rayner (1998), Lombardi and Godsill (2006). Non-Gaussian heavy-tailed effects are also important in the climatological sciences Katz and Brown (1992), Katz et al. (2002), in the medical sciences Chen et al. (2010) and for the understanding of sparse modelling/compressive sensing Unser et al. (2014), Unser et al. (2014), Unser and Tafti (2014), Amini and Unser (2014), Carrillo et al. (2016), Lopes (2016), Zhou and Yu (2017), Tzagkarakis (2009), Achim et al. (2010).

In this paper we study a very general class of non-Gaussian Lévy processes, the generalised inverse Gaussian process Barndorff-Nielsen et al. (2001), Eberlein and Hammerstein (2004), a process which can capture various degrees of heavy-tailed behaviour, including the gamma process, the inverse Gaussian and the reciprocal-gamma process, as well as processes that lie somewhere ‘in between’ these edge cases Eberlein and Hammerstein (2004). Our work also enables direct simulation of the mean–variance mixtures of GIG processes, leading to the generalised hyperbolic (GH) processes Eberlein (2001). Important sub-classes of the GH process include the asymmetric Student-t process which has previously been intractable for simulation and inference to our knowledge (this class is entirely distinct from the Student-t processes developed in the machine learning literature Shah et al. 2014). The use of asymmetric Student-t models in financial econometric applications is discussed in Zhu and Galbraith (2010).

Rosiński (2001) presents a generalised shot-noise representation of non-Gaussian Lévy processes, and it is this approach that we develop here for the GIG and GH processes. Our previous work using the shot noise representation has focussed on stable law processes and their applications in engineering, see Lemke and Godsillc (2015), Riabiz et al. (2020), Godsill et al. (2019) and references therein. There have been relevant developments in various special cases over recent years, including Barndorff-Nielsen (1997) who present the theory of normal-inverse Gaussian (NIG) processes, Rydberg (1997) who present approximate sampling methods for the NIG case, while Barndorff-Nielsen and Shephard (2001) give applications of series based methods for non-Gaussian Ornstein–Uhlenbeck (OU) processes. Zhang (2011b) provided rejection sampling and series based simulations for the GIG OU process using the concepts of Rosiński (2001), but these are applied to a different Background driving Lévy process (BDLP) than our work and indeed require evaluation of a generally intractable integral involving Bessel functions. In addition Zhang (2011a), Qu et al. (2021), Grabchak (2019) have provided exact simulation methods for the related class of tempered stable (TS) processes, while recent relevant literature on GIG Lévy fields can be found in Barth and Stein (2017). It should be noted that the shot noise method involves infinite series of decreasing random variables, and in practice the series must be truncated at some finite limit in simulation. This truncation is the only approximation involved in our methods, which are otherwise exact, and in most cases the series are found to converge rapidly as the number of terms increases.

The distribution of the GIG Lévy process at \(t=1\), the GIG distribution, possesses a three-parameter probability density function defined for positive real random variables as follows Eberlein and Hammerstein (2004):

$$\begin{aligned} f_{GIG}(x) = \left( \frac{\gamma }{\delta } \right) ^{\lambda } \frac{1}{2 K_{\lambda }(\delta \gamma )} x^{\lambda -1} e^{-\frac{1}{2} \left( \delta ^2 x^{-1} + \gamma ^2x \right) } {\mathcal {I}}_{x > 0} \end{aligned}$$
(1)

where \(\lambda \in {\mathbb {R}}\), \(K_{\lambda }(\cdot )\) is the modified Bessel function of the second kind and \({{\mathcal {I}}}\) is the indicator function. The support of parameters \(\gamma \) and \(\delta \) depends on the sign of \(\lambda \) such that:

$$\begin{aligned}&\delta \ge 0, \quad \gamma> 0, \quad \text {if} \, \lambda> 0, \\&\delta> 0, \quad \gamma> 0, \quad \text {if} \, \lambda = 0, \\&\delta > 0, \quad \gamma \ge 0, \quad \text {if} \, \lambda < 0. \end{aligned}$$

as discussed in Eberlein and Hammerstein (2004). Random variate generation algorithms for a GIG variable are studied in Devroye (2014), Hörmann and Leydold (2013).

It is shown in Barndorff-Nielsen and Halgreen (1977) that the GIG distribution is infinitely divisible and hence can be the distribution of a Lévy process at time \(t=1\). Furthermore, particular values of the parameters lead to special cases of the GIG distribution such as the inverse Gaussian \((\lambda = -1/2)\), Gamma (\(\delta =0\), \(\lambda >0\)) and the reciprocal-Gamma (\(\gamma = 0\), \(\lambda <0\)) distributions, and in these limit cases the normalising constant is replaced by the normalising constant of these well-known distributions.

The principal contribution of this paper is to provide a comprehensive suite of methods for simulation of GIG processes, without the need for evaluation of intractable integrals, beyond pointwise evaluation of the relevant Bessel functions. An auxiliary variables approach transforms the univariate GIG point process into a bivariate point process having the GIG process as its marginal by construction and requiring no explicit evaluation of integrals. We derive tractable dominating measures for the augmented GIG Lévy density, hence leading to a random thinning methodology for generation of jumps of the underlying marginal GIG process. The whole procedure is carried out through generation of familiar random variables (from the gamma family) and point processes (both gamma and tempered stable). In addition we are able to bound the average acceptance rates of the random variate generation and also to provide upper and lower bounds on the GIG Lévy density and the corresponding Jaeger integral. Finally the whole methodology is made accessible to practitioners through the publication of Matlab and Python codeFootnote 1 to implement the GIG sampling schemes.

Section 2 presents the necessary preliminaries for simulation of Lévy processes and their corresponding point processes, using a generalised shot-noise approach. Section 3 gives the specific forms for the GIG Lévy density and derives bounds on these densities, as well as a generic thinning method for tractable sampling of the underlying point processes, and presents in detail the simulation method for two different parameter ranges of the process. Section 4 presents example simulations, compared with exact simulations of the GIG random variable, and finally Sect. 5 discusses the application of GIG processes in simulation and inference for more general processes including the generalised hyperbolic process.

2 Shot noise representations, Lévy densities and thinning

In this section we present the required preliminaries about Lévy processes and shot-noise simulation of those processes, using the framework of Rosiński (2001). The characteristic function for a general non-negative-valued Lévy process W(t) having no drift or Brownian motion part is given by Kallenberg (2002), Corollary 15.8, as:

$$\begin{aligned}&E \left[ \exp (iuW(t)) \right] \nonumber \\&\quad = \exp \left( t \left[ \int _{(0,\infty )} (e^{iuw} -1)Q(dw) \right] \right) \end{aligned}$$
(2)

where Q is a Lévy measure on \([0,\infty )\) and satisfying \(\int _{(0,\infty )}(1 \wedge x)Q(dx)<\infty \) (Bertoin 1997, p.72). By the Lévy–Ito integral representation, we may express W(t) directly as:

$$\begin{aligned} W(t)&= \int _{(0,\infty )} w N([0, t], dw) \end{aligned}$$
(3)

Here N is a bivariate point process having mean measure \(Lebesgue \times Q\) on \([0,T]\times (0,\infty )\), which can be conveniently expressed as:

$$\begin{aligned} N=\sum _{i=1}^\infty \delta _{V_i,X_i} \end{aligned}$$

where \(\{V_i\in [0,T]\}\) are i.i.d. uniform random variables which give the times of arrival of jumps, \(\{X_i\in (0,\infty )\}\) are the size of the jumps and \(\delta _{V,X}\) is the Dirac measure centred at time V and jump size X.

Substituting N directly into (3) leads to

$$\begin{aligned} W(t)=\sum _{i=1}^\infty X_i{{{\mathcal {I}}}}_{V_i\le t},\,\,\,as \end{aligned}$$
(4)

The almost sure convergence of this series to \(\{W(t)\}\) is proven in Rosiński (2001). Most of the new material in this paper is concerned with generating the jump sizes \(\{X_i\}\) for the GIG case. We will thus usually refer to the point process as just the set of jump sizes, \(N=\{X_i\}\), but of course corresponding jump times \(\{V_i\}\) will need to be sampled as above in all cases in order to realise the process according to (4).

In principle it may be possible to simulate directly from the point process N, but in practice the potentially infinite number of jumps in any finite time interval make this impossible. If the jumps can be ordered by size, however, it may be possible to simulate all of the significant jumps and ignore or approximate the residual error from omitting the smallest jumps. It turns out that this can be done in a very convenient way, using the well-known approach of Ferguson and Klass (1972), Rosiński (2001), Wolpert and Ickstadt (1998). The starting point is the simulation of the epochs of a unit rate Poisson process \(\{{\Gamma }_i\}_{i\ge 1}\). The intensity function of this process is of course \(\lambda _{\Gamma }=1\), and the resulting realisations contain almost surely an infinite number of terms. We know how to simulate an arbitrarily large number of ordered terms from this process by repeatedly generating standard exponential random variables and calculating the cumulative sum of these to obtain an ordered \({\Gamma }_i\) sequence. Now, define the upper tail probability of the Lévy measure as

$$\begin{aligned} Q^+(x)=Q ([x,\infty )) \end{aligned}$$

and a corresponding non-increasing function h(.) as the inverse tail probability:

$$\begin{aligned} h(\gamma )={Q^+}^{-1}(\gamma )\,. \end{aligned}$$

Then, the following point process converges a.s. to N Rosiński (2001):

$$\begin{aligned} \sum _{i=1}^\infty \delta _{V_i,h({\Gamma }_i)} \end{aligned}$$

and the corresponding convergent representation of the Lévy process is (neglecting the compensator term \(c_i\), which is zero for all cases considered here):

$$\begin{aligned} W(t)=\sum _{i=1}^\infty h({\Gamma }_i){{{\mathcal {I}}}}_{V_i\le t},\,\, as \end{aligned}$$
(5)

Thus, to generate points directly from N by this method it is necessary to be able to compute \(h(\gamma )={Q^+}^{-1}(\gamma )\) explicitly. In the cases considered here it will not be possible to do this and instead an indirect approach is adopted, using thinning or rejection sampling Lewis and Shedler (1979), Rosiński (2001) of more tractable point processes for which \(h(\gamma )\) is directly available. In particular, we seek a ‘bounding’ process \(N_0\) having Lévy measure \(Q_0\) and satisfying \(dQ_0(x)/dQ(x)\ge 1\) \(\forall x\in (0,\infty )\); then, realisations of \(N_0\) are thinned with probability \(dQ(x)/dQ_0(x)\) in order to obtain samples from the desired process N. In all cases considered here the Lévy measure possesses a density function, which we also denote by Q(x) (using the minor abuse of notation ‘\(dQ(x)=Q(x)dx\)’) and the required bounding condition is then \(Q_0(x)\ge Q(x)\), with associated thinning probability \(Q(x)/Q_0(x)\).

Two bounding processes are used extensively in the methods of this paper, the tempered stable (TS) and gamma processes, which may be simulated by standard shot noise methods as follows:

2.1 Tempered stable point process

In the tempered stable (TS) case the Lévy density is, for \(\alpha \in (0,1)\) Shephard and Barndorff-Nielsen (2001) (see also Rosiński 2007)

$$\begin{aligned} Q(x) = Cx^{-1-\alpha } e^{-\beta x}, \quad \quad \quad x > 0 \end{aligned}$$
(6)

Several possible approaches to simulation of sample paths from tempered stable processes were proposed in Rosiński (2001), Rosiński (2007) and compared in Imai and Kawai (2011), which recommends the use of the inverse Lévy measure approach over thinning and rejection sampling methods. For the inverse Lévy measure approach the tail probability may be calculated in terms of gamma functions, but is not easily inverted, and numerical approximations are needed Imai and Kawai (2011). We thus adopt a thinning approach Rosiński (2001) in which the Lévy density is factorised into a positive \(\alpha \)-stable process with Lévy density \(Q_0(x)=Cx^{-1-\alpha }\) Samorodnitsky and Taqqu (1994) and a tempering function \(e^{-\beta x}\). The stable law process has tail mass \(Q_0^+(x)=\frac{C}{\alpha }x^{-\alpha }\) and hence \(h(\gamma )=Q_0^+(\gamma )^{-1}=\left( \frac{\alpha \gamma }{C}\right) ^{-1/\alpha }\). Having simulated the stable point process with rate function \(Q_0\), points are individually selected (thinned) with probability \(e^{-\beta x}\), otherwise deleted.

The recipe for generating the tempered stable process \(N_{TS}\) is then:

figure a

In practice i is truncated at some large value and an approximate realisation results.

2.2 Gamma process

The Lévy density for the Gamma process is:

$$\begin{aligned} Q(x)={C}{x^{-1}}e^{-\beta x}, \quad \quad \quad x > 0 \end{aligned}$$

Rosiński (2001) suggests four possible sampling schemes for this process. The tail probability is the exponential integral, which would require numerical inversion, see Wolpert and Ickstadt (1998). Instead we adopt the thinning (also called the ‘rejection method’ in Rosiński 2001) version of this in which a dominating point process is chosen as \(Q_0(x)=\frac{C}{x}(1+\beta x)^{-1}\). The tail probability for this is \(Q_0^+(x)=C\log \left( \beta ^{-1} x^{-1}+1\right) \) and hence \(h(\gamma ) = \frac{1}{\beta \left( \exp (\gamma / C) - 1 \right) }\). Points are then thinned with probability \((1+\beta x) \exp (-\beta x) \le 1\). As reported in Rosiński (2001), this thinning method is highly effective, with very few point rejections observed. The recipe for generating the gamma process \(N_{Ga}\) is then given in Algorithm 2.

figure b

As before i is truncated at some large value, yielding an approximate realisation of the gamma process.

3 Simulation from the generalised inverse Gaussian Lévy process

In this section the GIG Lévy density is presented and a general scheme is presented that will enable simulation from the GIG process.

The Lévy density for the generalised inverse Gaussian (GIG) process is given by Eberlein and Hammerstein (2004):

$$\begin{aligned} \frac{e^{-x\gamma ^2/2}}{x}\left[ \int _0^\infty \frac{e^{-xy}}{\pi ^2y|H_{|\lambda |}(\delta \sqrt{2y})|^2}dy\!+\!\text {max}(0,\lambda )\right] ,\,\,x\!>\!0\nonumber \\ \end{aligned}$$
(7)

where \(H_{\lambda }(z)=J_{\lambda }(z)+iY_{\lambda }(z)\) is the Hankel function of the first kind (z is always real-valued in the current context), \(J_{\lambda }(z)\) is the Bessel function of the first kind, and \(Y_{\lambda }(z)\) is the Bessel function of the second kind.

The GIG Lévy density comprises two terms: the initial integral, which we denote \(Q_{GIG}(x)\):

$$\begin{aligned} Q_{GIG}(x)=\frac{e^{-x\gamma ^2/2}}{x}\int _0^\infty \frac{e^{-xy}}{\pi ^2y|H_{|\lambda |}(\delta \sqrt{2y})|^2}dy,\,\,x>0 \end{aligned}$$

added to a second term, present only for \(\lambda >0\):

$$\begin{aligned} \frac{e^{-x\gamma ^2/2}}{x}\text {max}(0,\lambda ),\,\,x>0 \end{aligned}$$
(8)

This second term is a Gamma process that may be straightforwardly simulated using the standard methods of 2.2 and added to the simulation of points from the first term \(Q_{GIG}(x)\). Hence we will neglect this second term for now. It will be convenient to rewrite \(Q_{GIG}(x)\) using the substitution \(z=\delta \sqrt{2y}\) as:

$$\begin{aligned} Q_{GIG}(x)&=\frac{e^{-x\gamma ^2/2}}{x}\int _0^\infty \frac{e^{-xy}}{\pi ^2y|H_{|\lambda |}(\delta \sqrt{2y})|^2}dy \\&=\frac{2 e^{-x\gamma ^2/2}}{\pi ^2x}\int _0^\infty \frac{e^{-\frac{z^2x}{2\delta ^2}}}{z|H_{|\lambda |}(z)|^2}dz \end{aligned}$$

Note that this integral is known elsewhere as the Jaeger integral, which finds application in diffusive transport Freitas (2018). Beyond our direct interest in GIG processes, there is significant interest in approximating these integrals accurately, and our bounding approach is likely to provide accurate bounds and approximations which may be compared with those proposed in Freitas (2018).

Throughout this paper we will consider only positive and real-valued variables x, y and z, as the complex versions of these functions are not required in the present context.

Our general scheme is to consider the following intensity function associated with a bivariate point process on \((0,\infty )\times (0,\infty )\):

$$\begin{aligned} Q_{GIG}(x,z)=\frac{2 e^{-x\gamma ^2/2}}{\pi ^2x}\frac{e^{-\frac{z^2x}{2\delta ^2}}}{z|H_{|\lambda |}(z)|^2} \end{aligned}$$
(9)

which has by construction the GIG process as its marginal:

$$\begin{aligned} Q_{GIG}(x)=\int _{0}^\infty Q_{GIG}(x,z)dz \end{aligned}$$

We propose to simulate points directly from this bivariate process, hence avoiding any direct evaluation of the Jaeger integral. Since \(Q_{GIG}(x,z)\) is intractable itself for simulation, dominating processes \(Q^0_{GIG}(x,z)\ge Q_{GIG}(x,z)\) are constructed and points sampled from \(Q^0_{GIG}(x,z)\) are thinned with probability \(Q_{GIG}(x,z) / Q^0_{GIG}(x,z)\). Thus a significant part of our approach is in constructing suitable dominating functions that are tractable for simulation, and this is achieved by studying the properties of the Jaeger integral.

The first set of bounds are obtained from the basic properties of the Hankel function Watson (1944) and will lead in particular to a simulation algorithm for the case \(|\lambda |\ge 0.5\). Properties of the modulus of the Hankel function are mainly derived from the Nicholson integral representation Watson (1944),

$$\begin{aligned} |H_\nu (z)|^2=\frac{8}{\pi ^2}\int _0^\infty K_0(2z\sinh t)\cosh (2\nu t) dt \end{aligned}$$

where \(K_0\) is the modified Bessel function of the second kind. In particular this leads to an asymptotic (\(z\rightarrow \infty \)) series expansion (Watson 1944, Sect. 13.75):

$$\begin{aligned}&|H_{\nu }(z)|^{2} \sim \frac{2}{\pi z} \bigg ( 1+\frac{4\nu ^{2}-1}{2(2z)^{2}} + \frac{1\cdot 3(4\nu ^{2}-1)(4\nu ^{2}-9)}{2\cdot 4(2z)^{4}} \\&\quad + \frac{1\cdot 3\cdot 5(4\nu ^{2}-1)(4\nu ^{2}-9)(4\nu ^{2}-25)}{2\cdot 4\cdot 6(2z)^{6}}\,\,... \bigg ) \,, \end{aligned}$$

from which the well-known asymptotic value is obtained,

$$\begin{aligned} \underset{z\rightarrow \infty }{\text {lim}} z|H_{\nu }(z)|^{2}=\frac{2}{\pi } \end{aligned}$$

From the Nicholson integral the following important property can also be derived (Watson 1944, Sect. 13.74):

Property 1

For any real, positive z, \(z|H_{\nu }(z)|^{2}\) is a decreasing function of z for \(\nu >0.5\) and an increasing function of z when \(0<\nu <0.5\).

For the case \(\nu =0.5\), \(z|H_{\nu }(z)|^{2}=2/\pi \), i.e. a constant.

These basic properties are used to prove the following bounds:

Theorem 1

For any positive z and fixed \(|\lambda |\ge 0.5\), the following bound applies:

$$\begin{aligned} Q_{GIG}(x,z)\le \frac{ e^{-x\gamma ^2/2}}{\pi x}{e^{-\frac{z^2x}{2\delta ^2}}} \end{aligned}$$
(10)

and, for \(|\lambda |\le 0.5\),

$$\begin{aligned} Q_{GIG}(x,z)\ge \frac{ e^{-x\gamma ^2/2}}{\pi x}{e^{-\frac{z^2x}{2\delta ^2}}} \end{aligned}$$
(11)

with equality holding in both cases when \(|\lambda |=0.5\).

Proof

Property 1 and the limiting value of \(2/\pi \) lead to

$$\begin{aligned} z|H_\nu (z)|^2{\left\{ \begin{array}{ll} \ge \frac{2}{\pi },&{} \nu \ge 0.5\\ \le \frac{2}{\pi },&{}\nu \le 0.5\end{array}\right. } \end{aligned}$$
(12)

The results follow by substitution of (12) into \(Q_{GIG}(x,z)\) (9). \(\square \)

Corollary 1

The bound in (10), applicable for \(|\lambda |\ge 0.5\), can be rewritten as

$$\begin{aligned} \frac{ e^{-x\gamma ^2/2}}{\pi x}{e^{-\frac{z^2x}{2\delta ^2}}}=\frac{ {\delta {\Gamma }(1/2)}e^{-x\gamma ^2/2}}{\sqrt{2}\pi x^{3/2}}\sqrt{{Ga}}\left( z|1/2,\frac{x}{2\delta ^2}\right) \end{aligned}$$
(13)

where \(\sqrt{Ga}\) is the square-root gamma density, i.e. the density of \(X^{0.5}\) when \(X\sim Ga(x|\alpha ,\beta )\), having probability density function

$$\begin{aligned} \sqrt{\text{ Ga }}(x|\alpha , \beta ) = \frac{2 \beta ^{\alpha }}{{\Gamma }(\alpha )} x^{2\alpha -1} e^{-\beta x^2} \end{aligned}$$

Remark 1

It can be seen immediately that (13) corresponds marginally to a tempered stable process in x, and conditionally to a \(\sqrt{Ga}\) density for z, a tractable feature that will enable sampling from the dominating bivariate point process. In fact, this decomposition and that of Corollary 2 are the key point for our new GIG simulation methods. We are here decomposing the bivariate point process in (x, z) as a marked point process comprising a marginal Poisson process for x and a conditional mark random variable z (Cont and Tankov 2003, Sect. 2.6.4). The important point here is that the conditional distributions of z|x are integrable and tractable probability density functions and hence a marginal-conditional sampling procedure may be adopted constructively to simulate from the joint point process in (x, z).

Remark 2

Integration with respect to z leads to a simple upper bound on the GIG Lévy density:

$$\begin{aligned} \frac{ {\delta {\Gamma }(1/2)}e^{-x\gamma ^2/2}}{\sqrt{2}\pi x^{3/2}} \end{aligned}$$
(14)

which was also obtained by Zhang (2011b), and used in a rejection sampling procedure that requires a direct evaluation of the Jaeger integral, in contrast with our approach. For \(|\lambda |\le 0.5\), a similar argument leads to a lower bound with the same formula as (14).

Remark 3

The first bound for \(\nu \ge 0.5\) in (12) will be used shortly as an envelope function for rejection sampling in the case \(|\lambda |\ge 0.5\) using Corollary 1. The second bound in (12) for \(\nu \le 0.5\) will be adapted in Theorem 2 to develop more sophisticated bounds in this parameter range.

A second set of bounds can be stated as follows:

Theorem 2

Choose a point \(z_0\in (0,\infty )\) and compute \(H_0=z_0|H_{\nu }(z_0)|^2\). This will define the corner point on a piecewise lower or upper bound. Define now \(z_1 = \left( \frac{ 2^{1-2\nu }\pi }{{\varGamma }^2(\nu )}\right) ^{1/(1-2\nu )}\) and define the following functions:

$$\begin{aligned} A(z)={\left\{ \begin{array}{ll}\frac{2}{\pi }\left( \frac{z_1}{z}\right) ^{2\nu -1},&{}z < z_1 \\ \frac{2}{\pi },&{}z\ge z_1 \end{array}\right. } \end{aligned}$$

and

$$\begin{aligned} B(z)={\left\{ \begin{array}{ll}H_0\left( \frac{z_0}{z}\right) ^{2\nu -1},&{}z<z_0\\ H_0,&{}z\ge z_0\end{array}\right. } \end{aligned}$$

Then, for \(0<\nu \le 0.5\),

$$\begin{aligned} A(z)\ge z|H_{\nu }(z)|^2\ge B(z) \end{aligned}$$
(15)

and for \(\nu \ge 0.5\),

$$\begin{aligned} A(z)\le z|H_{\nu }(z)|^2\le B(z) \end{aligned}$$
(16)

with all inequalities becoming equalities when \(\nu =0.5\), and both A(z) bounds (left side inequalities) becoming tight at \(z=0\) and \(z=\infty \).

Proof

First note that \(z|H_{\nu }(z)|^2\) is an increasing function for \(0< \nu < 0.5\) and decreasing for \(\nu > 0.5\) as stated in Property 1. When combined with the tight bound of \(2/\pi \) as \(z\rightarrow \infty \) (12) this leads directly to the \(2/\pi \) components of A(z). Similarly the first part of the bound \(A(z)= z({\varGamma }(\nu )/\pi )^2\left( \frac{2}{z}\right) ^{2\nu }\), is justified since Freitas (2018) proves that \(z^{2\nu }|H_{\nu }(z)|^2\) is increasing when \(\nu >0.5\) and decreasing when \(\nu <0.5\), with a tight bound of \(({\varGamma }(\nu )/\pi )^2 2^{2\nu }\) as \(z\rightarrow 0\). Thus the left hand side of the inequalities, A(z), are established.

The right hand side of the inequalities, B(z), are proven by the monotonicity and sign of gradient for both \(z^{2\nu }|H_{\nu }(z)|^2\) and \(z|H_{\nu }(z)|^2\). We may choose an arbitrary corner point \((z_0,H_0)\) on \(z|H_{\nu }(z)|^2\) and monotonicity of \(z|H_{\nu }(z)|^2\) (Property 1) immediately implies that, for \(z\ge z_0\),

$$\begin{aligned} z|H_{\nu }(z)|^2{\left\{ \begin{array}{ll}\ge H_0,&{}0<\nu \le 0.5\\ \le H_0,&{}\nu \ge 0.5\end{array}\right. } \end{aligned}$$

Moreover, monotonicity and sign of gradient of \(z^{2\nu }|H_{\nu }(z)|^2\) Freitas (2018) imply that

$$\begin{aligned} z^{2\nu }|H_{\nu }(z)|^2{\left\{ \begin{array}{ll}\ge z^{2\nu -1}H_0,&{}0<\nu \le 0.5\\ \le z^{2\nu -1} H_0,&{}\nu \ge 0.5\end{array}\right. } \end{aligned}$$

from which the bounds B(z) are immediately obtained. \(\square \)

Fig. 1
figure 1

Plot of Bessel function bounds, \(\nu =0.3\). \(z_0\) set equal to \(z_1\)

Full size image
Fig. 2
figure 2

Plot of Bessel function bounds, \(\nu =0.8\). \(z_0\) set equal to \(z_1\)

Full size image

Remark 4

The proven bounds can be clearly visualised on a log-log plot, which highlights the asymptotic behaviour of the functions and bounds as \(z\rightarrow 0 \) and \(z\rightarrow \infty \), see Figs. 1 and 2.

Remark 5

The choice of \(z_0\in (0,\infty )\) is arbitrary. However, its value will impact the tightness of the bounding functions B(z) and will hence impact the effectiveness of our subsequent sampling algorithms and integral approximations. A suitable generic choice was found to be at the same z value as the corner point of the corresponding A(z) bounds, i.e. set \(z_0=z_1\), as plotted in Figs. 1 and 2, though further optimisation may be possible in applications.

Corollary 2

For the case \(|\lambda |<0.5\), the right hand bound in (15), B(z), can be substituted into \(Q_{GIG}(x,z)\) (9) and rearranged to obtain:

$$\begin{aligned} Q_{GIG}(x,z)&= \frac{2 e^{-x\gamma ^2/2}}{\pi ^2x}\frac{e^{-\frac{z^2x}{2\delta ^2}}}{z|H_{|\lambda |}(z)|^2}\nonumber \\&\le \frac{ 2 e^{-x\gamma ^2/2}}{\pi ^2x} {\left\{ \begin{array}{ll} \frac{z^{2|\lambda |-1}e^{-\frac{z^2x}{2\delta ^2}}}{H_0z_0^{2|\lambda |-1}},&{}0<z<z_0\\ \frac{e^{-\frac{z^2x}{\delta ^2}}}{H_0},&{}\text {otherwise} \end{array}\right. } \end{aligned}$$
(17)

Furthermore, the two piecewise sections of this bound may be factorised in terms of left- and right-truncated square-root gamma densities:

$$\begin{aligned}&Q_{GIG}(x,z) \le \frac{ e^{-x\gamma ^2/2}}{\pi ^2x^{1+|\lambda |}}\frac{(2\delta ^2)^{|\lambda |}\gamma (|\lambda |,z_0^2x/(2\delta ^2))}{H_0z_0^{2|\lambda |-1}} \nonumber \\ {}&\frac{{\Gamma }(|\lambda |)\sqrt{\text{ Ga }} (z||\lambda |,x/(2\delta ^2))}{\gamma (|\lambda |,z_0^2x/(2\delta ^2))}{{{\mathcal {I}}}}_{0<z<z_0}\nonumber \\ {}&\quad +\frac{ e^{-x\gamma ^2/2}}{\pi ^2x^{3/2}}\frac{(2\delta ^2)^{0.5}{\Gamma }(0.5,z_0^2x/(2\delta ^2))}{H_0} \nonumber \\ {}&\frac{{\Gamma }(0.5)\sqrt{\text{ Ga }} (z|0.5,x/(2\delta ^2))}{{\Gamma }(0.5,z_0^2x/(2\delta ^2))} {{{\mathcal {I}}}}_{z\ge z_0} \end{aligned}$$
(18)

where the truncated square-root gamma densities are, with their associated normalising constants:

$$\begin{aligned}&\frac{{\Gamma }(|\lambda |)\sqrt{\text{ Ga }} (z||\lambda |,x/(2\delta ^2))}{\gamma (|\lambda |,z_0^2x/(2\delta ^2))}{{{\mathcal {I}}}}_{0<z<z_0}\,, \\&\frac{{\Gamma }(0.5)\sqrt{\text{ Ga }} (z|0.5,x/(2\delta ^2))}{{\Gamma }(0.5,z_0^2x/(2\delta ^2))} {{{\mathcal {I}}}}_{z\ge z_0} \end{aligned}$$

and lower/upper incomplete gamma functions are defined in the usual way, for \(\text {Re}(s) > 0\), as:

$$\begin{aligned} \gamma (s,x) = \int _{0}^{x} t^{s-1} e^{-t} dt,\,\,\,\, {\Gamma }(s,x) = \int _{x}^{\infty } t^{s-1} e^{-t} dt \end{aligned}$$
Fig. 3
figure 3

Plot of optimised upper and lower bounds on \(Q_{GIG}(x)\), for various parameter settings. We use the following shorthand for the optimised \(Q^B\) bounds: \({Q^B_{GIG}}^{*}(x) = \underset{z_0\in (0,\infty )}{\mathrm {max/min}}\left\{ Q^B_{GIG}(x)\right\} \), where \(\text {min}\) applies for \(|\lambda |<0.5\) and \(\text {max}\) applies for \(|\lambda |\ge 0.5\). Overlaid also are the simple approximation (12) and the trends \(x^{-3/2}\) and \(x^{-(1+|\lambda |)}\)

Full size image

\(Q_{GIG}(x,z)\) can thus be split up into two tractable point processes for simulation: a first, \(N_1\), comprising a modified tempered \(|\lambda |\)-stable process with truncated \(\sqrt{Ga}\) conditional for z; and a second, \(N_2\), comprising a modified tempered 0.5-stable process with a further truncated \(\sqrt{Ga}\) conditional for z. We will later use this union of point processes as the dominating Lévy measure in simulation of the \(|\lambda |<0.5\) case.

Corollary 3

Integration of the bounding function (18) with respect to z allows a more sophisticated estimate for the Jaeger integral and therefore the Lévy density for the GIG process, which is an upper bound for \(|\lambda |<0.5\) and a lower bound for \(|\lambda |>0.5\), following the direction of the right hand inequalities in (18) and (16). A similar procedure inserting the left hand (A(z)) inequalities from (18) and (16) yields the corresponding lower (\(|\lambda |<0.5\)) and upper (\(|\lambda |>0.5\)) bounds. Define first the bounding functions \(Q^A_{GIG}(x)\) and \(Q^B_{GIG}(x)\), corresponding to the bounds A(z) and B(z) as follows:

$$\begin{aligned} Q^A_{GIG}(x)&= \frac{ e^{-x\gamma ^2/2}}{\pi x^{1+|\lambda |}}\frac{(2\delta ^2)^{|\lambda |}\gamma (|\lambda |,z_1^2x/(2\delta ^2))}{2z_1^{2|\lambda |-1}} \\&\quad +\frac{ e^{-x\gamma ^2/2}}{\pi x^{3/2}}\frac{(2\delta ^2)^{0.5}{\Gamma }(0.5,z_1^2x/(2\delta ^2))}{2} \\ Q^B_{GIG}(x)&= \frac{ e^{-x\gamma ^2/2}}{\pi ^2 x^{1+|\lambda |}}\frac{(2\delta ^2)^{|\lambda |}\gamma (|\lambda |,z_0^2x/(2\delta ^2))}{H_0z_0^{2|\lambda |-1}}\\ {}&\quad +\frac{ e^{-x\gamma ^2/2}}{\pi ^2x^{3/2}}\frac{(2\delta ^2)^{0.5}{\Gamma }(0.5,z_0^2x/(2\delta ^2))}{H_0}\end{aligned}$$

These are obtained by substituting the bounds A(z) or B(z) into the expression for \(Q_{GIG}(x,z)\) (9) and integrating with respect to z.

Noting that \(z_0\in (0,\infty )\) can be chosen arbitrarily, the \(Q^B_{GIG}\) estimates may be improved by optimising with respect to \(z_0\) to either maximise (\(|\lambda |< 0.5\)) or minimise (\(|\lambda |> 0.5\)) \(Q^B_{GIG}(x)\) at each point x.

Then, following the direction of the inequalities in (15) and (16) we obtain the following estimates of \(Q_{GIG}(x)\):

$$\begin{aligned}&Q^A_{GIG}(x)\le Q_{GIG}(x)\!\le \! \underset{z_0\in (0,\infty )}{\mathrm {min}}\left\{ Q^B_{GIG}(x)\right\} ,\,\,\,\,|\lambda |\le 0.5 \\&Q^A_{GIG}(x)\!\ge \! Q_{GIG}(x)\ge \underset{z_0\in (0,\infty )}{\mathrm {max}}\left\{ Q^B_{GIG}(x)\right\} ,\,\,\,\,|\lambda |\ge 0.5 \end{aligned}$$

with all inequalities becoming equalities at \(|\lambda |=0.5\). Optimisation to achieve the required maximum and minimum functions can be achieved by numerical search.

Example plots are given in Fig. 3 in which the optimised lower and upper bounds are plotted, showing very close agreement (often indistinguishable by eye) in their estimation of \(Q_{GIG}(x)\), over various parameter ranges and large x ranges (i.e. the bounds are visually quite tight). Overlaid are the \(x^{-3/2}\) and \(x^{-(1+|\lambda |)}\) trends and also the simple approximation of (14). This simple approximation diverges significantly from our proposed bounds, particularly when \(|\lambda |\) is not close to 0.5 (when \(|\lambda |=0.5\) all bounds are equal to the true function \(Q_{GIG}(x)\) for the inverse Gaussian (IG) process).

3.1 Simulation of GIG processes with \(|\lambda |\ge 0.5\)

In this section the specific algorithm applied for simulation in the case \(|\lambda |\ge 0.5\) is detailed. In this parameter range it has been found very effective to use Theorem 1 and Corollary 1, so that the dominating bivariate process is:

$$\begin{aligned} Q^0_{GIG}(x,z)=\frac{ e^{-x\gamma ^2/2}}{\pi x}{e^{-{(z^2x)}/{(2\delta ^2)}}} \end{aligned}$$

which is simulated as a marked point process having factorised (marginal-conditional) intensity function (Corollary 1) as follows:

$$\begin{aligned} Q^0_{GIG}(x,z)=\frac{ {\delta {\Gamma }(1/2)}e^{-x\gamma ^2/2}}{\sqrt{2}\pi x^{3/2}}\sqrt{{Ga}}\left( z|1/2,\frac{x}{2\delta ^2}\right) \end{aligned}$$

The thinning probability for points drawn from the dominating point process is then:

$$\begin{aligned} \frac{Q_{GIG}(x,z)}{Q^0_{GIG}(x,z)}=\frac{2}{\pi z|H_{|\lambda |}(z)|^2} \end{aligned}$$

The procedure for generation of points from the process with intensity function \(Q_{GIG}(x)\) is given in Algorithm 3. In the case of \(\lambda > 0\), points generated from the process with intensity function shown in Eq. (8) are added to the set of points coming from \(Q_{GIG}(x)\) to obtain jump sizes from the GIG process.

figure c

The procedure is completed as before by independently drawing an associated jump time \(v\sim {{{\mathcal {U}}}}[0,T]\) for each accepted jump size (‘point’) \(x\in N\). The value of the GIG Lévy process at t is given in Eq. (4).

Fig. 4
figure 4

Plot of upper and optimised lower bounds on \(\rho (x)\), for various \(|\lambda |>0.5\), showing that large jumps are rejected while with increasing probability small jumps are accepted. \(\delta =0.1\) in all cases and bounds do not depend on \(\gamma \)

Full size image

3.1.1 Acceptance rate

The mean acceptance probability for fixed x is:

$$\begin{aligned} \rho (x)&=E \left[ \frac{2}{\pi z|H_{|\lambda |}(z)|^2} \right] \\&=\int _{0}^{\infty } \sqrt{\text {Ga}}(z|1/2,x/(2\delta ^2))\frac{2}{\pi z|H_{|\lambda |}(z)|^2}dz \end{aligned}$$

which may be bounded using Theorem 2. For the current parameter setting (\(|\lambda |>0.5\)) we have \(A(z) \le z |H_{|\lambda |}(z)|^2 \le B(z)\). Lower and upper bounds are obtained by substituting A(z) and B(z) and then by direct integration to give:

$$\begin{aligned} \begin{aligned}&\frac{2}{\pi H_0} \left( \left( \frac{z_0^2x}{2\delta ^2}\right) ^{0.5-|\lambda |} \frac{\gamma (|\lambda |, \frac{z_0^2x}{2\delta ^2})}{{\varGamma }(0.5)} + \frac{{\varGamma }(0.5, \frac{z_0^2x}{2\delta ^2})}{{\varGamma }(0.5)} \right) \\&\quad \le \rho (x) \\&\quad \le \left( \frac{z_1^2x}{2\delta ^2}\right) ^{0.5-|\lambda |} \frac{\gamma (|\lambda |, \frac{z_1^2x}{2\delta ^2})}{{\varGamma }(0.5)} + \frac{{\varGamma }(0.5, \frac{z_1^2x}{2\delta ^2})}{{\varGamma }(0.5)} \end{aligned} \end{aligned}$$
(19)

As noted before, the corner point \(z_0\in (0,\infty )\) may be chosen arbitrarily, while \(z_1\) is fixed. Hence the lower bound may be optimised with respect to \(z_0\) at each x value to obtain a tighter lower bound:

$$\begin{aligned}&\underset{z_0\in (0,\infty )}{\text {max}}\bigg \{ \frac{2}{\pi H_0} \bigg ( \bigg (\frac{z_0^2x}{2\delta ^2}\bigg )^{0.5-|\lambda |} \frac{\gamma (|\lambda |, z_0^2x/(2\delta ^2))}{{\varGamma }(0.5)} \\&\quad + \frac{{\varGamma }(0.5, z_0^2x/(2\delta ^2))}{{\varGamma }(0.5)} \bigg ) \bigg \} \le \rho (x) \end{aligned}$$

The (numerically optimised) lower bound and (fixed) upper bound on the mean acceptance rate \(\rho (x)\) are visualised for different parameter settings in Fig. 4. It can be observed, as expected from the inequalities (19), that the bounds are tight as \(x \rightarrow 0\) and \(x \rightarrow \infty \), and that acceptance rates get lower as \(|\lambda |\) increases. In all cases though they indicate that only large jumps (large x values) will be rejected and that at some point all jumps below a certain x value are highly likely to be accepted on average.

In principle we can go a little further than the acceptance rate for fixed x and compute overall acceptance rates for the algorithm, and quantities such as the expected number of acceptances/rejections overall. In particular it may be informative to calculate the expected number of accepted points if the underlying Poisson process \(\{{\varGamma }_i\}\) is truncated to values of less than say c, i.e. the process is approximated using the random truncation \(\sum _{{\varGamma }_i\le c}h({\varGamma }_i)\). Since \(\{{\varGamma }_i\}\) is a unit rate Poisson process, the expected number of accepted points at a truncation level c is:

$$\begin{aligned} N_R(c)=\int _0^c \alpha \left( h({\Gamma }) \right) \rho \left( h({\Gamma }) \right) d{\Gamma } \end{aligned}$$
(20)

Here \(h({\varGamma })\) is the non-increasing function used to generate a particular point process (5), and in the case of Algorithm 3 this is the tempered stable process with \(h({\varGamma })=\left( \frac{\alpha {\varGamma }}{C}\right) ^{-1/\alpha }\). \(\alpha (x)\) is the acceptance probability in any accept–reject step prior to generating z; in this case this is the term \(\exp (-\beta x)\) from the accept–reject step from the tempered stable process, where \(\beta =\gamma ^2/2\) (Algorithm 1). Computing the integrand in (20) gives the point process intensity of acceptances as a function of the unit rate Poisson process’s time evolution. Figure 5 illustrates this by plotting the upper and lower bounds on this intensity for two values of \(\lambda <-0.5\). Note that on average many fewer points are accepted at the start of the series for the larger magnitude \(\lambda \), which lends some theoretical weight to the empirically observed property that the series are more slowly converging as \(\lambda \) becomes increasingly negative (and the process more light-tailed), see experimental simulations section for further results on this effect. We also plot in Fig. 6 the total number of expected rejections, computed as \(c-N_R(c)\) by numerically integrating (20) with the upper and lower bounds substituted for the integrand. This serves to illustrate again the significant impact of \(\lambda \) on total number of rejected points for a particular truncation limit c. Also overlaid on this figure are the actual mean number of rejections averaged over 1000 simulations of the relevant process with different truncation levels c, illustrating that the true expectation lies between the theoretical bounds.

Fig. 5
figure 5

Point acceptance intensity as a function of time evolution in the unit rate Poisson process \(\{{\varGamma }_i\}\). \(\gamma =0.2\), \(\delta =0.1\)

Full size image
Fig. 6
figure 6

Expected number of point process rejections as a function of the truncation level c of \(\{{\varGamma }_i\}\). \(\gamma =0.2\), \(\delta =0.1\), showing upper and lower bounds, plus average number of points rejected in 1000 sample paths generated by the shot noise method

Full size image

3.2 Simulation of GIG processes with \(0<|\lambda |<0.5\)

The simple bound of Theorem 1 cannot be straightforwardly applied to simulation the GIG process for \(|\lambda | <0.5\). Instead we use the more sophisticated piecewise bounds of Theorem 2 and Corollary 2, which give us the dominating point process (18) that can be split into two independent point processes \(N_1\) and \(N_2\) as follows:

$$\begin{aligned}&N_1:\,\,\,\,\, \frac{ e^{-x\gamma ^2/2}}{\pi ^2x^{1+|\lambda |}}\frac{(2\delta ^2)^{|\lambda |}\gamma (|\lambda |,z_0^2x/(2\delta ^2))}{H_0z_0^{2|\lambda |-1}} \\&\frac{{\Gamma }(|\lambda |)\sqrt{\text{ Ga }} (z||\lambda |,x/(2\delta ^2))}{\gamma (|\lambda |,z_0^2x/(2\delta ^2))}{{{\mathcal {I}}}}_{z<z_0}\\&N_2:\,\,\,\,\, \frac{ e^{-x\gamma ^2/2}}{\pi ^2x^{3/2}}\frac{(2\delta ^2)^{0.5}{\Gamma }(0.5,z_0^2x/(2\delta ^2))}{H_0} \\&\frac{{\Gamma }(0.5)\sqrt{\text{ Ga }} (z|0.5,x/(2\delta ^2))}{{\Gamma }(0.5,z_0^2x/(2\delta ^2))}{{{\mathcal {I}}}}_{z\ge z_0} \end{aligned}$$

Each may be simulated using a thinned tempered stable process for x and a truncated \(\sqrt{Ga}\) density for z. Having simulated each pair (x, z) they are accepted with probability equal to the ratio \(Q_{GIG}(x,z)/Q^0_{GIG}(x,z)\). Owing to the piecewise form of the acceptance probability the two processes may be treated independently, including accept reject steps and the point process union taken at the final step to achieve the final GIG samples. The sampling procedure for point process \(N_1\) is given in Algorithm 4.

figure d

Similarly, for \(N_2\), the procedure is given in Algorithm 5.

figure e

Finally, the set of points \(N=N_1\cup N_2\) is a realisation of jump sizes from the point process having intensity \(Q_{GIG}(x)\). The procedure is completed as before by generating independent jump times \(v\sim {{{\mathcal {U}}}}[0,T]\) for all points \(x\in N\).

3.2.1 Acceptance rates

Theorem 2 is used once again to find lower bounds on the expected acceptance rates. In this case the bound to apply is \(z|H_{\lambda }(z)|^2 \le A(z)\).

For \(N_1\) the average acceptance rate \(\rho _1(x)\) is

$$\begin{aligned} \rho _1(x)&= E \left[ \frac{H_0}{|H_{|\lambda |}(z)|^2 \left( \frac{z^{2 |\lambda | }}{z_0^{2|\lambda |-1}}\right) } \right] \\&= \int _0^{z_0}\frac{{\varGamma }(|\lambda |)\sqrt{\text {Ga}} (z||\lambda |,x/(2\delta ^2))}{\gamma (|\lambda |,z_0^2x/(2\delta ^2))} \\&\quad \frac{H_0}{z|H_{|\lambda |}(z)|^2} \left( \frac{z_0}{z}\right) ^{2|\lambda |-1} dz \end{aligned}$$

which may be lower bounded by substituting the bound \(z|H_{\lambda }(z)|^2 \le A(z)\) from Theorem 2 and by direct integration to give:

$$\begin{aligned} \rho _1(x) \ge \frac{H_0 z_0^{2|\lambda |-1} \pi ^2}{2^{2|\lambda |} {\varGamma }^2(|\lambda |)} \end{aligned}$$

Similarly, for \(N_2\) the average acceptance rate \(\rho _2(x)\) is

$$\begin{aligned} \rho _2(x)&= E \left[ \frac{H_0}{z|H_{|\lambda |}(z)|^2} \right] \\&= \int _{z_0}^{\infty } \frac{{\varGamma }(0.5)\sqrt{\text {Ga}} (z|0.5,x/(2\delta ^2))}{{\varGamma }(0.5,z_0^2x/(2\delta ^2))} \frac{H_0}{z|H_{|\lambda |}(z)|^2} dz \end{aligned}$$

and again lower bounds are obtained by substituting \(z|H_{\lambda }(z)|^2 \le A(z)\) from Theorem 2 and direct integration to give:

$$\begin{aligned} \rho _2(x) \ge \frac{\pi H_0}{2} \end{aligned}$$

We note that both of these bounds are in fact independent of the value of x. If we use B(z) to obtain an upper bound on the acceptance rates \(\rho _1\) and \(\rho _2\), both bounds are found to be 1, i.e. no useful information is gleaned from the upper bound.

3.3 Sampling from the marginal point process envelope

In steps 2 above we require simulation of the marginal point process for x in the dominating bivariate point process. For \(N_1\) this has intensity:

$$\begin{aligned} \begin{aligned} Q_1(x)&= \frac{e^{-x\gamma ^2/2}}{\pi ^2x^{1+|\lambda |}}\frac{(2\delta ^2)^{|\lambda |}\gamma (|\lambda |,z_0^2x/(2\delta ^2))}{H_0z_0^{2|\lambda |-1}} \\&= \frac{e^{-x\gamma ^2/2}}{\pi ^2x}\frac{\gamma (|\lambda |,z_0^2x/(2\delta ^2))}{(z_0^2 x/(2\delta ^2))^{|\lambda |} H_0 z_0^{-1}} \end{aligned} \end{aligned}$$
(21)

and we propose two possible methods for simulating this point process.

The first and most basic method uses the fact that the lower incomplete gamma function is upper bounded by the complete gamma function, i.e. \(\gamma (|\lambda |,z_0^2x/(2\delta ^2))\le {\varGamma }(|\lambda |)\), which allows generation from \(Q_1\) by thinning of a tempered stable process. We thus have the following upper bounding envelope as a tempered stable density:

$$\begin{aligned} Q_1(x)\le {\bar{Q}}_1^1(x) = \frac{e^{-x \gamma ^2/2} {\varGamma }(|\lambda |)}{\pi ^2 x^{1+|\lambda |}} \frac{(2\delta ^2)^{|\lambda |}}{H_0 z_0^{2|\lambda |-1}} \end{aligned}$$

This process may be routinely sampled by thinning of a positive tempered \(|\lambda |\)-stable process as described in 2.1.

Having generated points from the tempered stable envelope function the process is thinned with probability:

$$\begin{aligned} \frac{Q_1(x)}{{\bar{Q}}_1^1(x)}=\frac{\gamma (|\lambda |,z_0^2x/(2\delta ^2))}{{\varGamma }(|\lambda |)} \end{aligned}$$

Experimentally we found that the acceptance rates were quite low for this bound, meaning that fairly long tempered stable series had to be generated (but see note below about the case \(\gamma =0\)), so a more sophisticated bound was sought, as follows.

The second and more effective method employs the following bound on the lower incomplete gamma function (Neuman 2013, Theorem 4.1):

$$\begin{aligned} \frac{a\gamma (a,x)}{x^a} \le \frac{(1+ae^{-x})}{(a+1)} \end{aligned}$$

and so

$$\begin{aligned} \frac{\gamma (|\lambda |,z_0^2 x/(2\delta ^2))}{(z_0^2 x/(2\delta ^2))^{|\lambda |}} \le \frac{(1+|\lambda |e^{-z_0^2 x/(2\delta ^2)})}{|\lambda |(|\lambda |+1)} \end{aligned}$$
(22)

so that the point process intensity function is upper bounded by

$$\begin{aligned} Q_1(x)&\le \frac{e^{-x\gamma ^2/2}}{\pi ^2x}\frac{z_0(1+|\lambda |e^{-z_0^2 x/(2\delta ^2)})}{|\lambda |(1+|\lambda |)H_0}={{\bar{Q}}_1^{2a}(x)} \end{aligned}$$
(23)
$$\begin{aligned}&\le \frac{e^{-x\gamma ^2/2}}{\pi ^2x}\frac{z_0}{|\lambda |H_0}={{\bar{Q}}_1^2(x)} \end{aligned}$$
(24)

where (24) follows from \(0<e^{-z_0^2 x/(2\delta ^2)} \le 1\) for \(x \in [0, \infty )\). The bounding process \({\bar{Q}}^2_1(x)\) may be simulated as a single gamma process, having parameters \(a = \frac{z_0}{\pi ^2 |\lambda |H_0}\) and \(\beta =\gamma ^2/2\). It is then thinned with the following probability:

$$\begin{aligned} \frac{Q_1(x)}{{\bar{Q}}_1^2(x)}&= \frac{|\lambda | \gamma (|\lambda |, z_0^2 x/(2\delta ^2))}{\left( \frac{z_0^2 x}{2 \delta ^2} \right) ^{|\lambda |}} \end{aligned}$$
(25)

The whole procedure is given in Algorithm 6.

figure f

Note that a slightly more refined bounding procedure may be implemented by observing that (23) is the intensity function of the union of two gamma processes, which may be independently simulated and their union thinned with probability \(\frac{Q_1(x)}{{\bar{Q}}_1^{2a}(x)}\), at the expense of generating one additional gamma process.

Note also that the second method using Algorithm 6 does not work for \(\gamma =0\) since the bounding point process cannot be simulated in this case; instead the first method with bound \({\bar{Q}}_1^1(x)\) must be used for the \(\gamma =0\) case. This second method, using bounding function \({\bar{Q}}_1^2(x)\), is regarded as the superior approach for all parameter settings except for \(\gamma =0\) since it has increased acceptance probabilities and also relies on an underpinning Gamma process rather than a \(|\lambda |\)-stable process as its basis, and the series representation of the Gamma process is known to converge very rapidly to zero compared to the \(|\lambda |\)-stable process Rosiński (2001), Barndorff-Nielsen and Shephard (2001).

Similarly, for \(N_2\) we follow a thinning approach. The marginal dominating Lévy density for \(N_2\) is:

$$\begin{aligned} Q_2(x) = \frac{e^{-x\gamma ^2/2}}{\pi ^2x^{3/2}}\frac{(2\delta ^2)^{0.5}{\varGamma }(0.5,z_0^2x/(2\delta ^2))}{H_0} \end{aligned}$$
(26)

Using the complete gamma function as part of an upper bound for this density, we can sample \(N_2\) as a tempered stable intensity \(\frac{e^{-x\gamma ^2/2} (2\delta ^2)^{0.5} {\varGamma }(0.5)}{\pi ^2H_0 x^{3/2}}\), and apply thinning with probability \({\varGamma }(0.5, z_0^2 x/(2\delta ^2))/{\varGamma }(0.5)\). This procedure is found to work well, and the procedure is summarised in Algorithm 7.

figure g

Rejection rates could potentially be further improved by more tightly upper bounding the term \({\varGamma }(0.5, \frac{z_0^2 x}{2\delta ^2})\). The following simple bound is suitable, and valid for \(0<s<1\):

$$\begin{aligned} {\Gamma }(s,x) \le x^{s-1} e^{-x} \end{aligned}$$

This was not explored in the simulations as Algorithm 7 was found to perform well.

4 Simulations

In this section example simulations from the new method are generated up to \(t=T\), where \(T=1\), and compared with a random variable generator for the GIG distribution Devroye (2014), Statovic (2017) (a ‘ground truth’ simulation for comparing the distribution at \(t=T\), but unlike our method, not able to generate the entire path of the process). In our new shot noise case, \(M=1000\) terms are generated for each point process series representation and \(10^6\) random realisations of the GIG random variable are generated to produce Figs. 9-15, in which examples of the principal parameter ranges and edge case \(\gamma =0\) are presented. Note that M represents the total number of terms generated for each of the underlying tempered stable and gamma processes, rather than the final number of accepted terms. \(10^4\) random samples were generated in each case for our new shot noise method and \(10^6\) samples for the ‘ground truth’ method of Devroye (2014), Statovic (2017) at \(t=1\). In Figs. 16-19 example pathwise simulations are plotted for several parameter settings, drawing 30 independent realisations of the GIG process in each case, setting \(T=1\). For cases where \(\lambda >0\) the additional term in the Lévy density, \(\lambda \frac{e^{-x\gamma ^2/2}}{x}\), see (7), is generated as an independent Gamma process and the union of the these points with those from \(Q_{GIG}\) is taken to generate the final simulated process. Once the union of all the points has been taken, to form a set of jump sizes \(\{x_i\}\), independent and uniformly distributed jump times \(\{v_i\in [0,T]\}\) are generated and the process path may be realised as

$$\begin{aligned} \sum _i x_i{{{\mathcal {I}}}}_{v_i\le t},\,\,\,t\in [0,T] \end{aligned}$$

As previously discussed, the resulting process approximates a GIG Lévy process. By comparing the distribution of the Lévy process at \(t=1\) and the distribution of GIG random variables generated using Devroye (2014), Statovic (2017), the quality of this approximation may be evaluated, at least at the end-point of the interval \(t=T=1\). A possible statistical test in this case is the two-sample Kolmogorov–Smirnov (KS) test since in most cases the CDF of a GIG random variable is not available in closed form. The KS test checks whether the two samples come from the same distribution by reporting the KS statistic, computed as the supremum of the absolute difference between the empirical CDFs of the two samples Hodges (1958).

The results of such a KS test are reported in Fig. 7, for values of \(|\lambda |<-0.5\), varying shot noise series length M, and using random sample sizes of 10000. Remaining parameters were \(\gamma =0.1\), \(\delta =2\). It can be seen that the convergence of the KS statistic is slower as \(\lambda \) becomes more negative, as alluded to in earlier analysis. The final p-values at \(M=10000\) are typically well above 0.1, indicating that the hypothesis cannot be rejected for this M, at a confidence level of 0.1, with considerably earlier convergence for the smaller absolute values of \(\lambda \).

For \(0<|\lambda |<0.5\) we observed very rapid convergence of the KS statistics to acceptable levels as M increases, typically within a few dozen terms in the series, so we do not display these results.

For \(\lambda >0.5\), Fig. 8 gives the equivalent set of KS curves, showing more rapid convergence than their negative counterparts (presumably because of the additional gamma process that is mixed into the point process) and again showing p-values typically well above 0.1 for \(M=10000\).

Fig. 7
figure 7

KS statistics for \(\lambda < -0.5\) and various series lengths M

Full size image
Fig. 8
figure 8

KS statistics for \(\lambda > 0.5\) and various series lengths M

Full size image

Of course, we do recognise that our generated GIG processes are approximate, and so it should always be possible to reject the null hypothesis for sufficiently large sample sets, but nevertheless these tests do give a helpful indication of the new algorithm’s performance for different ranges of \(\lambda \). Similar behaviours were observed for other settings of parameters \(\gamma \) and \(\delta \).

Finally the computational burden of the algorithms is linear in M, the number of terms in the series. Moreover there are no random waiting times since the accept–reject steps are performed as one-off tests for each point generated in the series. The most significant contributors to the load are the calculation of the various acceptance probabilities and the sampling of the auxiliary variables z. The Matlab and Python code runs at similar speed on standard laptop platforms. On a Microsoft Surface (Intel(R) Core(TM) i7-1065G7 CPU @ 1.30GHz 1.50 GHz) the time for computation was very roughly \(10^{-6}M\)s per random process realisation with \(|\lambda |\ge 0.5\) and approximately double that for the more complex \(|\lambda |<0.5\) case. Note also that in Algorithm 3 some significant parallelisation is possible, since Steps 1.-3. do not depend on \(\lambda \). Thus it would be possible to simulate many different realisations for different \(\lambda \) by performing one set of steps 1.-3. and then performing step 4. independently for each \(\lambda \) value required. This could be of significant value in ensemble and Monte Carlo inference procedures, as well as for rapid visualisation of different \(\lambda \) scenarios.

Fig. 9
figure 9

Simulation comparison between the shot noise generated GIG process and GIG random variates, \(\lambda =-0.1\), \(\gamma =0.1\), \(\delta =2\). Left hand panel: QQ plot comparing our shot noise method (x-axis) with random samples of the GIG density generated using the ‘ground truth’ method. Right hand panel: Normalised histogram density estimate for our method compared with the true GIG density function

Full size image
Fig. 10
figure 10

Simulation comparison between the shot noise generated GIG process and GIG random variates, \(\lambda =-0.4\), \(\gamma =0.5\), \(\delta =1\), \(10^6\) random samples. Left hand panel: QQ plot comparing our shot noise method (x-axis) with random samples of the GIG density generated using the ‘ground truth’ method of Devroye (2014), Statovic (2017). Right hand panel: normalised histogram density estimate for our method compared with the true GIG density function

Full size image
Fig. 11
figure 11

Simulation comparison between the shot noise generated GIG process and GIG random variates, \(\lambda =-1\), \(\gamma =0.5\), \(\delta =4\). Left hand panel: QQ plot comparing our shot noise method (x-axis) with random samples of the GIG density generated using the ‘ground truth’ method of Devroye (2014), Statovic (2017). Right hand panel: normalised histogram density estimate for our method compared with the true GIG density function

Full size image
Fig. 12
figure 12

Simulation comparison between the shot noise generated GIG process and GIG random variates, \(\lambda =-0.3\), \(\gamma =0\), \(\delta =4\). Left hand panel: QQ plot comparing our shot noise method (x-axis) with random samples of the GIG density generated using the ‘ground truth’ method of Devroye (2014), Statovic (2017). Right hand panel: Normalised histogram density estimate for our method compared with the true GIG density function

Full size image
Fig. 13
figure 13

Simulation comparison between the shot noise generated GIG process and GIG random variates, \(\lambda =-1\), \(\gamma =0\), \(\delta =4\). Left hand panel: QQ plot comparing our shot noise method (x-axis) with random samples of the GIG density generated using the ‘ground truth’ method of Devroye (2014), Statovic (2017). Right hand panel: Normalised histogram density estimate for our method compared with the true GIG density function

Full size image
Fig. 14
figure 14

Simulation comparison between the shot noise generated GIG process and GIG random variates, \(\lambda =1\), \(\gamma =0.4\), \(\delta =4\). Left hand panel: QQ plot comparing our shot noise method (x-axis) with random samples of the GIG density generated using the ‘ground truth’ method of Devroye (2014), Statovic (2017). Right hand panel: normalised histogram density estimate for our method compared with the true GIG density function

Full size image
Fig. 15
figure 15

Simulation comparison between the shot noise generated GIG process and GIG random variates, \(\lambda =0.3\), \(\gamma =0.5\), \(\delta =2\). Left hand panel: QQ plot comparing our shot noise method (x-axis) with random samples of the GIG density generated using the ‘ground truth’ method of Devroye (2014), Statovic (2017). Right hand panel: normalised histogram density estimate for our method compared with the true GIG density function

Full size image
Fig. 16
figure 16

Pathwise simulations of the GIG process, \(\lambda =-0.1\), \(\gamma =0.5\), \(\delta =3\)

Full size image
Fig. 17
figure 17

Pathwise simulations of the GIG process, \(\lambda =-1\), \(\gamma =0.5\), \(\delta =3\) (log-scale)

Full size image
Fig. 18
figure 18

Pathwise simulations of the GIG process, \(\lambda =1\), \(\gamma =0.4\), \(\delta =4\)

Full size image
Fig. 19
figure 19

Pathwise simulations of the GIG process, \(\lambda =-0.4\), \(\gamma =0\), \(\delta =1\) (log-scale)

Full size image

5 Discussion

This paper has presented a generic simulation methodology for GIG processes. The methods are simple and efficient to implement and have good acceptance rates, which will make them of use in applications for practitioners. Moreover, we provide code in Matlab and Python in order to give researchers immediate access to the methods.

Simulation of GIG processes opens up the possibility of simulation and inference for many more complex processes; in particular, a direct extension takes the sampled GIG process and generates other processes within the generalised shot-noise methodology Rosiński (2001):

$$\begin{aligned} W(t)=\sum _{i=1}^{\infty } H(U_{i},{\Gamma }_{i}){{{\mathcal {I}}}}(V_{i}\le t) \end{aligned}$$

where \(U_i\) are iid uniform random variates and \(H(u,\gamma )\) a non-increasing function of \(\gamma \). Of particular interest will be the mean–variance mixture of Gaussians, with:

$$\begin{aligned} H(U,\gamma )\sim {{{\mathcal {N}}}}(\mu _{W} h(\gamma ),\sigma _{W}^{2} h(\gamma )) \end{aligned}$$

and \(J_{i}=h({\varGamma }_{i})\) is the ith ordered jump of the simulated GIG process. This approach, which is an exactly equivalent process to the time-changed Brownian motion description of Barndorff-Nielsen (1997), leads directly to a simulation method for the generalised hyperbolic (GH) process, and its conditionally Gaussian form will enable inference for these processes, using a conditionally Gaussian approach similar in spirit to Lemke and Godsillc (2015), Godsill et al. (2019). Our continued work on these processes will study these inferential approaches with the GIG/GH model and also their use as driving processes for stochastic differential equations, again extending the approach of Godsill et al. (2019) to these more general classes of process.