Fast simulation of tempered stable Ornstein–Uhlenbeck processes

Constructing Lévy-driven Ornstein–Uhlenbeck processes is a task closely related to the notion of self-decomposability. In particular, their transition laws are linked to the properties of what will be hereafter called the a-remainder of their self-decomposable stationary laws. In the present study we fully characterize the Lévy triplet of these a\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a$$\end{document}-remainders and we provide a general framework to deduce the transition laws of the finite variation Ornstein–Uhlenbeck processes associated with tempered stable distributions. We focus finally on the subclass of the exponentially-modulated tempered stable laws and we derive the algorithms for an exact generation of the skeleton of Ornstein–Uhlenbeck processes related to such distributions, with the further advantage of adopting procedures which are tens of times faster than those already available in the existing literature.


Introduction
The Lévy-driven Ornstein-Uhlenbeck (OU) processes have attracted considerable interest in recent studies because of their potential applications to a wide range of fields. As observed in Barndorff-Nielsen and Shephard (2001), these OU process are mathematically tractable and can be seen as the continuous-time analogues of the autoregressive AR(1) processes (see Wolfe 1982). They constitute indeed a rich and flexible class that can accommodate features such as jumps, semi-heavy tails and asymmetry which are well evident in the real physical phenomena as well as in the financial data. The energy and commodity markets exhibit for instance a strong meanreversion and sudden spikes which makes the use of the Lévy-driven OU-processes more advisable than the standard Gaussian framework. In addition, several approaches based on Lévy processes such as the Variance Gamma (VG) or the Normal Inverse Gaussian (NIG) have been proposed to overcome the known limits of the usual Black-Scholes model (see Madan andSeneta 1990 andBarndorff-Nielsen 1998): all these non Gaussian noises can of course be adopted as the drivers of OU processes. Example of their applications to mathematical finance can be found in Benth and Pircalabu (2018), Sabino (2020a) and Sabino and Cufaro Petroni (2021b) in the context of energy markets, in Bianchi and Fabozzi (Aug 2015) for the modeling of credit risk and in Barndorff-Nielsen and Shephard (2001) for stochastic volatility modeling.
The distributional properties of a non-Gaussian process of OU-type are closely related to the notion of self-decomposability (sd), because as noted in Barndorff-Nielsen and Shephard (2001) and in Taufer and Leonenko (2009), the stationary law of such a process must be sd. We recall here that a law with characteristic function (chf ) η(u) is said to be sd (see Sato 1999;Cufaro Petroni 2008) when for every 0 < a < 1 we can find another law with chf χ a (u) such that η(u) = η(au)χ a (u). (1) Of course a random variable (rv) X with chf η(u) is also said to be sd when its law is sd, and looking at the definitions this means that for every 0 < a < 1 we can always find two independent rv's-a Y with the same law of X , and a Z a with chf χ a (u)-such that in distribution Hereafter the rv Z a will be called the a-remainder of X and in general has an infinitely divisible distribution (id) (see Sato 1999). We will show in the following (see also Barndorff-Nielsen 1998;Cufaro Petroni 2021a, Sabino 2020a) that the transition law between the times t and t + t of a Lévy-driven OU process X (t) essentially coincides indeed with the law of the a-remainder of its sd-stationary distribution, provided that a = e −b t where b is the OU mean-reversion rate. It is therefore natural to investigate the properties of the a-remainder of a certain sd law even irrespective of its possible relation to the theory of the OU processes.
In this study we focus our attention on the class of the tempered stable (TS) distributions (see for instance Rosinski 2007 andGrabchak 2016) with finite variation, and we provide a general framework to derive their transition laws from their associated a-remainders. There are in fact two standard ways to associate a TS distribution to an OU process X (·): if its stationary law is a TS distribution we will say that X (·) is a TS-OU process; if on the other hand, X (·) is driven by a TS background noise we will say that X (·) is a OU-TS process.
Some of the results discussed in the following sections are not new: Masuda (2011, 2012) and Zhang (2011), for instance, have considered OU processes whose stationary law is an exponentially-modulated TS distribution, hereafter called classical TS (CTS), whereas (Bianchi et al. 2017) have taken into account the rapidly decreasing TS laws (RDTS), and (Grabchak 2020) finally merged all these laws in the larger class of the general TS distributions. Recently (Qu et al. 2021) have also studied both CTS-OU and OU-CTS processes.
We approach the theory of OU processes differently namely, from the point of view of the a-remainder's of the stationary law. To this end, the first contribution of this inquiry is the characterization of the Lévy triplet of the a-remainder of a sd law: this would of course constitute a crucial building block in the construction of a Lévy-driven OU processes. On the other hand, we remark that-by exploiting properties valid for every id distribution-an explicit knowledge of the Lévy triplet of the a-remainder makes very easy and straightforward the calculation of the cumulants of the transition law of a OU process. This turns out in particular, to be a remarkable asset in testing the efficiency of the simulation algorithms and can be adopted for the parameters estimation. It is also worthwhile remarking that the laws of the a-remainders are in fact id, and therefore they lend the possibility of producing an entire new class of associated Lévy processes as shown for instance in Gardini et al. (2021Gardini et al. ( , 2022a.
In this perspective, the second contribution of the present paper consists in harmonizing the know results relative to TS processes of OU type in the scheme of the a-remainders, and in showing the further advantages of this approach: some of the theorems in the aforementioned literature become indeed special cases of this proposed comprehensive framework. Also, compared to the existing literature we can test the simulation algorithms more thoroughly taking higher order cumulants instead of just the mean and the variance.
Finally, as done in Qu et al. (2021), we focus our attention on the case of the CTS related OU processes with their transition laws, both for the OU-CTS and the CTS-OU cases, but within the perspective of the a-remainders. The third contribution is the derivation of a few new algorithms intended to simulate the skeleton of such processes which are tens of times faster than those available in the existing literature. We find in particular that for the simulation of the CTS-OU processes our procedure is computationally more efficient and by far faster than that of Zhang (2011) because it does not rely on an acceptance rejection method (other than that required to draw from a CTS law), but rather on the inverse method (see Devroye 1986). For OU-CTS processes instead, we designed a new procedure based on an acceptance rejection method that, at variance with that of Qu et al. (2021), has the advantage of having an expected number of iterations before acceptance that can always be kept arbitrarily close to 1. More important, we focus on computational efficiency and propose a flexible approximation which can even skip such an acceptance rejection step which leads to a remarkably faster solution compared to that of Qu et al. (2021). Finally, our numerical experiments also suggest a simple strategy for the parameter estimation.
The paper is organized as follows: the Sect. 2 introduces the notations and the preliminary notions. In particular it presents the basic properties of non-Gaussian OU processes and their relation with the sd laws. In the Sect. 3 we then derive the Lévy triplet of the a-remainder of an arbitrary sd law, and in particular that of the a-remainder of a general TS distribution with finite variation. These results are instrumental to explicitly write down the transition law of CTS-OU process. In the subsequent Sect. 5 we focus on the OU-CTS processes, and in Sect. 6 we present the algorithms for the simulation of the skeleton of both the CTS-OU and the OU-CTS processes pointing out their differences and their advantages with respect to the solutions already existing in the literature. The Sect. 7 illustrates the effectiveness of our simulation schemes by comparing the true values and the Monte Carlo estimated values of the first four cumulants. We also consider some approximation schemes to further check the performance of our procedures and we propose a simple approach to the parameters calibration. Finally the Sect. 8 concludes the paper with an overview of future inquiries and of possible further applications.

