Point process simulation of generalised hyperbolic Lévy processes

Generalised hyperbolic (GH) processes are a class of stochastic processes that are used to model the dynamics of a wide range of complex systems that exhibit heavy-tailed behavior, including systems in finance, economics, biology, and physics. In this paper, we present novel simulation methods based on subordination with a generalised inverse Gaussian (GIG) process and using a generalised shot-noise representation that involves random thinning of infinite series of decreasing jump sizes. Compared with our previous work on GIG processes, we provide tighter bounds for the construction of rejection sampling ratios, leading to improved acceptance probabilities in simulation. Furthermore, we derive methods for the adaptive determination of the number of points required in the associated random series using concentration inequalities. Residual small jumps are then approximated using an appropriately scaled Brownian motion term with drift. Finally the rejection sampling steps are made significantly more computationally efficient through the use of squeezing functions based on lower and upper bounds on the Lévy density. Experimental results are presented illustrating the strong performance under various parameter settings and comparing the marginal distribution of the GH paths with exact simulations of GH random variates. The new simulation methodology is made available to researchers through the publication of a Python code repository.


Introduction
The statistical properties of randomly evolving phenomena in many real-world applications can be studied using stochastic differential equations.The random behaviour in such systems are typically characterised through a Brownian motion term which implicitly assumes that some version of the Central-limit theorem is valid for the observed random behaviour.However, numerous dynamical systems exhibit more heavy-tailed characteristics than the Gaussian, see for example applications in financial modelling ( [1]- [3]), communications ([4]- [9]), signal processing [10], image analysis ( [11], [12]), audio processing ( [13], [14]), climatological sciences ( [15], [16]), in the medical sciences [17] and for the understanding of sparse modelling/compressive sensing ( [18]- [26]).In such cases the stochastic driving term can be characterised by a Lévy process which generalises the sample paths and the marginal distributions of the term to include other parametric families such as the Poisson process and the α-stable process.
In general a continuous-time randomly evolving system may possess a continuous random Brownian motion component as well as sudden discrete random changes at random times ('jumps').General Lévy processes encompass both of these classes of random evolution such that the increments of the process are independent and stationary ( [27], [28]).Here we focus on a broad class of such processes that are composed purely of jumps.
While there is substantial theoretical and applied interest in simulation of Lévy processes per se, in our work we are ultimately concerned with modelling and inference for systems driven by non-Gaussian Lévy processes, such as the linear stochastic differential equation (SDE) model [29] dX(t) = AX(t)dt + hdW (t) where the more standard Brownian motion is replaced by a non-Gaussian Lévy process {W (t)} (the so-called Background-Driving Lévy process (BDLP) [30]), see for example our earlier work with Stable law Lévy processes ( [31]- [33]).The work presented here though is focussed purely on the generation of the underlying Lévy processes, with the extension of our point process methods to the SDE case being in principle straightforward, as demonstrated in ( [31], [32]) for the Stable law case.A publication in preparation will present convergence results from extending this paper to the linear SDE case.
In this paper, we study simulation methods for a very broad class of Lévy processes, the generalised hyperbolic (GH) process ( [34], [35]) (also known as generalised hyperbolic Lévy motion [36]), which captures various degrees of semi-heavy-or heavy-tailed behaviour such that the tails may be designed to be lighter than non-Gaussian Stable laws (which possess infinite variance [37]), but heavier than a Gaussian [38].Some important special cases include the hyperbolic process [39], the normal inverse-Gaussian process (NIG) [40] and the variance-gamma process [41], which were introduced in the context of modelling empirical financial returns, and the Student-t process (see also [42], [43] which was introduced as an extension to Gaussian processes in Machine Learning, although such processes are not equivalent to the Student-t Lévy processes simulated here).The GH distribution is defined as a normal variancemean mixture where the required mixing distribution is the generalised inverse-Gaussian (GIG) distribution.Our current work improves a point process simulation framework for GIG processes [44], extending it to the variance-mean mixture representation of the GH process, and in addition providing substantial modifications and improvements to the original methods.
The simulation of the sample paths of Lévy processes is a key area of research that enables the use of Lévy processes in inference and decision-making.In [45], Rosiński surveys generalised shot-noise series representations of Lévy processes and their relation with point processes, and this is the general framework adopted for the current paper (see [32], [44], [46], [47] and references therein, for our previous studies using this methodology).Other relevant developments include [48] which present the theory of NIG processes, [49] which present approximate sampling methods for the NIG case, and [50] which give applications of shot -noise series based methods for non-Gaussian Ornstein-Uhlenbeck (OU) processes.Exact simulation methods for the class of tempered stable (TS) processes are studied in ( [51]- [54]).In addition, approximate simulation methods for GH Lévy fields, which are infinite-dimensional GH Lévy processes, are studied in [55].
It is shown in [30] that the GH distribution is infinitely divisible and hence can be the distribution of a Lévy process at time t = 1.The GH distribution possesses a five parameter probability density function defined for random variables on the real line as follows ( [35], [36]) where K ν (•) is the modified Bessel function of the second kind with index ν.The parameter λ ∈ R characterises the tail behaviour, α > 0 determines the shape, 0 ≤ |β| < α controls the skewness, µ ∈ R is a location parameter and δ > 0 is the scale parameter.Alternative parametrisations of the probability density function in the limiting parameter settings are discussed in [35].
The three parameter probability density function f GIG (λ, δ, γ) of the GIG distribution may be linked to Eq. ( 1) via a variance-mean mixture of Gaussians [36].Using the parameterisation γ = α 2 − β 2 , the variance-mean mixture for the GH distribution may be expressed as where u is a GIG distributed random variable.Random variate generation algorithms for a GIG variable are studied in [56], [57], and their extension to GH distributed random variables are then obtained through the normal variancemean construction shown in Eq. (2).
GH processes are generally intractable for simulation since the Lévy density associated with the GIG process is expressed as an integral involving certain Bessel functions.Simulation methods based on generalised shot-noise representation of the GIG Lévy process are given in [44].These methods rely on the construction of dominating point processes that are tractable for simulation, followed by thinning methods derived from upper bounds on the intractable integrand.In the earlier work the series generated must be truncated to a finite number of terms which needs to be tuned by the user, and hence may be inefficient in some parameter regimes.In this paper we show for the first time a practical method for simulation of paths of the GH process, based on subordination with a GIG process.The paper provides several significant contributions.The first contribution is to provide improved simulation methods for the underlying GIG process based on tighter bounds for the construction of dominating processes and the corresponding thinning method, and a proof of convergence is provided for the first time, based on an earlier result by [58].Secondly, we derive adaptive truncation methods for approximating the infinite series involved in our representation which allow for the first time an automatic choice of the truncation level for the jumps of the GIG process.Once truncation has occurred we then approximate the residual error committed by adding an appropriately scaled Brownian motion term with drift, motivated by Central-limit theorem-style results for the residual error.Furthermore, the thinning (rejection sampling) methods are made significantly more computationally efficient through the introduction of 'squeezing' functions that both upper-and lower-bound the acceptance probabilities.Finally, acceptance probability bounds are derived and convergence properties of the novel simulation methods are compared against the methods introduced in [44].The simulation methodology is made available to researchers through the publication of a Python code repository1 .
The paper is organised as follows.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 introduces the specific form of the GIG Lévy density and derives various bounds on these densities as well as constructing dominating Lévy densities from the related bounds.Section 4 gives simulation algorithms for the GH Lévy process based on the simulation of the previously discussed dominating Lévy processes and associated thinning methods.Section 5, presents an adaptive truncation method for the infinite series involved in generalised shot-noise representations and a method of approximating the residual series.Section 6 gives a practical sampling algorithm based on squeezing functions for increasing the efficiency of simulation.Section 7 presents example simulations, comparing the distribution of the paths generated with exact simulations of GH random variates.

