Saddlepoint Approximation for the Generalized Inverse Gaussian L´evy Process

The generalized inverse Gaussian (GIG) L´evy process is a limit of compound Poisson processes, including the stationary gamma process and the stationary inverse Gaussian process as special cases. However, ﬁtting the GIG L´evy process to data is computationally intractable due to the fact that the marginal distribution of the GIG L´evy process is not convolution-closed. The current work reveals that the marginal distribution of the GIG L´evy process admits a simple yet extremely accurate saddlepoint approximation. Particularly, we prove that if the order parameter of the GIG distribution is greater than or equal to -1, the marginal distribution can be approximated accurately – no need to normalize the saddlepoint density. Accordingly, maximum likelihood estimation is simple and quick, random number generation from the marginal distribution is straightforward by using Monte Carlo methods, and goodness-of-ﬁt testing is undemanding to perform. Therefore, major numerical impediments to the application of the GIG L´evy process are removed. We demonstrate the accuracy of the saddlepoint approximation via various experimental setups.


Introduction
In the family of pure-jump increasing Lévy processes, both the gamma process and the inverse Gaussian process have wide applications. This is mainly because their marginal distributions, namely the gamma and the inverse Gaussian distributions, are convolution-closed and infinitely divisible. Therefore, these two Lévy processes can be easily extended to model non-stationary time-series data (see, e.g., Zhou et al., 2017;Cholette et al., 2019). This work introduces and studies a very general Lévy process, called the generalized inverse Gaussian (GIG) Lévy process, which includes the gamma and inverse Gaussian processes as special cases.
The GIG distribution was proposed byÉtienne Halphen in 1941 and popularized by Ole Barndorff-Nielsen in the 1970s. Barndorff-Nielsen et al. (1978) proved that any GIG distribution with a non-positive power parameter is the distribution of the first hitting time to level 0 for a timehomogeneous diffusion process with state space [0, ∞). This fact suggests the potential use of the GIG distribution as a lifetime distribution or the distribution for times between successive events in a renewal process (Embrechts, 1983). Halgreen (1979) further showed that the GIG distribution is self-decomposable. Therefore, all GIG probability density functions are unimodal (Yamazato, curate in a wide range of settings. In particular, it is often the case that relative errors of these approximations stay bounded in the extreme tails, a desirable property that is not shared by most other types of approximation. Saddlepoint approximations are constructed by performing various operations on the cumulant generating function of a random variable. For the development and discussion of saddlepoint methodology, see Daniels (1954) for details of the density approximation, Lugannani and Rice (1980) and Daniels (1987) for the discussion of a tail area approximation which has a uniform relative error, and Reid (1988) and Goutis and Casella (1999) for a review of saddlepoint techniques.
The main objective of this paper is to show that the marginal distribution of the GIG Lévy process can be well approximated by an analytical function and hence that the GIG Lévy process can be readily applied to model time series data. The remainder of this paper is organized as follows.
Section 2 introduces the GIG distribution and the GIG Lévy process. Section 3 gives a detailed explanation of the saddlepoint approximation and its uniqueness. Section 4 reveals that, although the saddlepoint approximation is fairly accurate, it is not exact, even after normalization. The saddlepoint density is then modified to provide an improved approximation. Section 5 addresses the problems of parameter estimation, random number generation and goodness-of-fit test. Section 6 explores the accuracy of the saddlepoint approximation via simulation. Finally, Section 7 concludes with a summary and remarks.

GIG Distribution & GIG Lévy Process
The density function of the GIG distribution is given by where a > 0, b > 0, and the order parameter λ ∈ R; K λ (·) is a normalizing constant (called modified Bessel function of the second kind): is an exponentially decaying function of v, diverges for all orders at v = 0, and has the property that K −λ (v) = K λ (v). Modified Bessel functions of the second kind of order {0, 1, 2, 3, 4, 5} are shown in Figure 1. We let GIG(λ , a, b) represent the GIG distribution (1). GIG distributions enjoy several nice probabilistic features. For example, if X follows the GIG distribution GIG(λ , a, b), then its reciprocal 1/X follows GIG(−λ , b, a). The GIG distribution includes as special cases the gamma distribution (b = 0 and λ > 0), the inverse gamma distribution (a = 0 and λ < 0), the inverse Gaussian distribution (λ = −0.5) and the hyperbolic distribution (λ = 0). Owning to its infinite divisibility, we can construct a Lévy process from the GIG distribution, herein called GIG Lévy process. We say that the process {X t ; t ≥ 0} is a GIG Lévy process, if the law of X 1 is the GIG distribution GIG(λ , a, b). A Lévy process can be fully determined by the characteristic function of X t which is given by the Lévy-Khintchine formula. In the manner of Dufresne et al. (1991), it is easy to prove that where Π(·) is called the Lévy measure. According to Barndorff-Nielsen and Shephard (2001), the Lévy measure of the GIG Lévy process is absolutely continuous with density is the Bessel function of the first kind, and Y |λ | (·) is the Bessel function of the second kind (see Chapter 9 of Abramovitz and Segun, 1970). For the GIG Lévy process, the arrival rate ∞ 0 Π(dx) is infinite (Morales, 2004). Hence, the GIG Lévy process is a limit of compound Poisson processes, composing of an infinite number of infinitesimal jumps.
The GIG Lévy process includes the stationary gamma process and inverse Gaussian process as special cases. If {X t ; t ≥ 0} is a stationary gamma process (resp., inverse Gaussian process), ∀ t > 0, X t always follows the gamma distribution (resp., the inverse Gaussian distribution). However, for the GIG Lévy process, only X 1 follows the GIG distribution. ∀ 0 ≤ s < t, if t −s = 1, then X t −X s does not follow the GIG distribution. This is because the GIG distribution is not closed under convolution. In other words, if two random variables Z 1 and Z 2 both follow the GIG distribution GIG(λ , a, b), then their sum Z 1 + Z 2 does not follow a GIG distribution. This undesirable feature has restrained the application of the GIG Lévy process in areas where the gamma process and inverse Gaussian process have been employed.
To make the GIG Lévy process a practical model, we need to formulate the density function of X t for any t > 0, which can be derived via the Fourier inversion formula. ∀ t > 0, let f X t (x;θ θ θ ) denote the density function of X t , where θ θ θ = (λ , a, b). The characteristic function of a GIG random variable X is where u ∈ R. Hence, we have where ϕ(u) is the logarithm of E[e iuX ]: Apparently, recovering the density function f X t (x;θ θ θ ) from its characteristic function is not possible explicitly. Hence, in the following section we introduce the saddlepoint method for constructing a closed-form approximation to f X t (x;θ θ θ ).

Brief Explanation
For readability, we introduce here the formal calculations to derive the saddlepoint approximation. Suppose X is a continuous random variable with density f (x). Let ψ(u) denote the moment generating function: ψ(u) = ∞ −∞ e ux f (x)dx. Via the Fourier transform, we have Substituting iu with u and applying the Closed Curve Theorem, we have where τ is within the interval of convergence for ψ(u) which we assume to contain the origin as In what follows, we let k (u, x) and k (u, x) respectively denote the first and second derivative w.r.t. u. Approximate k(u, x) by its Taylor expansion: whereû satisfies k (û, x) = 0 and k (û, x) > 0. Then we have in which, again availing the Closed Curve Theorem, we have set τ =û. Substituting u with iu and employing the idea of Laplace approximation, we have That is, which is the saddlepoint approximation for f (x). Note that log(ψ(u)) is called the cumulantgenerating function of the random variable X.

Saddlepoint Density for f X t (x;θ θ θ )
Let H λ (u) denote the cumulant generating function of X 1 (i.e., t = 1): for u < a 2 , Then the cumulant generating function of X t is tH λ (u). Following (3), the saddlepoint density approximation to f X t (x;θ θ θ ) is: where H λ (u) represents the second derivative of H λ (u) w.r.t. u, andû is the saddlepoint satisfying tH λ (û) − x = 0. Note thatû =û(x,t) is a function of x and t. The first derivative of H λ (u) is Hence,û is obtained by solving w.
The second derivative of H λ (u) is where The pseudo code in Algorithm 1 summarizes the steps for evaluating the saddlepoint density functionf X t (x;θ θ θ ) at any point x.

Existence and Uniqueness of the Saddlepoint
The feasibility of the saddlepoint approximation depends on the existence and uniqueness of the solutionû to tH λ (u) = x and onû satisfying H λ (û) > 0. In this section we discuss the existence and properties of the real root of the equation Proposition 1. There is no real root of the equation H λ (u) = ξ whenever ξ ≤ 0.
Proof. Define a function M(u, ξ ) as which exists only for u < a 2 . Taking partial derivative of M(u, ξ ) w.r.t. u, we have that and that Here and in Proposition 2, we implicitly utilize the dominated convergence theorem to exchange derivatives and integrals. It is clear that the integrand e u(x−ξ ) f (x; λ , a, b) and its partial derivative w.r.t. u are integrable functions of x.
From Equation (10) we know that, when ξ ≤ 0, M (u, ξ ) > 0 for any u < a 2 . Hence, from Equation (9) we know that, when ξ ≤ 0, H λ (u) − ξ > 0 for any u < a 2 ; that is, H λ (u) − ξ = 0 has no real root when ξ ≤ 0. Proof. Note that M (u, ξ ) is strictly increasing with u, because Then for any rootû of the equation When ξ > 0, we rewrite Equation (10) into: For the first integration, we have For the second integration, let m denote the maximum value of f (x; λ , a, b) over the interval then the equation H λ (u) − ξ = 0 has one and only one root.
Equation (9) indicates that lim Proof. We note that, for large values of v, the asymptotic approximation of Then it follows from (6) that: According to Proposition 2, for any ξ > 0, if there exists a root of the equation H λ (u) = ξ , then the root is unique and satisfies H λ (û) > 0. To prove that H λ (u) is a strictly increasing function, we only need to prove that H λ (u) > 0 for any u < a 2 . Let ξ 0 be the point at which lim On one hand, for any 0 < ξ < ξ 0 , we have and therefore there is a unique root to the equation H λ (u) = ξ for any 0 < ξ < ξ 0 . On the other all u < a 2 , and therefore H λ (u) − ξ 0 < 0 for all u < a 2 ; that is, we have 0 < H λ (u) < ξ 0 for any u. Therefore, we can claim that H λ (u) is a bijection from (−∞, a 2 ) to (0, ξ 0 ). Combining with the fact that H λ (û) > 0 everywhere, we can conclude that H λ (u) is a strictly increasing function of H λ (u) = +∞, then with ξ increasing from 0 to +∞, the unique rootû increases • when λ ≥ −1, there is a unique simple root for any ξ > 0; • when λ < −1, there is a unique simple root for any 0 When λ < 0 and λ = −1, define v = (a − 2u)b and we have For small values of v, the asymptotic For small values of v, the asymptotic approximation of To conclude, we have Theorem 1 also explains why the inverse Gaussian distribution, i.e. λ = −0.5, can be (exactly) approximated by the saddlepoint density (Daniels, 1980).

Approximation Error
In many applications, the saddlepoint density does not integrate to one, and hence needs to be normalized. We here point out that generally the saddlepoint approximation for f X t (x;θ θ θ ) is not exact, even after normalization.
We prove by calculating the ratiof X t (x;θ θ θ ) f X t (x;θ θ θ ) . If the ratio is not 1, then we can conclude that the saddlepoint approximationf X t (x;θ θ θ ) is not exact. Moreover, if the ratio changes with x, then we can conclude that even the normalized saddlepoint approximation is not exact. Note that, on one hand, only f X 1 (x;θ θ θ ) has an explicit expression. On the other hand, if the normalizedf X 1 (x;θ θ θ ) is not exact, then for any t > 0, the normalizedf X t (x;θ θ θ ) is not exact either. Therefore, we only need to examine the ratiof X 1 (x;θ θ θ ) f X 1 (x;θ θ θ ) for x > 0. In Figure 2 we plot the ratio for different values of (a, b), with λ fixed at value 2. In each row, b  takes a value from {0.1, 1, 10}, while in each column, a takes a value from {0.1, 1, 10}. The x-axis represents the value of x, and the y-axis represents the value of the ratiof X 1 (x;θ θ θ ) f X 1 (x;θ θ θ ) . Figure  f X t (x;θ θ θ ), even after normalization, is not exact. Moreover, in Figure 2, all the ratios are above 1, while in Figure 3, all the ratios are below 1. All the ratio curves are bounded. Note that the relative difference betweenf X 1 (x;θ θ θ ) and f X 1 (x;θ θ θ ), as measured by the ratio, increases with x; however, with x increasing, the true value f X 1 (x;θ θ θ ) quickly converges to 0. Hence, if we plot the two density functions, they are visually the same (see Section 6). In fact, via exhaustive numerical study, we find that the saddlepoint approximationf X t (x;θ θ θ ) has the following property: Proposition 4. The saddlepoint approximationf X t (x;θ θ θ ) is exact only when λ = −0.5. When λ = −0.5, it is not exact even after normalization. When λ > −0.5, the ratiof X t (x;θ θ θ ) f X t (x;θ θ θ ) is larger than one for any x > 0, and increases with x. When λ < −0.5, the inverse ratio f X t (x;θ θ θ ) f X t (x;θ θ θ ) is larger than one for any x > 0, and increases with x.
Remark 1. We can relate the GIG density (1) to the inverse Gaussian density f (x; −0.5, a, b) by and the expectation is taken w.r.t. f (x; −0.5, a, b). The saddlepoint density can exactly approximate f (x; −0.5, a, b), and E[X λ +0.5 ] is independent of x. Therefore, the approximation error is introduced by the exponentiation x λ +0.5 . When λ increases from -1 to infinity, the change of the differencef X t (x;θ θ θ ) − f X t (x;θ θ θ ) from negative to positive may be caused by the exponentiation x λ +0.5 , which changes from a decreasing function to an increasing function.

A Modified Approximation
Proposition 4 indicates that if we want to improve the approximationf X t (x;θ θ θ ), we have to multiply it by a non-constant factor. Following Section 3.1, let ψ 0 (w) denote the moment generating function of an appropriate distribution that admits an analytic density function f 0 (x). Define k 0 (w, x) = log(ψ 0 (w)) − wx, and letŵ be the unique root of k 0 (w, x) = 0. Instead of approximating k(u, x) by a truncated Taylor expansion, Ait-Sahalia and Yu (2006) approximated it by Both k(u, x) and k 0 (w, x) are strictly convex. Hence, u =û if and only if w =ŵ. Now it is clear that the appropriateness of the benchmark density f 0 (x) means that ψ 0 (w) is defined on a non-trivial interval, andŵ exists wheneverû exists. Moreover, the two local functions, k(u, x) − k(û, x) around u and k 0 (w, x) − k 0 (ŵ, x) aroundŵ, are expected to behave alike.
Now we can treat u as a function of w. By differentiating twice the above equation and setting w to beŵ, we have Then we have or, equivalently, We notice that For the GIG Lévy process, one candidate of f 0 (x) for f X t (x;θ θ θ ) is given by which simply replaces b in Equation (1) with bt 2 . Hence, the quantitiesŵ, ψ 0 (ŵ) and k 0 (ŵ, x) can be readily calculated following Section 3.2. Letf X t (x;θ θ θ ) denote the modified saddlepoint approximation for f X t (x;θ θ θ ). According to Section 3.2, we havē When t = 1, f 0 (x) is identical to f X 1 (x;θ θ θ ), and hencef X 1 (x;θ θ θ ) is exact. Note that, when λ = −0.5, f 0 (x) is the marginal density function of the inverse Gaussian process.
We examine the similarity between f 0 (x) and f X t (x;θ θ θ ), and for ease of exposition we let t = m be an integer. Let i.i.d. random variables {Y 1 , . . . ,Y m } follow the GIG distribution (1), and i.i.d.
random variables {X 1 , . . . , X m } follow the inverse Gaussian distribution f (x; −0.5, a, b). Then the summand ∑ m i=1 Y i follows the distribution f X m (x;θ θ θ ) with the moment generating function given by .
The corresponding f 0 (x) for f X m (x;θ θ θ ) is Let Y be a random variable from f 0 (x), and X a random variable from f X m (x; −0.5, a, b) -a special case of f X m (x;θ θ θ ) and also the density function of ∑ m i=1 X i . The moment generating function of Y is .
We observe that the moment generating function of f X m (x;θ θ θ ) and that of f 0 (x) are alike. Hence, it is expected that the modified saddlepoint densityf X t (x;θ θ θ ) will yield a better approximation.

Parameter Estimation
We now fit the GIG Lévy process to a time series dataset and estimate the unknown parameters.
Letθ θ θ denote the maximum likelihood (ML) estimate of θ θ θ , and Θ a restricted parameter set: Θ = Here we presume that the true value of θ θ θ falls in Θ, as otherwise the saddlepoint method will fail. Although Θ is a restricted set, it is still larger than the parameter sets of the gamma process and the inverse Gaussian process. Starting from time 0, suppose we have data {x 0 , x 1 , x 2 , · · · , x m } collected at time points 0 = t 0 < t 1 < t 2 < · · · < t m . Define ∆x i = x i − x i−1 and ∆t i = t i − t i−1 for i = 1, . . . , m. Then the likelihood for observing ∆x i shall be f X ∆t i (∆x i ;θ θ θ ) which, according to Equation (2), is difficult to evaluate.
We notice from Figure 2 that the approximation error is quite small and uniformly bounded. In Section 6.1, we will corroborate that the saddlepoint density approximates the true density nearly exactly. In fact, due to the round-off error, numerical integration off X t (x;θ θ θ ) over the interval (0, +∞) even gives the value 1. Hence we might perform ML estimation by directly maximizing Here, we use the notationû(x,t) to highlight that the root is a function of x and t. Maximizing the above log-likelihood function is undemanding: the rootû(∆x i , ∆t i ) can be quickly found using, e.g., the bisection method, because H λ (u) is a strictly increasing function of u.

Random Number Generation for f X t (x;θ θ θ )
We generate data fromf X t (x;θ θ θ ) and treat them as sampled from f X t (x;θ θ θ ). Though we can evaluatê f X t (x;θ θ θ ) for any x, the rootû(x,t) does not have an analytic expression. Hence, there is no method available to directly draw i.i.d. samples fromf X t (x;θ θ θ ). We here propose to adopt the Markov chain Monte Carlo (MCMC) technique to generate a sequence of dependent samples, denoted by {X i , i = 1, 2, · · · }, which is a Markov chain with the equilibrium distributionf X t (x;θ θ θ ). Then the Nth element (with N being sufficiently large), X N , can be used as a random sample fromf X t (x;θ θ θ ).
Readers are referred to the excellent texts of Chen et al. (2012) and Meyn and Tweedie (2012) and review papers by Tierney (1994) and Andrieu et al. (2003) for more information on MCMC. We herein develop a sampler based on the Metropolis-Hastings (MH) algorithm.
Let X i = x denote the current state of the Markov chain. The MH sampler is composed of three steps: (1) Generate a proposal sampleẍ from a proposal distribution g(ẍ|x).
To fully develop an MH sampler, we need to specify the proposal distribution g(ẍ|x). Here we work with the Gaussian distribution centered at the current value x with standard deviation σ (> 0).
The value of σ should be subjectively determined to maintain the acceptance rate of proposals in a reasonable range. Note that the target distribution,f X t (x;θ θ θ ), does not have full support, while the Gaussian proposal distribution does. Hence, we need to work with a slightly different proposal distribution -the truncated Gaussian distribution: are respectively the density function and cumulative distribution function of the standard normal distribution. Then the acceptance probability is simply α = min{1,f X t (ẍ;θ θ θ ) Remark 2. If the value of λ inf X t (x;θ θ θ ) is not very large, we can apply the importance sampling with the proposal distribution being the inverse Gaussian distribution (λ = −0.5) or the hyperbolic distribution (λ = 0). Importance sampling is faster than the MH sampler, as the samples are independent. Godsill and Kndap (2021) recently developed a data-generation technique for the GIG Lévy process. They first constructed a bivariate point process having the GIG Lévy process as its marginal, and then developed an acceptance-rejection sampling method for the bivariate point process. By contrast, our simulation method is much simpler.

Goodness-of-Fit Test
To test the goodness of fit, we propose to employ empirical-distribution-function test statistics, e.g., the Kolmogorov-Smirnov test. The idea is to invoke the probability integral transformation and calculate y i = F X ∆t i (∆x i ;θ θ θ ) for i = 1, . . . , m. To calculate {y 1 , · · · , y m }, we need to be able to approximate F X t (x;θ θ θ ) for any t > 0. Again, this can be accomplished by employing the saddlepoint method (Lugannani and Rice, 1980;Daniels, 1987): We might assume {y 1 , · · · , y m } are in ascending order. Denote byF m (y) the empirical distribution function of the data {y 1 , · · · , y m }. The Kolmogorov-Smirnov statistic is defined bŷ the Cramér-von Mises statistic is defined bŷ and the Anderson-Darling (AD) statistic is defined bŷ We then employ the parametric bootstrap technique (Stute et al., 1993) to calculate p-values: 3. Compute y * i =F X ∆t i (x * i ;θ θ θ * ) for i = 1, · · · , m, and then compute the values of the test statistics.
4. Repeat steps 1 to 3 for a large number of times to obtain the corresponding p-values.
6 Numerical Study

Performance of the Saddlepoint Approximation
We examine the performance of the saddlepoint approximation by approximating f X 2 (x;θ θ θ ), i.e., t = 2. To simulate a random value from f X 2 (x;θ θ θ ), we simulate two random values from f X 1 (x;θ θ θ ) and then add them together to obtain a realization of X 2 . Repeat in this manner to simulate a dataset with size 100,000 from f X 2 (x;θ θ θ ). Then the kernel density plot of the simulated data provides an accurate graphical representation of f X 2 (x;θ θ θ ). To examine the accuracy off X 2 (x;θ θ θ ), we just need to plotf X 2 (x;θ θ θ ) within the kernel density plot, which is illustrated in Figures 4, 5 and 6. In Figures 4-7 and 12, the black curve represents the kernel density estimate, the red curve : Density plots of f X 2 (x;θ θ θ ) (black),f X 2 (x;θ θ θ ) (red) andf X 2 (x;θ θ θ ) (green), when b increases.
We now investigate the impact of t on the performance of the saddlepoint approximation by fixing (λ , a, b) at (2, 2, 6). Likewise, to simulate an observation of X t , we simulate t observations of X 1 and then add them together. To plot the kernel density estimate, we simulate 100,000 data points for each value of t. Figure 7 shows the results. Again, it is observed that the saddlepoint density  Figure 7: Density plots of f X t (x;θ θ θ ) (black),f X t (x;θ θ θ ) (red) andf X t (x;θ θ θ ) (green), when t increases.
f X t (x;θ θ θ ) is fairly accurate for each value of t. Moreover, numerical integration off X t (x;θ θ θ ) over the interval (0, +∞) gives the value 1 (due to round-off error). Figure 12 covers more exhaustive parameter settings. Figures 4-7 and 12 verify the competence of the saddlepoint approximation, which shall greatly simplify the inference procedure of the GIG Lévy process.

Parameter Estimation
In this section, we examine the feasibility of directly maximizing ∑ m i=1 log(f X ∆t i (∆x i ;θ θ θ )) for ML estimation. The time series data with size m = 100 are simulated as follows. The time increments {∆t 1 , ∆t 2 , · · · , ∆t m } are randomly sampled with replacement from the set {1, 2, ..., 9, 10}. Then randomly generate ∆t i (i = 1, . . . , m) samples from the GIG distribution GIG(λ , a, b) and set their sum as the realization of the increment ∆x i . A variety of parameter settings are examined: a ∈ {0.5,1,3,5,7,9},b ∈ {0.5,1,3,5,7,0.5,1,3,5,7,9}. Hence,there are in total 288 different parameter settings. For each parameter setting, we repeat for 5000 times and hence obtain 5000 ML estimates of the parameter vector (a, b, λ ). The relative error of every ML estimate is calculated, which is the difference (between the estimate and the true value) divided by the true value.
We remark that, with any one of the three parameters {a, b, λ } being known, the other two parameters can be accurately estimated by directly maximizing ∑ m i=1 log(f X ∆t i (∆x i ;θ θ θ )). We corroborate this statement via the relative-error box plots in Figures 8, 9 and 10. In each column, λ takes in turn a value from {1, 3, 5, 7, 9}. In each row, b takes in turn a value from {1, 3, 5}. In each panel, a takes in turn a value from {0.5, 1, 3, 5, 7, 9}.
• The box plots in Figure 8 characterize the variation of the relative errorâ −a a , assuming that is observed that, for every box plot in Figure 8, the distance between the 1st quantile and 3rd quantile is quite small. For each combination of the values of the three parameters, with 5000 repetitions, the 5000 ML estimates scatter symmetrically at the two sides of the true parameter value, and the median of the 5000 relative errors is virtually zero.
• The box plots in Figure 9 describe the variation of the relative errorb −b b , assuming that the (This is because, as verified in Figure 11, when a is small and λ is large, the value of b has little impact on the GIG density function). However, the median of the 5000 relative errors is still close to zero, implying that the ML estimateb is unbiased.
• The box plots in Figure 10 describe the variation of the relative errorλ −λ λ , assuming that the true value of a is known. For the five panels in each column, a takes in turn a value from {1, 3, 5, 7, 9}. For the three panels in each row, b takes in turn a value from {1, 3, 5}. For the six box plots in each panel, λ takes in turn a value from {0.5, 1, 3, 5, 7, 9}. Again, the 5000 ML estimatesλ scatter symmetrically at the two sides of the true parameter value, the distance between the 1st quantile and 3rd quantile is quite small, and the median of the 5000 relative errors is virtually zero.
Figures 8-10 reveal that, with any one of the three parameters {a, b, λ } being known, the ML estimates of the other two parameters obtained by directly maximizing ∑ m i=1 log(f X ∆t i (∆x i ;θ θ θ )) are unbiased and efficient. Therefore, by fixing any of the parameters {a, b, λ } at an arbitrary value within its domain, we can obtain a new two-parameter Lévy process which, via the saddlepoint technique, can be readily applied to practical problems. In other words, the set of applicable purejump increasing Lévy processes has been significantly enriched.
In the general case, when all the three parameters {a, b, λ } are unknown, the ML estimatê θ θ θ undoubtedly will have a larger variance. Through comprehensive numerical study, we found that, when λ > 1, the ML estimateθ θ θ obtained by directly maximizing ∑ m i=1 log(f X ∆t i (∆x i ;θ θ θ )) is unbiased. However, when λ ≤ 1, the ML estimateλ tends to be larger than the true value, and accordingly the ML estimateb tends to be smaller than the true value. We will tackle this problem through felicitous modifications of the saddlepoint approximation, which is left for future work.

Conclusions
We uncovered the simplicity of the GIG Lévy process by proving that, when λ ≥ −1, the marginal distribution of the GIG Lévy process admits an explicit form, which is a highly accurate approximation. The availability of the analytic and accurate approximation greatly simplifies the problems of parameter estimation, goodness-of-fit testing, random number generation, etc. Particularly, if any one of the three parameters is known, or if λ > 1, the unknown parameters can be accurately and efficiently estimated by directly maximizing the saddlepoint-approximation log-likelihood function. Due to the generality of the GIG Lévy process, the set of practicable pure-jump increasing Lévy processes has been significantly enriched. Our continued work on this process will propose a well-grounded modified saddlepoint approximation.