Notations and preliminary remarks
Before proceeding, we introduce some notation relative to the various distributions that we consider in the paper. We write G(α, β) to denote the gamma distribution with shape parameter α > 0 and rate parameter β > 0. We write U([0, 1]) to denote a uniform distribution in [0, 1] and P(λ) to denote the Poisson distribution with mean λ > 0. Finally, we use the shortcuts pdf, lch and iid for probability density function, logarithmic characteristic and independent and identically distributed, respectively.
Take a-possibly non-Gaussian-one-dimensional Lévy process L(·), and the OU process X (·) which is solution of the stochastic differential equation (SDE) to wit Hereafter L(·) will be called background driving Lévy process (BDLP) and L(t) will also represent its stationary increment of width t that completely define the process, but for an arbitrary initial condition. It is known that the chf ϕ L (u, t) of these increments, and their lch ψ L (u, t) = ln ϕ L (u, t) are retrievable from a given id law with chf ϕ L (u) = e ψ L (u) according to Based on the Lévy-Khinching representation theorem (see Sato 1999 pp. 37-38) ψ L (u) can be characterized in terms of the generating triplet (γ g , σ, ν) g : g is called the truncating function. In this paper, unless explicitly mentioned, we consider g(x) = 1 |x|≤1 and write the relative generating triplet (γ, σ, ν). Following a Barndorff-Nielsen and Shephard (2001) convention, if D is the stationary distribution we will say that X (·) is a D-OU process; when on the other hand the rv L(1) is distributed according to the id law D we will say that X (·) is an OU-D process.
A well-known result (see for instance Cont andTankov 2004 or Sato 1999) states that a distribution D can be the stationary law of some OU-D process if and only if D is selfdecomposable (sd, see more below). In addition, just by taking an arbitrary degenerate initial condition X 0 = x 0 , P-a.s.-and if we can manage to retrieve the distribution of its second, integral term Z (t)-from the pathwise solution (4) it is also apparently possible to deduce the transition law of the Markov process X (·), and therefore all its distributional details in an explicit form. It is appropriate to point out moreover that in the Eq. (4) we provided the solution in terms of the original BDLP L(t) rather than of the dimensionless time BDLP L(b t) as done in Barndorff-Nielsen and Shephard (2001): therefore a few results of ours will turn out to be explicitly dependent on the parameter b. The differences between these two equivalent representations are also discussed in Barndorff-Nielsen (1998), Barndorff-Nielsen and Shephard (2003) or Schoutens (2003) page 48.
Going back now to the SDE (3) it is possible to see (see also Barndorff-Nielsen et al. 1998) that the solution process (4) is stationary if and only if its chf ϕ X (u, t) is constant in time and coincides with the chf ϕ X (u) of the (sd) invariant initial distribution that turns out to be decomposable according to where now, at every given t, ϕ Z (u, t) = e ψ Z (u,t) denotes the chf of the rv Z (t) in (4). This last statement apparently means that the law of Z (t) in the solution (4) coincides with that of the a-remainder of the sd, stationary law ϕ X provided that a = e −b t , and that moreover we have It turns out therefore that studying the transition law of an OU process essentially amounts to find first its stationary law, and then the law of its a-remainder (5): it is easy indeed to see from (4) that the chf of the time homogeneous transition law with a degenerate initial condition X (0) = x 0 , P-a.s. is As a consequence one can focus on the properties of the a-remainders of these sd distributions in order to deduce the transition pdf of the associated OU processes. A number of additional relations between the distribution of the stationary process and that of the rv L(1) are known: it is possible to show for instance, that between the lch's ψ X (u) = ln ϕ X (u) of the stationary distribution, and ψ L (u) = ln ϕ L (u) of L(1) the following relation holds (see Taufer andLeonenko 2009 andSchoutens 2003) On the other hand, assuming for simplicity that the Lévy measure of L(1) and that of the stationary process admit densities-respectively denoted as ν L (x) and ν X (x)and supposing that ν X (x) is differentiable, it also results (see Sato 1999;Cont and Tankov 2004) Taking advantage finally of (6), (7) and (8) it is easy to see that the transition lch of the OU process can also be written in terms of the corresponding ψ L (u) of L(1) in the form As a consequence, as already observed in Sabino and Cufaro Petroni (2021a) and Sabino (2020a), we can also calculate the cumulants c X ,k (x 0 , t), k = 1, 2, . . . of X (t) for X 0 = x 0 from the cumulants c L,k of L(1) according to On the other hand, according to (6), the said cumulants c X ,k (x 0 , t), k = 1, 2, . . . for X 0 = x 0 can also be derived from those of the stationary law here denoted c X ,k These quantities can be used both as benchmarks to test the performance of the simulation algorithms, and to carry out an estimation procedure based on the generalized method of moments.
In the following sections we focus our attention on the tempered stable processes (TS; see Rosinski 2007) with finite variation (for details see Cont and Tankov 2004), and we will analyze both the TS-OU and the OU-TS processes. To this end, we recall that the finite variation TS processes have Lévy densities of the form where the tempering term q(x) with q(0) = 1 is monotonically decreasing and q(+∞) = 0 for x > 0, and monotonically increasing and q(−∞) = 0 for x < 0. We do not adopt however the full characterization of Rosinski (2007) because we will focus on one-dimensional laws only, and mainly on the exponentially modulated TS, also known as (bilateral) classical tempered stable laws (CTS), where q(x) = e −β 1 x for x ≥ 0, while q(x) = e β 2 x for x < 0 with β 1 , β 2 > 0.