Generalised shot-noise representations
In this section series representations of Lévy processes given in ( [45], [59]) that enable their simulation are reviewed.Let W (t) be a Lévy process on some time interval of interest t ∈ [0, T ] having no drift or Brownian motion part, and hence containing purely jumps; then the characteristic function (CF) is given by ( [59], Corollary 13.8), as Under this definition W (T ) is a random variable whose distribution is infinitely divisible.
We will also require a more restricted class of non-negative, non-decreasing Lévy processes X(t), the subordinator process, whose CF is given by: and which has the more restrictive requirement that Q X (dx) defines the density of jumps for {X(t)} such that the expected number of jumps of size In order to generate sample paths from the GH process, we will use the so-called variance-mean mixture representation of its Lévy measure, which is the normal mixture representation of the GH Lévy measure, analogous to the normal mixture representation of its probability density (2), and where Q GIG is the Lévy measure of a Generalised inverse Gaussian (GIG) subordinator process (see [60], [61] and the following section for further detail).
It is first required to simulate a realisation {x i } ∞ i=1 of the jumps from the underlying GIG subordinator process and to use a generalised shot-noise representation [45] to simulate from Q GH : where uniform random variables representing the arrival time of jumps, and independent of the jump sizes {W i } ∞ i=1 , which are independently distributed as: Note that [45] proves the almost sure convergence of such series to W (t) for x i non-increasing, i.e. jumps of X(t) are generated in order of decreasing size.The terms c i are centering terms which we may take as zero for the GH class of processes as a result of the condition in Eq. ( 3).
The task remaining is to generate ordered realisations of the jumps in the subordinator, {x i } ∞ i=1 .Here the Lévy-Itô integral representation of X(t) may be invoked: where N is a bivariate point process with mean measure Leb.× Q on [0, T ] × R 0 which may be represented with Dirac functions as where again {V i ∈ [0, T ]} ∞ i=1 are i.i.d.uniform random variables independent of {X i } which represent the arrival time of jumps, {X i } ∞ i=1 are the jump sizes, and T is the duration of the time interval considered.If we substitute N into Eq.( 6) we obtain a series representation of X(t) as: The classical method to generate such a subordinator process, with Lévy measure Q ( [45], [61]- [63]) simulates jumps of decreasing size by an appropriate transformation of the epochs of a unit rate Poisson process.Briefly, an arbitrarily large number of epochs {Γ i } i=1,2,... is randomly simulated from a unit rate Poisson process.These terms may be transformed into the jump magnitudes of the corresponding subordinator process by calculating the upper tail probability of the Lévy measure A corresponding non-increasing function h(•) is then defined as the inverse tail probability, h(γ) = inf x {x; Q + (x) = γ} which assigns a non-increasing jump value to each of the ordered Poisson epochs, {X i = h(Γ i )}.Thus, small Γ i values correspond to large jumps h(Γ i ) and vice versa.It can be seen from this definition that E[#{X i ; X i ≥ x}] = Q + (x): this procedure is essentially following an analogous formulation to the standard inverse CDF method for random variate generation, but applied here to a point process intensity function instead of a probability distribution.Formally, the mapping theorem [64] ensures that the resulting transformed process is a Poisson process having the correct Lévy density Q(x).
Since there is an infinite number of jumps in the series representation (8), the simulation is in practice truncated at a finite number of terms and the remaining small jumps are ignored or approximated somehow [65], a topic that is addressed in Section 5 of this paper.
The generic method reviewed here requires the explicit evaluation of the inverse tail measure h(γ) which is not tractable for the GIG process.An alternative approach was devised in [44], simulating from a tractable dominating point process N 0 having Lévy measure Q 0 such that dQ 0 (x)/dQ(x) ≥ 1, ∀x ∈ (0, ∞) for which h 0 (γ) is directly available.The resulting samples from N 0 are then thinned with probability dQ(x)/dQ 0 (x) as in ( [45], [66]) to obtain the desired jump magnitudes {x i } of the subordinator process.The generic procedure is given in Alg. 1 for a point process Q(x) having dominating density Q 0 (x) ≥ Q(x) and h 0 (γ) = inf x {x; Q + 0 (x) = γ}.
Algorithm 1 Generation of the jumps of a point process having Lévy density Q(x) and dominating process Q 0 (x) ≥ Q(x).
Generate the epochs of a unit rate Poisson process, {Γ i ; i = 1, 2, 3...}, 3.For i = 1, 2, 3..., Note that our work here will later require partial simulation of such point processes on measurable sets A on jump magnitudes, i.e.Q A (x) = I A (x)Q(x), and typically A will simply be an interval (a, b], b ≤ ∞.This partial simulation is straightforwardly achieved by replacing Step 2) in Alg. 1 with the steps provided in Alg. 2.
Algorithm 2 Generation of Poisson process epochs corresponding to jump magnitudes x i ∈ (a, b] where a > b. and Exp(1) is the unit mean exponential distribution, and noting that the While loop in Alg. 2 may in practice be replaced with a draw from M ∼ P oisson(Q + 0 (b) − Q + 0 (a)) followed by M iid draws for the (unordered) Γ i terms from a uniform distribution U (Q + 0 (b) − Q + 0 (a)).A rejection sampling procedure such as Alg. 1 may be viewed within the generalised shot-noise framework of [58] in which the process is expressed as a random function of the underlying Poisson epochs {Γ i } as follows where H(γ, .) is a non-increasing function of γ, and e i are random variables or vectors drawn independently across i.
[58] Th. 4.1 proves the a.s.convergence of such series under mild conditions.In particular the conditions of Th. 4.1 (A) are satisified.First the distribution of H(., .) is expressed as a probability kernel σ(γ, A) for measurable sets A: and it follows from the Marking Theorem [64] that the resulting point process has Lévy measure Applying this to verify Alg. 1, take H(γ, e) = h(γ)e and e ∈ {0, 1} binomial with P{e = 1} = Q(h 0 (γ))/Q 0 (h 0 (γ)).We will consider only non-zero jumps, since 'jumps' of size zero do not impact the point process, and indeed Lévy measures Q(dx) are not defined for x = 0. Then it follows for measurable sets and hence the resulting Lévy measure is as required.Here we have made the substitution x = h 0 (γ), so γ = Q + 0 (x) and dγ=Q 0 (x)dx.While the procedure of Alg. 1 is well known to be valid, see e.g.[58], we include the sketch proof here since we will use more sophisticated versions of it to prove validity of our own algorithms for GIG and GH process simulation in subsequent sections of the paper.
Simple and well known examples of the procedures in Algs. 1 and 2 are the tempered stable and Gamma processes, which we will require as part our sampling procedures for the GIG process later in the paper.The corresponding Lévy densities and thinning probabilities for these cases are given in [44] (Section 2.1 and 2.2).The associated sampling algorithms are repeated here for reference purposes in Algorithms 3 and 4.
Algorithm 3 Generation of the jumps of a tempered stable process with Lévy density Q T S (x) = Cx −1−α e −βx (x ≥ 0) where 0 < α < 1 is the tail parameter and β ≥ 0 is the tempering parameter.
Algorithm 4 Generation of the jumps of a gamma process with Lévy density Q Ga (x) = Cx −1 e −βx (x ≥ 0) where C > 0 is the shape parameter and β > 0 is the rate parameter.