The TS distributions and their a-remainders
In the forthcoming sections we will discuss both the TS-OU and the OU-TS processes looking in particular to the properties of the a-remainder of their stationary laws, and we will focus our attention chiefly on the CTS subfamily. It is worthwhile noticing first that the study of these processes has extensively been carried on in the literature and that several types of TS laws have been investigated. For instance, Masuda 2011, 2012) and (Zhang 2011) have considered OU processes whose stationary marginal law is a CTS distribution, whereas (Bianchi et al. 2017) assume a rapidly decreasing TS law (RDTS) and (Grabchak 2020) finally harmonizes all these types of laws considering the larger class of general TS distributions. Albeit many results can consequently be found in the literature cited so far, we will nevertheless elaborate a little on this topic also to show how the proofs of the propositions can be carried out in a simple way by taking advantage of the properties of the a-remainders and of their Lévy triplets. For this purpose let us remember in particular that, as it is well-known, the sd laws constitute a subclass of the class of the id distributions having an absolutely-continuous Lévy measure with density where k(x) is increasing in (−∞, 0) and decreasing in (0, +∞) (see Cont and Tankov 2004, Proposition 15.3). Remark then that every TS distribution satisfying (16) also is sd. The law of the a-remainder of a sd law is id too (see Sato 1999) and the following proposition characterizes it in terms of its Lévy triplet. (γ, σ, ν), then for every 0 < a < 1 the law of its a-remainder has Lévy triplet (γ a , σ a , ν a ):

Proposition 1 Consider a sd law with Lévy triplet
Proof The Lévy-Khintchin representation of the lch of our sd law can be given in two equivalent ways taking 1 |x|≤1 and 1 |x|≤ 1 a as truncating functions, respectively Therefore, using both the representations (20) with the change of variable x = ax in the second integral, the lch of the a-remainder Due to the properties of k(x) it turns out that k(x) − k( x / a ) > 0 for every x and for every 0 < a < 1; and, as it also happens that with ν a (x) defined in (19) we have it is easy to see that ν a (x) qualifies as a Lévy density, then (γ a , σ a , ν a ) represents the Lévy triplet of the law of the a-remainder of the sd law.
In the context of the OU processes where the law of Z (t) in (4) at the time t is the aremainder of the stationary law for a = e −bt , the Proposition 1 along with the Eq. (9) enables us to connect the Lévy density ν Z (x, t) of the id rv Z (t) at the time t to ν X (x) and ν L (x), the Lévy densities respectively of the sd stationary law and of the BDLP L(1) at time t = 1: According to the previous representations we can therefore adopt one of two possible strategies to study the properties of the transition law of a Lévy-driven OU process: the first based on the stationary law and more suitable for a D-OU process; the second using the distribution of the BDLP and more suitable for an OU-D process.
A Lévy process is said to be of finite variation when its trajectories are of finite variation with probability 1, and it is possible to prove (Cont and Tankov 2004) that this happens if and only if its characteristic triplet (γ, σ, ν) satisfies the conditions Consequently, the lch of a process L(·) with finite variation and with triplet (γ , 0, ν) can represented as and it is called drift. Note that the Lévy triplet is not given by (d, 0, ν) but by (γ , 0, ν). An important subclass of such processes is that of subordinators that are Lévy processes with almost surely non-decreasing sample paths: in this event their Lévy triplet must satisfy the conditions As a matter of fact it would also be easy to see from the Lévy-Khintchin characterization theorem that any process of finite variation can be written as the difference of two independent subordinators-for instance the Variance Gamma (VG) processes can be represented as the difference of two Gamma processes-and therefore, without loss of generality, we can focus our attention on subordinators only. Without presuming now that the actual processes involved are of finite variation, subordinators or even Lévy processes, we will say hereafter for short that an id law with Lévy triplet (γ, σ, ν) is of finite variation when it satisfies the conditions (22), and is a subordinator when it satisfies the conditions (23). Moreover, for a TS subordinator we consider tempering functions which can be written as where p > 0 and Q(·) is a probability measure. Of course, taking p = 1 we get the set of completely monotone functions as in Rosinski (2007)). Finally, in the rest of the paper we will consider only subordinators with zero drift and we will denote CT S(α, β, c) a CTS law having a Lévy density (16) with q(x) = 0, x < 0, q(x) = e −βx , x ≥ 0, therefore with p = 1 amd Q(ds) = δ β (ds), where δ represents the Dirac delta.