The generalised inverse Gaussian Lévy process
In this section, the GIG Lévy process and its Lévy measure are defined.Tractable bounds on this Lévy measure are required in order to simulate the GIG (and hence the GH) process, and in this section we provide improved bounds compared with those in [44].These improved bounds are proven in the following section to have higher acceptance rates for the rejection sampling procedures that underlie the algorithms.
The density of the Lévy measure of a GIG process ( [35], Eq. 74), following a change of variables as in [44], is given by is the Bessel function of the third kind, also known as the Hankel function of the first kind, which is defined in terms of J λ (z), the Bessel function of the first kind, and Y λ (z), the Bessel function of the second kind.The presence of an integral involving the Bessel function makes the simulation of such processes intractable except for certain edge cases.
Naturally, the GIG Lévy density can be divided into two terms as and a second term, present only for λ > 0 as which is the Lévy density of a gamma process with shape parameter λ and rate γ 2 /2.It is straightforward to simulate from this second term using Alg.4, thus our attention is directed towards simulation of the point process with Lévy density Q GIG (x).
In order to avoid any direct calculation of the integral in Q GIG (x), the general approach proposed in [44] is to consider a bivariate point process Q GIG (x, z) on (0, ∞) × (0, ∞) which has, by construction, the GIG Lévy density as its marginal, i.e.
Thus, joint samples {x i , z i } are simulated from the point process with intensity function Q GIG (x, z), from which the samples {x i } are retained as samples from Q GIG (x).However, simulation from Q GIG (x, z) is still intractable because of the presence of the Bessel function.This is overcome by constructing tractable bivariate dominating point processes with intensity function Q 0 GIG (x, z) and thinning with probability Q GIG (x, z)/Q 0 GIG (x, z) to yield samples from the desired process Q GIG .
The generic approach proposed here will involve a marginal-conditional factorisation of both point processes: where Q 0 GIG (z|x) and Q GIG (z|x) are proper probability densities, i.e.
Thus z may be interpreted as a marking variable and (x, z) ∈ (0, ∞) × (0∞) form a bivariate Poisson process [64].The generic algorithm for sampling Q GIG (x) is then given below, followed by its proof of validity under the generalised shot noise approach.