Proposition 2 Consider a TS law with
Lévy measure ν such that ν((−∞, 0]) = 0 and Then, the Lévy density of its a-remainder is Proof As already remarked every TS law is sd and, if confined to x > 0, it is easy to see that (23) is satisfied so that the law is a subordinator and hence also of finite variation. As a consequence there is an a-remainder and from Proposition 1 we have is just a rescaled form of the original TS Lévy measure, we find that ν 2 (x) > 0 is also normalizable. From (24) and from Gradshteyn and Ryzhik (2007), 3.434.1 we have indeed that and therefore the integral in (28) turns out to be finite.
Under the conditions of the Proposition 2 we have found that ν 2 (x) is an integrable non-negative function therefore, taking a = ∞ 0 ν 2 (x)dx g a (x) = ν 2 (x)/ a can be interpreted a full fledged pdf. As a consequence ν 2 (x) can be considered as the Lévy density of a compound Poisson law of parameter a and jump length distributed according to the pdf g a (x). Since on the other hand ν 1 (x) is the Lévy density of the original TS law but rescaling the parameter c to c(1 − a α ), according to the previous proposition we can claim that the a-remainder Z a of a TS law with Lévy density (25) is-in distribution-the sum Y a +C a of two independent rv's: a TS Y a of the type (25) with parameters α, c(1 − a α ) and measure Q(·), and a compound Poisson rv where N s is a Poisson rv of parameter a and J a,k , k = 1, 2, . . . are a sequence of iid rv's with pdf g a (x) As can be seen from (7) in the Sect. 2, the transition law of an OU process (4) directly follows from the a-remainder of its sd stationary distribution when we take a = e −b t . Adding in the results of the Proposition 2 we can therefore fully display the Lévy measure of the said transition laws when the stationary distribution is a TS subordinator with Lévy density (25). This class of laws is not without merits in itself and is a relevant one-dimensional subfamily of the general tempered stable laws discussed in Grabchak (2020). It includes on the other hand the CTS subordinators with 0 ≤ α < 1, while in fact the Theorem 1 in Zhang and Zhang (2008) for Inverse Gaussian-OU processes (IG-OU) and the Theorem 1 in Zhang (2011) for TS-OU are special cases of the Proposition 2. The proposition also covers the p TS-OU introduced in Grabchak (2016), and for p = 2 the case of Rapidly Decreasing TS-OU (RDTS-OU) discussed in Kim et al. (September 2010) and Bianchi et al. (2017), the Modified TS-OU (MTS-OU) studied in Kim et al. (2009), and the Bessel TS-OU (BTS-OU) discussed in Chung (2016). It is expedient to notice at this point that, although for a fixed time t the law of Z (t) defined in (4) is id and coincides with that of the aremainder of the stationary law taking a = e −b t , this process is not Lévy because a changes in time.
The particular case of a TS-OU process with α = 0 (namely a G-OU process if the stationary law is a CTS) is also noteworthy: in this event indeed the BDLP turns out to be just a compound Poisson process. In fact it is easy to see that, for x > 0 and α = 0, from (9), (21) and (25) Therefore the BDLP will be a compound Poisson process where now N (t) is a Poisson process with intensity λ = c b, and J k are iid jumps with pdf (remember that q(x) is supposed to be non increasing). Then the pathwise solution (4) of the OU Eq.
(3) becomes now where τ k represent the jumping times of the Poisson process N (t). Of course this representation is valid for any BDLP compound Poisson and not only for that in (29). This type of Z (t) has also interesting financial applications beyond the context of OU processes: it can indeed describe random cash-flows occurring at random timeto-maturities with a rate of return equal to b. Remark moreover that-as observed by Lawrance (1980) in the context of Poisson point processes-for every t > 0 we have

Finite variation CTS-OU processes
In this section we focus on CTS-OU processes with finite variation, namely the subclass of OU processes with stationary distribution CT S α, β, c with 0 ≤ α < 1. Such a subclass is especially manageable and in this particular case the Proposition 2 entails indeed the following result.

Proposition 3
The CTS-OU process X (t) with initial condition X (0) = X 0 , P-a.s., and with stationary distribution CT S α, β, c whose Lévy density (16) with 0 ≤ α < 1, x > 0 has the tempering function q(x) = e −β x , β > 0, can be represented as andJ i are iid jumps, independent from N a and distributed according to the pdf Proof The law of Z (t) in the pathwise solution (4) coincides with that of the aremainder Z a of the stationary law CT S α, β, c whose Lévy density is From Proposition 2 taking p = 1 and Q(ds) = δ β (ds) we then have Then, ν 2 (x) is associated to the law of a compound Poisson rv with parameter a and jumps distributed according to f J (x) = ν 2 (x)/ a . On the other hand, since we can also write and this concludes the proof.

Remark 1
The law with pdf f J (x) coincides with the DTS distribution of Zhang (2011): from Eq. (36) we see that it is a mixture of a gamma law G(1 − α, βV ) with a random V distributed according to the pdf It is easy to verify moreover that where U ∼ U(0, 1) is a uniform rv, and therefore its simulation can be based on standard routines. In particular, when α = 1/2, X (·) turns out to be an IG-OU process and the simulation of its skeleton no longer requires now the acceptance-rejection methods adopted in Zhang and Zhang (2008) and in Qu et al. (2021), but can be based on the method illustrated in Michael et al. (1976) (see also Devroye (1986) page 148).

Finite variation OU-CTS processes
In this section we will consider Lévy-driven OU processes whose BDLP is a CTS process with q(x) = e −βx . In this specific case, the stationary law is not known in an explicit form (see for instance Table 2 in Barndorff-Nielsen and Shephard (2003)), but according to our discussion in the Sect. 2 the transition law of the solution (4) of the Eq. (3) can nevertheless be retrieved through the formula (7) if the law of Z (t) is known (remember that Z (t) remains the same for every initial condition). On the other hand we have shown that the formula (21) enables us to deduce the Lévy measure density of Z (t) at a given t directly from the Lévy measure density of the BDLP. We will show thus in the present section that the transition law of our OU-CTS process is the convolution of a CTS law (with parameters different from that of the BDLP) and a compound Poisson law.
(3) with L(1) ∼ CT S α, β, c and with X (0) = X 0 , P-a.s. is in distribution the sum of three independent rv's where X 1 is distributed according to the law CT S α, β a , c 1−a α α b , while and J k are iid rv's with pdf Remark 2 The pdf f J (x) (41) can be seen as a mixture of the gamma laws G(1−α, βV ) with a random rate parameter V distributed according to the pdf which is correctly normalized. Proposition 4 covers the case of a OU-gamma process illustrated in Qu et al. (2019) when α tends to zero: we have indeed lim Proof Based on Eq. (21) and with the change of variable y = wx, the Lévy density of the term Z (t) in the pathwise solution (4) of an OU-CTS process is (remember that the coefficient a = e −bt is time dependent) The first term apparently is the Lévy density of a CTS law CT S α, β a , c 1−a α α b because it is easy to see that On the other hand ν 2 (x) > 0 for every x > 0 because e −βwx − e − β a x > 0 when 0 < w < 1 / a , and moreover with v = aw we find (see 3.434.1 in Gradshteyn and Ryzhik (2007) where apparently 0 < a < +∞. As a consequence is a valid pdf and then, ν 2 (x) represents the Lévy density of a compound Poisson law with parameter a and jumps distributed according to the pdf f J (x). It would be possible to show now (see 3.381.3 in Gradshteyn and Ryzhik (2007)) that where (γ , z) is the incomplete gamma function. For later computational convenience however, we prefer to give an alternative representation of this jumps distribution. Since and with an exchange in the order of the integrations, the pdf (43) becomes (41) and is a mixture of gamma laws G(1 − α, βV ) with a random rate parameter distributed according to the pdf (42), as stated in the Remark 2.

Simulation algorithms
The simulation of CTS laws have been widely discussed in several studies (see for instance Masuda (2011, 2012), Zhang (2011) and Grabchak (2019Grabchak ( , 2020 and the references therein) and several software packages are available for such purpose. Therefore, in this section we will only illustrate how to exactly simulate the CTS-OU and the OU-CTS processes that, although similar in names, are two rather different objects as explained in the previous sections. As far as the CTS-OU processes are concerned, our contribution is the enhancement of the simulation performance by taking advantage of the representation (36) for f J (x) that, in contrast to Zhang (2011), does not require any acceptance-rejection procedure. On the other hand, with regard to OU-CTS processes, we design a new simulation procedure for the drawings from the mixture with pdf (42). At variance with the approach of Qu et al. (2021), this algorithm is based on an acceptance-rejection method whose expected number of iterations before acceptance however can be made arbitrarily close to one and is therefore more efficient. More important, the dominating pdf is monotone and convex and can be made arbitrary close to the pdf of the mixture, and thus, the acceptance rejection step can be skipped leading to a very fast solution. In our numerical experiments we consider a time grid t 0 , t 1 , . . . , t M , t m = t m − t m−1 , m = 1, . . . , M with M steps.

CTS-OU processes
The simulation procedure for the generation of the skeleton of CTS-OU process is based on the Proposition 3 and is summarized in the Algorithm 1. We remark that when Generate an independent Poisson rv with a in (35) 6 u i ← U i ∼ U (0, 1), i = 1, . . . , n Generate n iid uniform rv's . . , n Generate n independent gamma rv's all with the same scale 1 − α and random rates 10 x 2 ← n i=1 j i 11 X (t m ) ← a X(t m−1 ) + x 1 + x 2 . 12 end for α = 0 a CTS-OU process is a compound Poisson process with a gamma stationary law whose efficient exact simulation can be found in Sabino and Cufaro Petroni (2021a).

OU-CTS processes
The simulation steps for the skeleton of a OU-CTS process are then summarized in the Algorithm 2. The sampling from a CTS law has been widely studied by several Generate an independent Poisson rv with a in (35) Generate n iid rv's with pdf given by Equation (42) Generate n independent generalized gamma rv's all with the same p, scale p − α and random rates 9 x 2 ← n i=1 j i 10 X (t m ) ← a X(t m−1 ) + x 1 + x 2 . 11 end for authors (see for instance Devroye (2009) and Hofert (2012)), and here the only nonstandard step is the fifth one in the Algorithm 2, namely that allowing the generation of the jumps of the compound Poisson process of Proposition 4. On the other hand, as mentioned in Remark 2, these jump sizes are iid rv's following a gamma law with shape 1 − α and a random rate. Therefore the unique remaining task is to sample from a law with the pdf f V (v) (42). Since however this pdf is not monotonic in [1, 1 / a ] for every value of its parameters, we first define the new rv It is straightforward to check then that f W (w) is monotonic and convex in [0, 1] and hence one can rely on the inversion-rejection algorithm illustrated in Devroye (1986) page 355. The solution that we propose here is very similar to that, and in effect consists in replacing the steps required for the sequential search of the inversion part with the method of partitioning the densities into intervals (see once again Devroye 1986 page 67). We remark indeed that f W (w), besides being monotonic and convex, also has the following upper bound namely it is dominated by a linear function where G(a, α) is the area under g(x). We could therefore devise a simple acceptance-rejection procedure where G(a, α) should be as close to 1 as possible because it roughly represents the number of iterations needed in the rejection algorithm. While however G(a, α) → 1 + when a → 1 − , unfortunately it is G(a, α) → +∞ for a → 0 + . Taking therefore a = e −b t , this latter limit means that the generation of a OU-CTS process with a large b t might either have a heavy computational cost, or potentially require a large number of simulations.
In principle we could consider only small time steps, but on the other hand the acceptance-rejection sampling can be easily improved, via the modified decomposition method elucidated in Devroye (1986) where we can also write Apparently theḡ (w), = 1, . . . , L turn out to be piecewise linear pdf 's, while the p constitute a discrete, normalized distribution. Increasing the number L of the intervals, with given 0 < a < 1 and 0 < α < 1, G L (a, α) can be made arbitrary close to 1 because it measures the trapezoidal approximation of 1 0 f W (w)dw = 1. On the other hand the random drawing from the laws with pdf 'sḡ (w), = 1, . . . , L is very simple and can be implemented via the standard routines and is extremely fast. Based on this last observation and on the trapeziodal rule, one can select a relatively high number L, approximate f W (w) with g L (w) and avoid running the acceptance-rejection procedure. The numerical experiments illustrated in Sect. 7 show that such an approximation returns cumulants which are consistent with their theoretical value. Moreover it is easy to implement and brings a remarkable reduction of the computational time, especially for programming languages like Python, R and MATLAB because for loops can be avoided and the code can be 'vectorized'.
Denoting now with S a rv with distribution P {S = } = p , = 1, . . . , L, and with Y a rv with pdfḡ (w), the Algorithm 3 summarizes the instructions needed to implement the sixth step in the Algorithm 2. We remark finally that an alternative procedure, leading to similar results, might have been some shrewd decomposition of f W (w) rather than of its dominating curve. However Devroye (1986) at the page 70 nicely spell out the reasons why the procedure here adopted is in principle preferable.

Algorithm 3
Remark 3 It is worthwhile mentioning that an alternative procedure relying on a different acceptance-rejection strategy has been proposed in Qu et al. (2021). In contrast to this last approach, however, in our algorithm G L (a, α) can be made arbitrary close to 1 irrespective of the value of a (and of the size of the time-step), and therefore our approach turns out to be computationally more efficient.
On the other hand we can also take advantage of the interplay between the OU-CTS and the CTS-OU processes to gain an insight into the possible benefits of the different simulation strategies: from (10) we find indeed that the Lévy density of the BDLP L(t) for a CTS-OU process is where the first term apparently provides the BDLP L 1 (t) of an OU-CTS, whereas the second term corresponds to a compound Poisson L 2 (t) (see Cont and Cont and Tankov 2004 page 132). Therefore the path-wise solution (4) of our CTS-OU process is now where Z 1 is a OU-CTS (with Z 1 (0) = 0) and, as shown in (31), Z 2 is a compound Poisson that can easily be simulated based on (32). On the other hand, according to the Proposition 4, the OU-CTS process Z 1 (t) is in its turn the sum of a (time-dependent) CTS rv X 1 and of a (time-dependent) compound Poisson rv X 2 ; so that ultimately a CTS-OU process with a degenerate initial condition X 0 = x 0 -beyond being of the form (33) presented in the Proposition 3-can now be seen also as the sum of four random terms: one distributed according to a CTS law, two compound Poisson rv's and a degenerate summand. Of course the two representations coincide in distribution and, as a matter of fact, the four-terms representation reproduces again that of Qu et al. (2021); but the possible alternative simulation algorithms stemming from the four term representation, although perfectly correct, would require now the generation of three rv's and the use of acceptance-rejection methods in addition to that needed for the sampling of a CTS distributed rv, and therefore they would turn out to be rather less efficient than the Algorithm 1.
Finally, based on Remark 2 and on the results of Qu et al. (2019), we notice that for α = 0 the simulation of V is here much easier because no acceptance-rejection method is required.

Numerical experiments
In this section, we will assess the performance and the effectiveness of our algorithms through extensive numerical experiments. All the simulation experiments in the present paper have been conducted using Python with a 64-bit Intel Core i5-6300U CPU, 8GB. The performance of the algorithms is ranked in terms of the percentage error relative to the first four cumulants denoted err % and defined as err % = true value − estimated value true value Finally, we focus on the computational performance of our solutions in comparison to that of the algorithms in the existing literature.

CTS-OU processes
Since the Lévy density on [0, +∞) of the stationary law CT S (α, β, c) of a CTS-OU process with finite variation and parameters α, β, c is its cumulants are (see for instance Cont and Tankov 2004, Proposition 3.13) and therefore from (14) and (15) we obtain the cumulants of X ( t) with the degenerate initial condition In our numerical experiments we consider a CTS-OU process with x 0 = 0 and parameters (b, c, β) = (10, 0.8, 1.4) whose trajectories with α ∈ {0.3, 0.5, 0.7, 0.9} are displayed in Fig. 1 where of course the case α = 0.5 is that of an IG-OU process. We remark that the sampling from an IG law can be performed via the many-to-one transformation method of Michael et al. (1976), and therefore no acceptance-rejection procedure is required in Algorithm 1 to generate the skeleton of an IG-OU process.
The Tables 1 and 2 compare then the true values of the first four cumulants c X ,k (0, t) with their corresponding estimates from 10 6 simulations respectively with t = 1/365 and t = 30/365. We can conclude therefrom that the proposed Algorithm 1 produces cumulants that are very close to their theoretical values. Table 3 focuses on the computational times required by each term in Eq. (33) and highlights the difference between our methodology and that proposed in Zhang (2011). The numerical results clearly show that our solution is at least ten times faster because, in contrast to Zhang (2011), it avoids the acceptance rejection step required for the simulation of the compound Poisson component X 2 .
For the sake of brevity, we do not report the additional results obtained with different parameter settings that anyhow bring us to the same findings. Overall, from the numerical results reported in this section, it is evident that the Algorithm 1 proposed above can achieve a very high level of accuracy as well as a conspicuous efficiency.   Comparing the first four true cumulants with their corresponding MC-estimated values (multiplied by 10) obtained with R = 10 6 replications and t = 30/365 Computational times in seconds of X 1 and X 2 in Eq. (33). Algorithm 1, steps 5 to 10, vs Algorithm in Zhang (2011) The fact that we can easily compute the cumulants of an OU process substantiates the advantages of focusing our treatment on the law of the a-remainder of its stationary distribution. In addition to both the simple derivation of the transition pdf and the detailed testing of its statistical properties, we could indeed also conceive a parameter estimation procedure based on the generalized method of moments (GMM). We remark finally that the law of an a-remainder always is id, and therefore a simple modification of the simulation procedure presented in the Algorithm 1 could be adopted for the generation of a Lévy process whose law at time t = 1 is that of the a-remainder of a CTS distribution.

OU-CTS processes
Here too we will benchmark the results of the numerical experiments against the true values of the first four cumulants of OU-CTS process at time t with X (0) = 0. From the formula (46) for the cumulants of a CT S (α, β, c) distribution, and from (12) and (13) we first recover indeed the cumulants of X ( t) with the degenerate initial For our simulations we consider then the same parameter settings of the previous section-x 0 = 0, (b, c, β) = (10, 0.8, 1.4)-adapted to an OU-CTS process, and with α ∈ {0.3, 0.5, 0.7, 0.9} we get the sample trajectories displayed in the Fig. 2, where of course the case α = 0.5 is that of an OU-IG process.
In order to generate a OU-CTS process, we adopt four strategies. This first choice is Algorithm 2 where X 2 is generated with plain Algorithm 3 taking L = 100 equally spaced subintervals of [0, 1]. As observed in Sect. 6, the pdf f W (w) is increasing and convex in [0, 1] which can be partitioned into a sufficiently large number of sub-intervals such that f W (w) can be approximated by a simple step-wise linear pdf. Based on this observation, the second one is a very good approximation still based on Algorithm 2 where X 2 is simulated taking L = 2000 equally spaced subintervals of [0, 1] and approximating the law with pdf f W (w) with that with pdf g L (w) (properly normalized); the acceptance rejection steps are then skipped. Other than reducing the number of computations, this strategy has the advantage to avoid for loops in programming languages like Python, MATLAB and R whose computational performance otherwise normally suffers some limitations.
In addition to these two configurations for Algorithm 2, we also consider here two other approximate procedures: the first boils down to simply neglect X 2 in the Proposition 4; the second-in the same vein of  dealing with the normal inverse Gaussian-driven OU processes-takes advantage of the approximation of the law of Z (t) in (4) with that of e −k t L(t) where L(t) ∼ CT S α, β / a , c t .
Tables 4 and 5 compare then the true values of the first four cumulants c X ,k (0, t) with their corresponding estimates for 10 6 simulations with t = 1/365 and t = The Tables 4 and 5 clearly show that the cumulants returned by Algorithm 2 are equally close to the theoretical values irrespective to the choice of t and to the configuration. Hence, the trapezoidal approximation of f W (w) returns an efficient approximation as already noted in Devroye (1986) pages 67 and 68. Moreover, we can now conclude that our Algorithm 2 has the lowest percent errors, but nevertheless, in some practical situations, the errors of two other approximations could be deemed acceptable taking also into account that their computational cost is lower. In particular, the second alternative procedure outperforms the third one and its percent errors are not much higher than those of the exact method. When however the time step is larger, or equivalently when a = e −b t is close to 0, the three procedures give radically different outcomes and, as it is shown in the Table 5, the two approximate methods return completely biased results. Conversely, the exact method and the related approximate method continue to be reliable and its percent errors remain small even for the higher cumulants.
The previous state of affairs for an OU-CTS is due to the fact that X 2 in Proposition 4 produces only a second order effect when t → 0 + : using indeed a Taylor expansion  we find and therefore the compound Poisson X 2 has a relevant impact only when t is not too small. Notice instead that this is not how a CTS-OU process behaves because from the Proposition 3 we see that for t → 0 + so that X 2 results in a first order effect and cannot be neglected even for small t. As mentioned above, to tackle the simulation of the random rate, Qu et al. (2021) have proposed an alternative solution that is based on an acceptance rejection method again, and whose expected number of iterations before acceptance tends to 1 for small time steps ( t → 0 + ; a → 1 − ); but unfortunately this value somehow deteriorates and tends to 2 (50% of acceptance) for large time steps ( t → +∞; a → 0 + ). From our previous findings we know instead that neglecting X 2 is a fair approximation for finer time grids so that the impact on the computational cost of the acceptance rejection is rather restricted. On the other hand, no matter how large the time step t is, with our algorithm the expected number of iterations before acceptance can be kept as close to 1 as possible because it depends on the accuracy of a trapezioidal approximation. Therefore, recalling also that the computational cost to generate a simple discrete rv is very low, our approach turns out to be computationally more efficient. However, from the computational point of view, one may argue that it may be nevertheless faster to generate a trajectory with 30 points with t = 1/365 instead of one single point with t = 30/365. To this end, Table 6 summarizes the computational times in seconds of X 1 and X 2 in Eq. (39) and for the latter component, it highlights the difference among the plain Algorithm 2 taking L = 100, Algorithm 2 with the approximation of f W (w) with g L (w) taking L = 2000, and the methodology of Qu et al. (2021). The computational times of Algorithm 2 with the aforementioned approximation are tens of times smaller than those returned by the other solutions. For instance, the computational times of a OU-IG process is now a matter of seconds and not of tens of seconds anymore with R = 10 6 replications. We remark that the computational time of the plain Algorithm 2 is also faster than that of Qu et al. (2021) because the former one has a lower number of operations before acceptance. We also stress the fact that we paid attention to an efficient implementation of the methodology of Qu et al. (2021) indeed, the computational times in Table 6 are smaller than those reported in the original paper although we used a computer with less computational power. We also remark that the computational time of X 2 with L = 2000 taking f W (w) ≈ g L (w) and t = 30/365 is less than 30 times that of just X 1 with t = 1/365 in contrast to plain Algorithm 2 and the methodology of Qu et al. (2021). We conclude then that our methodology provides a remarkable computational advantage and is the preferable one.
These observations could also lead to a convenient strategy combining parameters estimation and exact simulation of the OU-CTS processes. Assuming that the data could be made available with a fine enough time-granularity (e.g. daily t = 1/365), we could base the parameters estimation on the likelihood methods by approximating the exact transition pdf with that of a CTS law CT S α, β a , c(1−a α ) α b . However, to avoid being forced to always simulate the OU-CTS processes on a fine time-grid allowing the approximations (for instance if one needs to simulate it at a monthly granularity t = 30/365), the generation of the skeleton of such processes will be preferably based on the exact method of the Algorithm 2.

Conclusions
In this paper we studied the transition laws of the tempered stable related OU processes with finite variation from the standpoint of the a-remainders of sd distributions: in fact, Table 6 OU-CTS Computational times in seconds of X 1 and X 2 in Eq. (39). Algorithm 1, steps 5 to 10, vs Algorithm in Qu et al. (2021) the transition law of any OU process essentially coincides with the distribution of the a-remainder of its stationary sd distribution. To this purpose, we first derived the Lévy triplet of the a-remainder of a general sd law which is then instrumental to find the representation of the transition law of tempered stable related OU processes with finite variation. We thereafter focused our attention on the CTS-OU and the OU-CTS processes: respectively those whose stationary law is a CTS distribution, and those whose BDLP is a CTS process. As already done in Zhang andZhang (2008, 2009), Kawai and Masuda (2011) and Qu et al. (2021), we showed that their transition law coincides with the distribution of the sum of a CTS distributed rv (with scaled parameters), of a suitable compound Poisson rv and of a degenerate term: we accordingly also derived their path-generation algorithms.
As for the simulation of the skeleton of CTS-OU processes, we focused on the computational efficiency indeed, our proposed procedure amounts to a remarkable improvement with respect to the existing solutions presented in Zhang and Zhang (2008), Zhang (2011), Kawai and Masuda (2011): it does not rely on additional acceptance rejection methods other than that required to generate a CTS distributed rv. Our computational times are by far faster and are suitable for real-time applications. On the other hand, the simulation procedure for a OU-CTS process is based on an acceptance rejection approach more efficient than that described in Qu et al. (2021), because here the number of iterations before acceptance can be made arbitrarily close to 1 no matter how fine we choose the time grid of the skeleton. We also proposed an approximation that even avoids running the acceptance-rejection step bringing a remarkable reduction of the computational time especially if the solution is implemented in programming languages like Python, R and MATLAB. Our numerical experiments show that this last procedure returns correct cumulants and computational times which are tens of times smaller than the other solutions.
Although we considered in the present paper only the CTS distributions restricted on the positive real axis, the results can be easily extended to the bilateral case and the simulation of the relative processes would be simply obtained by running twice the proposed algorithms. A further object of our future inquiries will be instead the possible extension to the p -TS related OU processes (see Grabchak Grabchak (2021b)) combined with the application of the algorithms recently proposed in Grabchak (2021a) to draw samples from p -TS laws. We remark moreover that, due to the fact that the laws of the a-remainders are id, our approach is also suited to build and simulate new Lévy processes via the subordination of a Brownian motion with the Lévy process generated by the a-remainder of a gamma and IG law, respectively (as done for instance in Gardini et al. (2021Gardini et al. ( , 2022a).
All these algorithms would finally be especially useful for a simulation-based statistical inference, and for some financial applications like as the derivative pricing and the value-risk calculations. To this end, a possible future research line could be the study of the time reversal simulations in the spirit of some recent papers by Pellegrino and Sabino (2015) and Sabino (2020b) relatively to the time-changed OU processes introduced in Li and Linesky (2014).
Funding Open Access funding provided by University of Helsinki including Helsinki University Central Hospital.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.