Generate Poisson epochs {Γ i } and corresponding ordered jump sizes {x
We now proceed to prove the convergence of Alg. 5 using [58] Th. 4.1 (A).Note that the proof is here presented for the first time.
Lemma 1.The Generalised Shot noise process X(t) = i H(Γ i , e i )I(V i ≤ t) generated according to Alg. 5 converges a.s. to the Poisson point process with Lévy density Q(x).
Proof.Alg. 5 generates, for each point x i , an auxiliary marking variable z i ∼ Q 0 (z i |x i ), and an acceptance variable where once again A 0 are all measurable sets excluding 0. Hence the resulting Lévy measure is as required.The remaining conditions in [58] Th. 4.1 (A) are simply that Q() is a Lévy measure, which is true by construction (it is a subordinator and hence satisfies (3)), and a second technical condition that is always satisfied by subordinators, see [58] Remark 4.1.
A new set of bounds on z|H ν (z)| 2 is now given in Theorem 1 below, which will be used in bounding the overall function (10).The bounds are graphically illustrated for the two distinct parameter ranges in Figs. 1 and 2. The proof of the theorem follows a similar scheme to Theorem 2 in [44] and is hence only briefly stated:

This will define the corner point on a piecewise lower or upper bound. Choose now any
and define the following functions: . .
Proof.As in Theorem 2 of [44], but now we allow for a range of values 0 . The proof is an obvious extension of the original theorem since for any value of z 1 less than , the function A(z) is shifted to the left and lies below of those plotted in Fig. 1 and above of those plotted in Fig. 2, hence providing a valid but less tight bounding function.

Remark. Choice of any z
leads to a poorer (less tight) bound on the true function z|H ν (z)| 2 .In particular, as z 1 → 0 we obtain the crude and well known ( [67], Section 13.75) bound A(z) = 2 π , which was employed in a first version of our method for |λ| > 0.5 ([44], Theorem 1).This asymptotic bound forms the loosest bound A(z) and all other valid choices of z 1 yield increasingly tight bounds on the required Lévy density Q GIG (x, z), which we employ in our new and improved sampling algorithms for |λ| > 0.5.
Corollary 1.1.For any positive z and fixed |λ|, the following bounds are obtained by replacing z|H ν (z)| 2 with A(z) and B(z) in the definition of Q GIG (x, z) (10).For the case |λ| ≥ 0.5 we have: and for 0 < |λ| ≤ 0.5 we have: with equality being achieved in both cases for |λ| = 0.5.Here Q A GIG (x, z) is defined as: and Q B GIG (x, z) defined as: Remark.Setting z 0 = z 1 , it can be clearly seen that the ratio a constant independent of the value of x and z.This fact, which can be clearly visualised in Figs. 1 and 2 (note the log-scale), is utilised later in our development of a retrospective 'squeezed' sampler, see Section 6.
Corollary 1.2.The bound in Eq. (15a) can be rewritten in factorised form as where is a conditional right-truncated square-root gamma density2 with its associated normalising constant.The marginal term Q A N1 (x) is a modified tempered |λ|-stable process3 .Corollary 1.3.The bound in Eq. (15b) can be rewritten in a similar way as is a conditional left-truncated square-root gamma density with its associated normalising constant.The marginal term Q A N2 (x) is a modified tempered 0.5-stable process.
In Corollaries 1.2 and 1.3, the point process intensities correspond marginally to a (modified) tempered stable process in x, and conditionally to a truncated √ Ga density for z.This feature enables sampling from the dominating bivariate point process Q A GIG (x, z) by first sampling x and then, conditional on the value of x, sampling a corresponding z value.Here, for the parameter range |λ| ≥ 0.5, we are extending the approach previously derived only for the parameter range 0 < |λ| < 0.5 [44], (and here summarised in Section 4.3).Notice that allowing z 1 → 0, as discussed in the Remark following Theorem 1, will result in our previous crude bounding function for the parameter range |λ| > 0.5 ([44], Corollary 1).
Since the functions Q A GIG (x, z) and Q B GIG (x, z) form both upper and lower bounds on the target point process Q GIG (x, z), see ( 13) and ( 14), we are able to construct effective sampling algorithms for all parameter ranges |λ| > 0 based upon rejection sampling ideas, see sections 4.1 and 4.3.Furthermore, in Section 6 we use the corresponding lower bounds to create efficient 'squeezed' versions of these algorithms.In the case of 0 < |λ| ≤ 0.5, the new algorithm is a direct development of that presented in [44], while in the case |λ| > 0.5 the algorithm now follows the same structure as for the other parameter range, in contrast with the previous approach from [44] that uses a cruder bound.In all parameter ranges we propose significant improvements over the previous work, including better bounds for simulation of the marginal process, adaptive truncation, simulation of a Gaussian residual term and squeezed rejection sampling.
In the next section we show how to move from simulations of the underlying GIG process towards our ultimate aim here, which is simulation of the GH process.

Simulating GH processes
In this section simulation algorithms for the GH Lévy process are presented.Our approach relies on the definition of a GH process as a subordinated Brownian motion where the subordinator is a GIG process ( [30], [40]).In this approach, jumps sizes {x i } are first generated from the underlying GIG process Q GIG , see previous sections for full details.Then, jumps for the corresponding GH process are obtained as for some β ∈ and σ > 0.
The conditional simulation of the GH process is common to all parameter regimes and is presented in Alg. 6.Given jump times and magnitudes (V i , W i ) the corresponding value of the GH Lévy process at t is given in Eq. ( 5).
Algorithm 6 Simulation of GH process.
1. Generate a large number of points x i from the GIG Lévy process Q GIG , see Algorithms 8 and 10 for |λ| ≥ 0.5, or Algorithms 12 and 14 for 0 < |λ| ≤ 0.5, 2. If λ > 0, draw an additional set of points {x j } from the gamma process λe −xγ 2 /2 x , x > 0, see ( 9) and take the union {x i } ← {x i } ∪ {x j }, 3.For each point x i , draw an independent and identically distributed random variate u i ∼ N (0, 1), 4. The corresponding GH jump sizes w i are obtained as: For each jump w i , draw an independent and identically distributed jump time v i ∼ U(0, T ).
6.The process W (t) is then obtained as: We now detail the methods for generation of the underlying GIG process.
4.1 The case for |λ| ≥ 0.5 Here a new algorithm is presented for simulation in the parameter range |λ| ≥ 0.5, based on the bound Q A GIG (x, z) derived in previous sections, which is an improved bound compared with that in [44] Algorithm 3. The process associated with the Lévy density Q A GIG (x, z) can be considered as a marked point process split into two independent point processes N 1 and N 2 having factorised (marginal-conditional) intensity functions as given in Corollaries 1.2 and 1.3, respectively.
Both N 1 and N 2 correspond to a marginal modified tempered stable process for x and a conditional truncated √ Ga density for z.Each simulated pair (x, z) is accepted with probability equal to the ratio Q GIG (x, z)/Q A GIG (x, z).As a result of the piecewise form of Q A GIG (x, z), the accept/reject steps for N 1 and N 2 may be treated independently and the union of points from the two processes forms the final set of GIG points.The thinning probabilities for points drawn from Q A N1 and Q A N2 are then: Due to the presence of upper and lower incomplete gamma functions in the marginal point process envelopes and direct simulation from these Lévy densities are still intractable and hence dominating processes and associated thinning methods are required.
For Q A N1 (x) the following bound is used to formulate a tractable dominating process ( [68], Theorem 4.1): so that the dominating process can be expressed as: Notice that the point process associated with Q A,d N1 (x) may be considered as the union of two independent gamma processes.Points are then independently accepted with probability . The corresponding algorithm is given in Alg. 7.

Generate a gamma process N 1
Ga having parameters a 1 = z1 2π|λ|(1+|λ|) and For each x ∈ N accept with probability otherwise reject and delete x from N .
Having simulated the x values from the marginal point process associated with Q A N1 (x), the corresponding z values are simulated from a right-truncated square-root gamma density as in Corollary 1.2 and accept-reject steps are carried out to obtain samples from the N 1 point process.The complete procedure is outlined in Alg. 8.
2. Simulate x i from the marginal point process associated with Q A N1 (x) as given in Alg. 7, 3.For each x i , simulate a z i from a truncated square-root gamma density For the simulation of Q A N2 (x) in ( 21), a bound on the term Γ(0.5, z 2 1 x/(2δ 2 )) is established by using the well-known equivalence Γ(0.5, x) = √ π erfc( √ x) where erfc(•) is the complementary error function.Two valid upper bounds on the gamma function are then obtained directly from [69] as While the first inequality is tighter, using it requires the simulation of two TS processes instead of a single process.Hence in the current implementation the right hand bound in Eq. ( 24) is chosen and the associated dominating point process envelope is then given by which can be simulated as a TS process and for each x i the probability of acceptance is Γ(0.5, ).The corresponding algorithm is given in Alg. 9.
1. Generate a tempered stable process N M T S with parameters C = δ √ 2π , α = 0.5 and β = 2. For each point x ∈ N M T S , accept with probability Γ(0.5, ) ), otherwise reject and delete x from N M T S .
Using the simulated values x i , the corresponding z i values are generated from the conditional left-truncated square-root gamma density Q A N2 (z|x) and the whole procedure is outlined in Alg.10.
1. N 2 = ∅, 2. Simulate x i from the marginal point process associated with Q A N2 (x) as shown in Algorithm 9, 3.For each x i , simulate a z i from a truncated square-root gamma density Γ(0.5) √ Ga(z|0.5, x i /(2δ 2 )) Γ(0.5, z 2 1 x i /(2δ 2 )) Note that the bound shown in Eq. ( 24) is a significant improvement over the bound based on the complete gamma function as used in Alg.7 of [44].
Finally, the set of points N = N 1 ∪ N 2 is a realisation of jump magnitudes corresponding to a GIG process having intensity function Q GIG (x).The associated GH process may be obtained using Alg.6.
Remark.Note that whenever λ > 0, the set of points generated from the process with intensity function in Eq. ( 9) is added (by taking a union operation) to the set of points coming from Q GIG (x) in Algorithms 8 and 10 for |λ| ≥ 0.5, or Algorithms 12 and 14 for 0 < |λ| ≤ 0.5, to obtain the full set of jumps {x i } from the required GIG process.
Remark.Notice that for γ = 0, the gamma process N 1 Ga in Alg.7 is not well-defined since its rate parameter is zero and hence N 1 cannot be simulated.In this case, set z 1 = 0 and then only samples from N 2 are required.In this case the conditional density for z becomes the complete square-root gamma density instead of a truncated one and the marginal modified tempered stable process for x reduces to a standard TS process that does not require any further thinning.The simulation algorithm then becomes equivalent to Alg. 3 in [44].For all other values of γ we , in order to achieve the improved bound according to Theorem 1.

Acceptance rates for simulation from
We now analyse the acceptance rates for the new procedure.This will enable a quantitative comparison with our previous methods in [44].The acceptance probabilities for the two point processes N 1 and N 2 are obtained from ( 18) and ( 19) as The expected value of the acceptance rates for fixed x may be evaluated w.r.t. the sampling densities for z, i.e.Q A N1 (z|x) (20) and Q A N2 (z|x) (21): dz However, the presence of the term z|H |λ| (z)| 2 in both integrals makes them intractable.The expected acceptance rates may then be bounded by using the same functions A(z) and B(z), introduced in Theorem 1, to replace z|H |λ| (z)| 2 .Note that only the lower bound on the acceptance rates are of interest here since upper bounding both expectations using A(z) leads to a trivial upper bound of 1 on the acceptance rates.
The acceptance rates associated with the N 1 and N 2 processes can be lower bounded using Theorem 2 and the proof is provided in the Appendix.
For any fixed x and |λ| ≥ 0.5, the following lower bounds on E [ρ 1 (x, z)] and E [ρ 2 (x, z)] apply: |λ|−0.5 γ(0.5, Γ(0.5, Note that the corner point z 0 ∈ (0, ∞) may be chosen arbitrarily, while z 1 is considered to be fixed.Therefore the lower bound may be optimised for both N 1 and N 2 w.r.t.z 0 for each x value.For N 1 the optimised lower bound may be expressed as |λ|−0.5 γ(0.5, , z 0 ∈ (0, z 1 ) Similarly for N 2 the optimised lower bound may be expressed as Note that the acceptance rates for our previous algorithm Section 3.1.1 of [44] can be considered as a limiting case of the new procedure corresponding to z 1 = 0, and hence the new acceptance rates are at least as large as the previous rates for each x > 0. In Fig. 3, the optimised lower and upper bounds on the acceptance rates for the new procedure are shown.Additionally for comparison, the lower bounds on the previous algorithm (z 1 = 0) are plotted in Fig. 3.This illustrates that the new procedure is a significant improvement in terms of the acceptance rate for N 2 (i.e.E[ρ 2 (x)]) for small |λ| values, and a slight improvement for larger |λ| values.

4.3
The case of 0 < |λ| ≤ 0.5 This parameter range was covered in [44] and we review the basics here for completeness.We provide several improvements to this approach, including the more efficient sampling of point process N 1 as two gamma processes, improved bounds on the incomplete gamma functions, as well as the adaptive truncation method, residual approximation and the squeezed sampling methods detailed in subsequent sections.
The process associated with the upper bounding Lévy density Q B GIG (x, z) for this parameter range, see (16a) and (16b), can be considered once again as a marked point process split into two independent point processes N 1 and N 2 having factorised intensity functions as given in Corollary 2 of [44] as Again, N 1 and N 2 correspond to a marginal modified tempered stable process for x and a conditional truncated √ Ga density for z.The upper and lower incomplete gamma functions in the marginal point process envelopes Q B N1 (x) and Q B N2 (x) require the use of dominating processes and thinning methods similar to section 4.1.Using the bound in Eq. ( 22) the density Q B N1 (x) can be transformed into two independent gamma processes and the methodology is summarised in Alg.11.This simulation algorithm improves upon the method shown in Algorithm 4 of [44] by transforming the problem of simulating a tempered stable process into that of simulating two gamma processes, which are known to converge rapidly in terms of the number of points required in Eq. ( 5).The convergence of these sums are discussed further in Section 5. Having simulated points x i from the marginal density Q B N1 (x), z i values are simulated from a right-truncated square-root gamma density and the accept-reject step for each x i is performed as shown in Alg.12.

Generate a gamma process N 2
Ga having parameters a 2 = z0 π 2 H0(1+|λ|) and otherwise reject and delete x from N .
Algorithm 12 Generation of N 1 for 0 < |λ| ≤ 0.5 2. Simulate x i from the marginal point process associated with Q B N1 (x) as shown in Algorithm 11, 3.For each x i , simulate a z i from a truncated square-root gamma density For Q B N2 (x), the incomplete gamma function is upper bounded by the complete gamma function to produce a process that is tractable for simulation and the corresponding algorithm is given in Alg. 13.For each x i value simulated from the marginal density, a corresponding z value is sampled from the conditional left-truncated squareroot gamma density and the whole methodology is outlined in Alg.14.
Note that the bound defined in Eq. ( 24) could also be used in this parameter regime to obtain a more efficient sampling algorithm; however, the acceptance rates using the simpler bound are found to work well.
Algorithm 14 Generation of N 2 for 0 < |λ| ≤ 0.5 1. N 2 = ∅, 2. Simulate x i from the marginal point process associated with Q B N2 (x) as shown in Algorithm 13, 3.For each x i , simulate a z i from a truncated square-root gamma density Γ(0.5) √ Ga(z|0.5, x i /(2δ 2 )) Γ(0.5, z 2 0 x i /(2δ 2 )) Lastly, the set of points N = N 1 ∪ N 2 is a realisation of jump magnitudes corresponding to a GIG process having intensity function Q GIG (x) and once again the corresponding GH process may be obtained using Alg.6.

Adaptive truncation and Gaussian approximation of residuals
The shot noise methods discussed in previous sections involve infinite series of decreasing random variables and in practice must be truncated after a finite number of terms.In this section we propose novel methods for adaptive determination of the number of terms required in the truncated series.This adaptive truncation can both save substantially on the computational burden of generating very long shot noise series, and also ensure (probabilistically) a specified error tolerance.Furthermore, we provide lower bounds on the mean and variance of the GIG and GH residual sequences which will be used in order to approximate the residual term as Brownian motion once the adaptive truncation procedure is terminated.
The adaptive truncation and residual approximation methods studied in this section are designed to match the moments of a realisation from a GH process at time t to its theoretical moments.For a subordinator Lévy process X(t) with Lévy measure Q(dx) and finite first and second moments, the mean and variance of the subordinator may be expressed as [70] E and where the distribution of the process at time T , denoted as X(T ), defines the associated random variable.These integrals are intractable for the GIG case.Hence upper and lower bounds on these integrals are studied.
For normal variance-mean processes W (t), such as the GH process, the associated Lévy measure can be expressed as a function of the subordinator Lévy measure Q GIG (dx) as in (4).Hence the mean and variance of the GH process can be found in terms of the moments of a GIG process as and where β and σ are the skewness and scale parameters as defined in (17).
In Section 5.1 we study adaptive truncation of infinite series for subordinator processes and provide associated algorithms for the GIG case.In principle the adaptive truncation scheme can also be described in terms of the GH moments using (32) and (33).This transformation to the GH moments is explicitly required for the Gaussian approximation of residual moments of the GH process and lower bounds on these residuals are studied in Section 5.2.The Gaussian approximation we use is motivated by proofs of the convergence of normal variance-mean mixture residuals to a Brownian motion in all cases of the GH process except for the normal-gamma process, and these results will be presented in a forthcoming publication.

Adaptive truncation of shot noise series
The shot noise series for a subordinator X(t) with its jumps truncated at ε may be defined as The difference of X(t) and the truncated series X ε (t) characterises the residual error caused by truncation and may be expressed as a random process R ε (t) such that where ε is the value at which the jump magnitude sequence {x i } are stochastically truncated (i.e. the truncated series X ε has a random number of terms with x i greater than or equal to ε).
The statistical properties of the residual error R ε (t) as a function of ε can be used to study the convergence of the truncated series to the Lévy process X(t).The number of terms used in the approximation of X(t) may be dynamically adjusted depending on the particular realisations of {x i } and the required precision of approximation.Theorem 3 below and its Corollary describes the construction of a probabilistic bound on the residual error caused by truncation of a subordinator process in terms of upper bounds on its residual moments and provide the residual mean and variance for the tempered stable and gamma processes which are used as dominating processes for sampling the GIG process.
Note that the bound in the Theorem below is a pointwise bound at a particular time t, whereas ideally a pathwise bound might be desired that applies across all times.Martingale inequalities can be used in principle to achieve this, see [71], although these may not be directly applicable here as we do not in general have an exact characterisation of the residual mean and variance for the GIG process, only upper and lower bounds on these.Theorem 3.For the residual error R ε (t) associated with truncation of a subordinator process, the following probabilistic bound applies for any E > με and truncation level ε > 0: where με ≥ µ ε and σε ≥ σ ε are upper bounds on , and E is a threshold that may depend on the random realisation X ε (t).
Proof.The mean and variance of a subordinator process X(t) are given in (30) and (31).Similarly, the mean and variance of the truncated process residual R ε (T ) can be found as and where both of these integrals are well-defined and finite for any valid subordinator and 0 < ε < ∞, by condition (3).Now, using the expected value and standard deviation of R ε (t) we may bound the residual error using concentration inequalities.Specifically, Chebyshev's inequality states that for a random variable R ε (t) with finite expected value µ ε and finite non-zero variance We require here only the right tail probability mass corresponding to the event R ε (t) − µ ε ≥ kσ ε , and this is clearly less than or equal to the probability of the event |R ε (t) − µ ε | ≥ kσ ε ; hence rearranging we arrive at Now, if we have instead upper bounds µ ε ≤ με and σ ε ≤ σε , it is clear that µ ε + kσ ε ≤ με + kσ ε and so Finally, a simple rearrangement with E = με + kσ ε leads to the theorem as stated.
Corollary 3.1.Our simulation algorithms for the GIG process involve thinning/rejection sampling operations in order to generate point processes N 1 and N 2 from gamma and tempered stable dominating processes, see Algs. 7, 9, 11 and 13.By construction, the resulting thinned processes have Lévy density Q(x) strictly less than or equal to that of the dominating process in each case, say Q 0 (x), i.e.Q(x) ≤ Q 0 (x).Hence the means and variances of the truncation error, calculated using ( 37) and (38), are strictly less than or equal to those of the corresponding dominating gamma and tempered stable processes, and we may thus take the means and variances of the underlying TS or gamma processes as the upper bounds με and σε required in Theorem 3.
For the TS process the expected value and variance of the residual process is found as and In limit as β → 0 we obtain the stable subordinator whose moments can be obtained either directly or as limits of the above TS case using γ(s,x) x s → 1 s as x → 0, giving: and Similarly for the gamma process the expected value and variance of the residual process can be found as and Take, for example, generation of N 2 in Algs.9 and 10.The starting point is generation of a TS process with parameters C = δ √ 2π , α = 0.5 and β = 2 , implemented using Algorithm 3. The mean and variance for the truncated residual of this process are obtained from ( 40) and (41).Algs.9 and 10 then perform random thinning on the TS points.Hence the resulting process N 2 has truncated residual with mean and variance no larger than those of the corresponding TS process.Thus Theorem 3 applies, using the TS mean and variance as the upper bounds με and σ2 ε .The other cases of N 1 and N 2 simulation follow a similar argument, using the appropriate gamma or TS process to generate upper bounds on the moments required for Theorem 3.
Finally, a probabilistic upper bound on the GIG residual may be obtained by adding the upper bounds on the means and variances for N 1 and N 2 , since the two point processes are independent.Corollary 3.2.Improved residual errors and corresponding bounds on these are available if the mean µ is available, as then the improved estimate X (t) = X (t) + µ may be formed as proposed in [65].In our GH case however we only have upper and lower bounds µ ≤ µ ≤ μ established in Theorems 3 and 4 (see below).In this case it would seem appropriate to take a conservative line and substitute the lower bound µ ε in place of µ .Then Theorem 3 can be modified in step (39) as follows, To justify this, use Chebyshev directly to give But we have and hence Pr(A) ≤ Pr(B) from which (46) follows.
Rearranging this expression with E = με − µ ε + kσ ε a new expression is obtained, valid for E + µ ε − με > 0: This expression would then be recommended for practical use in adaptive truncation schemes, yielding always smaller probabilities of exceedance than Theorem 3 for fixed threshold E since E + µ − με ≥ E − με .
Using the probabilistic bounds given in Theorem 3 and its Corollaries, an adaptive truncation scheme can be devised to determine a suitable value for ε for each generated realisation of the process.As ε decreases, so we accumulate sequentially the realised value of X ε (t) (or its mean-adjusted version from Corollary 3.2) according to Eq. (34).A tolerance E = τ X ε (t) is chosen, where 0 < τ 1, which is designed to truncate the series once the predicted residual has become very small in comparison with the series realised to level ε.Then a probability threshold, p T 1 is chosen for comparison with P r(R ε (t) ≥ E), in order to decide when to terminate the simulation.A generic adaptive truncation scheme is outlined in Alg. 15 for a point process N associated with a subordinator Lévy process X(t) having Lévy density Q(x).
Algorithm 15 Simulation of N with adaptive truncation using tolerance τ and probability threshold p T for a point process with Lévy density Q() having moment bounds μεn and σεn .

Define a decreasing truncation schedule ε
• Else B = True.
x > 0, having moment bounds μk εn and σk εn with adaptive truncation using tolerance τ and probability threshold p T .1. Define a decreasing truncation schedule 2. While B = False The GH simulation algorithms studied in this work and [44] are made up of two independent point processes N 1 and N 2 .An adaptive truncation algorithm such as Alg.15 can be applied separately to each process to obtain the resulting jumps.This misses a trick however, since either series could in principle be truncated even earlier once its residual error is very small relative to the accumulated sum of both N 1 and N 2 .There are many workable schemes based around this idea and one possible such approach is presented in Alg.16.It is presented in a general form that can apply to the parallel simulation and adaptive truncation of K independent subordinator point processes N k having Lévy densities Q k (x) and overall Lévy density For N 1 in the most general settings of Alg. 8 and 12, the dominating process is made up of two independent gamma processes N 1 Ga and N 2 Ga .In the case of N 2 in both settings, Alg. 10 and 14, a single dominating tempered stable process is required.Hence an efficient method of simulation is running Alg.16 on these 3 independent dominating processes.The residual means and variances of the tempered stable and gamma processes required by Alg. 15 and Alg.16 are shown in Corollary 3.1.It is worth noting that the convergence of a gamma process is typically significantly faster than a tempered stable process and so the N 1 process tends to terminate much sooner than N 2 .
Remark.In the edge parameter setting λ < 0 and γ = 0, the marginal point process simulation methods shown in Alg.7 and 11 are not valid as a result of N 1 Ga becoming undefined for γ = 0.For this setting, Alg. 3 of [44] may be used together with the adaptive truncation and residual approximation methods introduced in this section.
For this parameter setting the tempered stable process defined in Alg. 3 of [44] becomes a stable process since the tempering parameter β is equal to 0. The associated residual moments of a stable process are presented in Corollary 3.1 which should be used to implement Step 2) of Alg. 15.

Gaussian approximation of residual errors
Here we present a Brownian motion approximation method for the residual error R ε (t) of a GIG or GH process caused by the truncation of a shot noise series as defined in Eq. (35).Such an approach is well known from previous work, see e.g.[65], but here we propose an intermediate solution in which a Brownian motion is injected whose drift and variance are lower bounds compared with the exact result, which is intractable in general for the GIG and GH processes.
The theoretical mean and variance of the residual error for the GH case can be found as a function of the mean and variance of an associated GIG residual error R ε (t) using Eqs.( 32) and (33).Hence similar to Section 5.1, we provide lower bounds µ ε , σ 2 ε on the mean and variance of a residual GIG process as a function of the truncation level ε.Together with the upper bounds discussed in Corollary 3.1, these lower bounds characterise the residual error of truncating the infinite shot noise series for the GIG and GH cases.Using the residual approximation module the series representation of the GIG Lévy process X(t) can be expressed as where X ε (t) is computed in the usual way as {i:xi≥ε} x i I Vi≤t , B(t) is an independent standard Brownian motion term and µ ε , σ 2 ε are lower bounds on the mean and variance of the residual error R ε (T ) given a truncation level ε.According to Eqs. ( 32) and ( 33), the lower bounds on the moments of the residual error R ε W (t) of the associated GH process W (t) can be obtained as Hence for the GH process the same procedure is adopted to obtain an approximation as The approximation of the residual error is in addition to the adaptive truncation methods described in Section 5.1 that provide a specified level of truncation ε for the simulation of sample paths from a GH process.In Theorem 4 below we provide the required lower bounds on the mean and variance of a residual GIG process with Lévy density Q GIG (x) as a function of ε.The derivation of these bounds are presented for a specific T as it is straightforward to scale the moments according to time.The lower bounds are then used to evaluate the mean and variance of the residual error in the GH process simulation and approximate the contribution of the residual small jumps according to Eq. (49).A Brownian motion approximation to the residual is known to hold for many shot noise series, see [65], with the gamma process being a well-known exception that does not converge to a Gaussian as ε → 0. In our own work we have proven convergence of the shot noise series for the GH process to a Brownian motion in all cases except the normal-gamma, and these results will be presented in a future publication.
Note that the lower bounds presented in Theorem 4 are valid only for the λ < 0 setting.The additional mean and variance components present for λ > 0 are associated with a gamma process with Lévy density (9) and residual moments given by Eqs. ( 44) and (45).These terms are simply added to the lower bounds µ ε and σ 2 ε in order to obtain the final bounds on the GIG process.
Remark.For the special case of |λ| = 0.5, the GIG density Q GIG (x) has a functional form equivalent to a TS process with intensity function Thus both the simulation algorithms in Section 4, and the adaptive truncation and residual approximation methods presented in this Section are significantly simplified.It is straightforward to simulate a TS process using Alg. 3 and the residual moments shown in Eqs. ( 40) and ( 41) provide exact expressions for the residual moments of the equivalent GIG process.
The true mean and variance of the residual moments replace the upper and lower bounds required for determining the adaptive truncation level ε and the associated moments of the approximating Brownian term.In this case Eq. ( 49) exactly matches the first and second moments of the residual process.

Squeezed rejection sampling
In this section, we provide a practical extension to the sampling algorithms discussed in Section 4 that is designed to increase the efficiency of simulation.The above methods for sampling of N where we have used the squeezing inequality (13) and where the final equality applies when we set z 0 = z 1 , see Remark 3.This implies that we may carry out a retrospective sampling step: draw a uniform random variate W i on [0, 1], test whether it is less than or equal to 2 πH0 , and if so, accept x i with no Steps 3 and 4 required.If W i > 2 πH0 , carry out steps 3 and 4, using the same realised W i to carry out the test in Step 4 5 .The modified version of Alg. 8 using squeezed rejection sampling is presented in Alg. 17.
An exactly similar modification applies for generation of N 2 in Alg. 10 for |λ| > 0.5.In the other parameter range |λ| < 0.5, the bounds are reversed, see (14) and so Step 4 in the squeezed sampler is replaced with 'if w i ≤ πH0 2 '; otherwise squeezed procedures for N 1 and N 2 in this parameter range are modified exactly as in Alg. 17.
We can see that the savings arising from this method could be substantial in cases where 2 πH0 is close to unity, saving almost all of the heavy computations in the original Steps 3 and 4.This will occur when |λ| is close to 0.5, and improvements will lessen as |λ| moves away from this value.The fraction of saved computation is shown as a function of |λ| are shown in Fig. 4, exhibiting a broad range of λ values for which savings are useful.

Simulations
In this section realised sample paths from the simulation algorithms in Section 4 together with the modifications discussed in Section 5 are presented.The distribution of the realised GH processes generated up to T = 1 is compared against samples from the GH distribution which are generated using exact GIG random variable samplers ( [56], [73]).While these methods are able to generate samples for a specific t = T , our method is able to generate the entire path of the process up to t = T and hence information about the dynamics within the interval (0, T ) is made available.Furthermore, the accuracy and efficiency of the proposed simulation algorithm and modifications are compared with the method presented in [44].
Since there is no known closed-form cumulative distribution function (CDF) for the GH distribution, the accuracy of the simulated sample paths are measured using the two-sample Kolmogorov-Smirnov (KS) test which compares the empirical distribution functions of the sample paths at T = 1 and samples from the GH distribution.These tests involve 10 6 independent samples from each method and the results are shown in Tables 1, 2 and 3 for various parameter settings.The probability threshold p T = 0.05 for the adaptive truncation algorithm is found to work well and is used for all cases throughout the section.In Table 1 the results are shown for the simulation procedure using only adaptive determination of the number of jumps and in Table 2 the improvements in convergence by approximating the residual small jumps, as studied in 5.2, are shown.There is a clear trade-off between the accuracy of the distribution at T = 1 and the time required per sample for a given λ value.This can be observed by comparing the KS statistic and time per sample for different tolerance parameters τ , e.g. in Alg.16.A large tolerance results in worse convergence but allows faster simulation, and hence the tolerance parameter may be adjusted depending on the requirements of the application.Furthermore as |λ| increases both the time required per sample and the KS statistic increases as a result of the reduced acceptance rates of the simulation algorithm as plotted in Fig. 3.
Comparing Tables 1 and 2, it can be seen that there is no significant increase in time per sample when applying residual approximation methods and there is an improvement in the KS statistic for all parameter settings.Particularly, the results suggest that the improvement in large tolerance parameter cases where τ = 0.1 are substantial.In order to compare our new methods with the method studied in [44], the adaptive truncation algorithm in Alg.16 is limited to producing a certain maximum number, M = 10 4 , of jump magnitudes per sample path while running the experiments for Tables 1 and 2. The KS statistics and the time required per sample for the same parameter settings and a fixed M number of jumps from the method in [44] are reviewed in Table 3. Comparing these results to Table 2, it can be seen that improved convergence results can be obtained by the novel algorithms introduced in this work while also providing significantly reduced time complexity.
To present a more intuitive understanding of the accuracy of our methods, QQ (quantile-quantile) plots of sample paths generated up to T = 1 and samples from a random variable generator are shown in addition to histograms with the true probability density function overlaid.The parameter values in the examples are selected to reflect the different characteristic behaviour of the GH process as well as edge cases such as the normal-inverse Gaussian process (λ = −1/2), and the Student-t process (γ = 0, λ ≤ 0).The normal-inverse Gaussian (NIG) distribution, or the distribution of NIG processes at T = 1, forms an exponential family and hence all of its moments have analytical expressions [60].As a result of its tractable probabilistic properties, the NIG process finds application in modelling turbulence and financial data ( [48], [49]).The NIG process is in fact a special case of the GH process where λ = −0.5.This results in the bounds given in Corollary 15b being exactly equal to the GIG density and thus the acceptance rate of points simulated from Alg. 3 of [44] is 1.0.The QQ plot, density estimate and sample paths for this parameter setting are shown in Figs. 13 and 14.Another special case of the GH process is the Student-t process where λ < 0, γ = 0 and δ 2 = −2λ.The Student-t distribution is parameterised using a single parameter ν, called the degrees of freedom, which is a positive real number.This parameter is related to the usual parameters of a GH process such that λ = −ν/2 and δ = √ ν.
As remarked in Section 5, it is not possible to simulate a Student-t process using Algorithms 7, 8, 9 and 10 because the corresponding gamma processes are not well-defined in this case.Instead, for this parameter setting Algorithm 3 of [44] is used together with the adaptive truncation and residual approximation methods presented in Section 5 to produce the jumps from the subordinator GIG process.The QQ plot, density estimate and sample paths for the resulting samples from the Student t process are shown in Figs. 15 and 16.Note that removing the condition δ 2 = −2λ still results in a well defined Lévy process with its marginal distribution parameterised by Eq. (3.11) in [35].It is common to apply the Student-t distribution to financial data sets as an alternative to the Gaussian distribution to account for the heavier tails observed in asset returns.There has been numerous studies that suggest that the asymmetric Student-t distribution, β = 0, provides a better fit to financial data sets compared to its symmetric counterpart ([74]- [76]).To the best of our knowledge, we present the sample paths of an asymmetric Student-t process for the first time in Fig. 17 together with the resulting marginal QQ plot and density estimate in Fig. 18.The marginal density of this limiting case of asymmetric GH processes are given in Eq. (3.9) of [35].

Conclusions
The point process representation of a generalised hyperbolic process and the generalised shot noise methods developed in this work provide the the first complete methodology for simulation of generalised hyperbolic (GH) Lévy processes, giving a unified framework for a very broad range of heavy-tailed and semi-heavy-tailed non-Gaussian processes.The continuous time formulation, simulating directly in continuous time path space, can be employed for accurate uncertainty propagation, path visualisation, modelling and inference, especially for irregularly sampled time series datasets.
The presented methods are based on the subordination of a Brownian motion by the generalised inverse Gaussian (GIG) process, and we have here provided novel improvements in GIG process simulation compared with our previous work [44].We have proved also that these series representations are almost surely convergent, verifying the conditions presented in [58].In addition to these improvements we present a novel scheme for adaptive truncation of the random shot noise representation based on probabilistic exceedance bounds, relying on new upper and lower bound expressions for the moments of truncated residual of the shot noise process; these truncations methods are shown to reduce computational burden dramatically without noticeable compromising of accuracy.Further computational savings are made through the use of squeezed rejection sampling, again based on our lower and upper moment bounds.We plan to use the new GH process simulators as a fundamental building block for modelling of stochastic differential equations (SDEs) driven by GH processes, in a spirit similar to [32], where the conditionally Gaussian form of our models will be of great benefit in inference for states and parameters for GH process-driven SDEs, finding application in spatial tracking, finance and vibration data modelling, to list only a few possibilities.These models and their applications will be further studied in future publications.where β 0 > 1 and x ≥ 0. The free parameter β 0 can be adaptively optimised across the whole domain.The bound can be derived using the well-known equivalence Γ(0.5, x) = √ π erfc( √ x).
Hence a lower bound on Q B N2 is obtained as which is the density of a tempered stable process and the associated residual moments are previously derived in the proof of Theorem 3 and shown in Eq. ( 40) and (41).The lower bound on the residual mean and variance of Q B N2 for a truncation level ε are then obtained directly as For the case 0 < |λ| ≤ 0.5, the marginal Q GIG (x) is bounded as Similarly for n ∈ {1, 2}, lower bounds on the moments may be found as The lower bounds on the incomplete gamma functions established for the previous parameter setting can be directly substituted into Q A N1 (x) and Q A N2 .The resulting densities are again characterised as a gamma process and a tempered stable process respectively such that Hence the associated residual moments of these processes can be directly found again using Eqs.( 44), ( 45) and ( 40), (41).The corresponding lower bound on the residual mean and variance of the GIG process for the case 0 < |λ| ≤ 0.5 and truncation level ε are then found as and the number of jumps is a Poisson random variable with mean µ [a,b] .We will be dealing with infinite activity processes for which ∞ 0 Q X (dx) → ∞ and hence there are almost surely an infinite number of jumps in time interval [0, T ].

Figure 4 : 3 .
Figure 4: Fractional computational saving for the retrospective sampling procedure, as a function of |λ|.

For
the case |λ| > 0.5, Figs. 8, 10, 12 shows the QQ plot and histogram for different parameter settings and Figs.7,9,11 presents sample paths from our new adaptive simulation algorithm with residual approximation.The corresponding processes are simulated using Algorithms 7, 8, 9 and 10.Similarly, Figs.5,6 presents an example in the 0 < |λ| < 0.5 setting and the corresponding simulation methods are described in Algorithms 11, 12, 13 and 14.The number of samples from each simulation method is N = 10 6 and the sample path plots show randomly selected 50 paths.

Figure 6 :
Figure 6: Simulation comparison between the shot noise generated GH process and GH random variates, λ = −0.4,γ = 0.1, δ = 1, β = 0.The adaptive truncation parameters are p T = 0.05 and τ = 0.01.Left hand panel: QQ plot comparing our shot noise method (y-axis) with random samples of the GH density generated using a random variate generator (x-axis).Right hand panel: Normalised histogram density estimate for our method compared with the true GH density function.

Figure 8 :
Figure 8: Simulation comparison between the shot noise generated GH process and GH random variates, λ = −0.8,γ = 0.1, δ = 1, β = 0.The adaptive truncation parameters are p T = 0.05 and τ = 0.01.Left hand panel: QQ plot comparing our shot noise method (y-axis) with random samples of the GH density generated using a random variate generator (x-axis).Right hand panel: Normalised histogram density estimate for our method compared with the true GH density function.

Figure 12 :
Figure 12: Simulation comparison between the shot noise generated GH process and GH random variates, λ = −10, γ = 0.1, δ = 1, β = 0.The adaptive truncation parameters are p T = 0.05 and τ = 0.1.Left hand panel: QQ plot comparing our shot noise method (y-axis) with random samples of the GH density generated using a random variate generator (x-axis).Right hand panel: Normalised histogram density estimate for our method compared with the true GH density function.

Figure 13 :
Figure 13: Simulation comparison between the shot noise generated NIG process and NIG random variates, λ = −0.5, γ = 0.1, δ = 1, β = 0.The adaptive truncation parameters are p T = 0.05 and τ = 0.01.Left hand panel: QQ plot comparing our shot noise method (y-axis) with random samples of the NIG density generated using a random variate generator (x-axis).Right hand panel: Normalised histogram density estimate for our method compared with the true NIG density function.

Figure 15 :
Figure 15: Simulation comparison between the shot noise generated student-t process and student-t random variates, λ = −2.5, γ = 0, δ = √ 5, β = 0.The adaptive truncation parameters are p T = 0.05 and τ = 0.1.Left hand panel: QQ plot comparing our shot noise method (y-axis) with random samples of the student-t density generated using a random variate generator (x-axis).Right hand panel: Normalised histogram density estimate for our method compared with the true student-t density function.

2
The density Q B N2 is a modified tempered stable characterised by the upper incomplete gamma function, which can be lower bounded as ([78], Th. 2)

2 Finally
the lower bound on the residual mean and variance of the GIG density can be expressed in terms of the residual means and variances of the lower bounding gamma and TS processes as µ ε = µ B N1 (ε) + µ B N2 ( 1and N 2 in both parameter settings, given in Algs.8,10 and 12,14, involve a computationally expensive pair of steps in the sampling of a truncated gamma random variate (Step 3) and a pointwise evaluation of the Hankel function (Step 4).Notice however, that Theorem 1 provides both lower and upper bounds on the term z|H ν (z)| 2 , which allows us to specify squeezing functions on Q GIG (x, z),4as given in (15a)-(16b).This allows for a labour-saving retrospective sampling procedure in which, for a fixed fraction of points x i , we may replace the simulation of a conditional random variable z and rejection sampling based on its value (Steps 3 and 4) with a simple 1-step accept/reject and no requirement to sample z or evaluate H |λ| .Considering first the case |λ| ≥ 0.5 and the sampling of process N 1 , see Alg. 8, we have generated at Step 2 a single point realisation x i from the process Q A N1 (x).Now consider Step 4, which accepts x i with probability

Table 1 :
The results of two-sample KS tests for adaptive simulation algorithm.

Table 2 :
The results of two-sample KS tests for adaptive simulation with residual approximation algorithm.

Table 3 :
[44]results of two-sample KS tests for the simulation algorithm in[44].