Exact simulation of normal tempered stable processes of OU type with applications

We study the Ornstein-Uhlenbeck process having a symmetric normal tempered stable stationary law and represent its transition distribution in terms of the sum of independent laws. In addition, we write the background driving Lévy process as the sum of two independent Lévy components. Accordingly, we can design two alternate algorithms for the simulation of the skeleton of the Ornstein-Uhlenbeck process. The solution based on the transition law turns out to be faster since it is based on a lower number of computational steps, as confirmed by extensive numerical experiments. We also calculate the characteristic function of the transition density which is instrumental for the application of the FFT-based method of Carr and Madan (J Comput Finance 2:61–73, 1999) to the pricing of a strip of call options written on markets whose price evolution is modeled by such an Ornstein-Uhlenbeck dynamics. This setting is indeed common for spot prices in the energy field. Finally, we show how to extend the range of applications to future markets.


Introduction
Recent studies have shown a growing interest in the modeling with non-Gaussian Lévy processes of Ornstein-Uhlenbeck (OU) type typically for financial applications. Indeed, for instance, interest rates and in particular the spot prices of energy and commodity markets exhibit mean-reversion and spikes. Such features are hard to capture with plain Lévy processes and even with Gaussian OU processes. As observed in a series of papers by Barndorff-Nielsen (1998) and Barndorff-Nielsen andShephard (2001, 2003) the family of tempered stable (TS) processes represent a mathematically manageable class suitable for several financial applications.
The views, opinions, positions or strategies expressed in this article are those of the author and do not necessarily represent the views, opinions, positions or strategies of, and should not be attributed to E.ON SE.
Moreover plain TS processes, under certain parameter conditions, can be taken as subordinators to then time-change standard Brownian motions yielding another class of Lévy processes, named normal tempered stable processes (NTS), where the well-known variance gamma (VG) and normal inverse Gaussian (NIG) processes are particular cases. These TS subordinators are built by tempering the Lévy density of stable subordinators by a given tempering function. This is typically done with a classical exponentially decreasing function, but the approach can be conveniently extended taking completely monotone tempering functions á la Rosiński (2007) and can be generalized as done by Bianchi et al. (2017), Kim et al. (September 2010) and finally by Grabchak (2016Grabchak ( , 2019. On the other hand, when moving to the OU types the transition law is in general available in terms of characteristic functions or of Lévy measures which makes the parameter estimation and the path simulation non trivial. To this end, for TS-driven OU processes some new and efficient approaches have been proposed, for instance in Qu et al. (2019Qu et al. ( , 2021, Cufaro Petroni (2021b, 2022), Grabchak (2020) and Grabchak and Sabino (2022) covering the construction of Lévy-driven OU processes from either the stationary law or from the backward driving Lévy process (BDLP). When however, one considers NTS process of OU type, the state of art is not fully complete and the aim of this study is to fill this gap for symmetric NTS distributions.
There are in fact two standard ways to associate a NTS distribution to an OU process X (·): if its stationary law is a NTS distribution we will say that X (·) is a NTS-OU process; if on the other hand, X (·) is driven by a NTS background noise we will say that X (·) is a OU-NTS process. Indeed, a common convention states that if D is the stationary distribution we will say that an OU process X (·) is a D-OU process; when on the other hand, the background noise is distributed according to the infinitely divisible law D we will say that X (·) is an OU-D process. Even though these mathematical objects are labeled by similar names they are instead very different. Well-known examples of NTS distributions are the VG and NIG laws which then according to the aforementioned naming convention, lead to NIG-OU and VG-OU processes in the former case, and to OU-NIG and OU-VG processes in the latter one.
As far as we are aware of, in the current state of affairs, the transition density and efficient simulation methods are available for VG-OU processes (see Sabino and Cufaro Petroni (2021b)), for OU-VG processes (see Cummins et al. (2017) and Sabino (2020)) and for OU-NTS process with a symmetric distribution (see ) whereas, they are still outstanding for the combination NTS-OU which motivates this investigation. Actually, the special case of NIG-OU processes, will full asymmetric law, was already addressed in Barndorff-Nielsen (1998) without however providing an explicit information of the transition law but rather the BDLP only.
The first contribution of this work is the decomposition of the transition law of a symmetric NTS-OU process in terms of independent and manageable distributions. We also consider the more general tempering functions in the spirit of Grabchak (2016) that covers completely monotone and exponential functions as well. The second contribution is the calculation of the BDLP and its representation as the sum of independent processes: a NTS and a compound Poisson. Here instead, we take only a completely monotone tempering. These two results entail two alternate procedures for the simulation of the skeleton of a NTS-OU processes. More precisely, since one of the components of the BDLP is a NTS process one can rely on the findings of  in the context now of a OU-NTS process and apply them to the original problem. This solution however, is slower than the one which starts from the stationary law since it involves more computational steps as we also illustrate with extensive numerical experiments.
Having a fast and unbiased simulation algorithm is also instrumental for the parameter estimation which can accordingly be accomplished via Monte Carlo-based methods other than via moment (or cumulant) matching.
The simple decomposition in terms of the sum of independent random variables allows a simple computation of the characteristic function and the cumulative generating function of the transition law which is then an additional contribution. In the context of derivative pricing these quantities are instrumental for the derivation of the risk-neutral conditions and the application of FFT-based techniques (see Carr and Madan (1999)). We present an application to the pricing of a strip of daily European call options written on a market whose spot dynamics follows a NTS-OU process which is typical for energy and commodity markets. Of course, the range of applications can be enlarged to different contracts, for instance to Asian options using the Fourier techniques like those in Zhang and Oosterlee (2013), to storage facilities or swing options with the least-squares Monte Carlo method (see Boogert and de Jong (2008)), or even to the future markets in order to capture the volatility smile and the Samuelson effect (see for instance Benth et al. (2019), Latini et al. (2019) and Piccirilli et al. (2021)). Moreover, these FFT and MC approaches can be employed in model calibration from quoted options in which case the parameters are obtained by minimizing the distance between the quoted values and the model option prices computed either via FFT or via Monte Carlo simulations.
Finally, additional applications could be interest rate or credit modeling similar to those studied in Bianchi and Fabozzi (Aug 2015) and even physical modeling like collective motion in the charged particle accelerator beams as discussed in Cufaro Petroni (2007).
The paper is structured as follows. Sect. 2 sets up the scene and defines the notation. Sect. 3 illustrates the main results whereas, Sect. 4 presents the numerical experiments and the potential applications. Finally, Sect. 5 concludes the paper.

Notations and preliminary remarks
In this section we present the concept of self-decomposability, its relation to non-Gaussian OU processes and normal variance-mean mixtures. We also introduce the notation and the shortcuts that will be used throughout the paper.

Notation
We write N (μ, σ 2 ), σ > 0, to denote the Gaussian distribution with mean μ and variance σ 2 , U([0, 1]) to denote the uniform distribution in [0, 1] and P(λ) to denote the Poisson distribution with parameter λ > 0. We label GGa (d, p, β), where (x) represents the Euler gamma function. Of course, for p = 1 we retrieve the gamma distribution with scale d and rate β denoted (d, β) and if X ∼ GGa (d, p, β) then X p ∼ (d/ p, β p ). Moreover, W κ,μ (z) is the Whittaker function and K ν (x) is the modified Bessel function of the second kind (see Gradshteyn and Ryzhik (2007)). We use the shortcut sd for self-decomposable distributions and use the shortcut rvfor random variable and iid for independently and identically distributed. Finally, we use chf, lch, cgf, cdf and pdf as shortcuts for characteristic function, logarithmic characteristic function, cumulant generating function, cumulative distribution function and density function, respectively.

Self-decomposable laws and OU processes
A law with 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 Of course, a rvX with chf η(u) is also said to be sd when its law is sd, and this means that for every 0 < a < 1 we can always find two independent rv's, Y with the same law of X , and Z a with chf χ a (u) such that in distribution Hereafter the rvZ a will be called the a-remainderof X and is infinitely divisible (see Sato (1999)). Several works (see for instance Barndorff-Nielsen (1998), Cufaro Petroni (2021b, 2022) and Sabino (2020)) have pointed out that the transition law between the times t and t + t of a Lévy-driven OU process X (·) essentially coincides 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. To see that, take a non-Gaussian one-dimensional Lévy process L(·), the BDLP, and the OU process X (·) with meanreversion rate b > 0, which is solution of the stochastic differential equation (SDE) namely, Assuming the notation introduced in Section 1, 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 sd. In addition, 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 sd according to where now, at every given t, ϕ Z (u, t) = e ψ Z (u,t) and ψ Z (u, t) denote the chf and the lch of the rvZ (t) in (4), respectively. This last statement apparently means that the law of Z (t) in the solution (4) coincides with that of the a-remainderof the sd, stationary law provided that a = e −b t , as mentioned before.
A number of relations between the distribution of the stationary process and that of the rvL(1) are known, we present here only those ones relevant for this study. 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 results (see Sato (1999), Cont and Tankov (2004)) One can also explicitly calculate the cumulants c X ,k (x 0 , t), k = 1, 2, . . . of X (t) for X 0 = x 0 either from the cumulants c L,k of L(1) according to or from those of the stationary law here denoted c X ,k where δ i j denotes the Kronecker delta. 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.

Tempered stable subordinators
In the following sections we will focus our attention on the OU processes whose stationary law is a NTS law, that can be seen as the one of a Brownian motion (BM) at time t subordinated by a TS subordinator with no drift. To this end, we recall that a TS subordinator S(·) with no drift is a finite variation Lévy process with characteristic exponent and 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 q(x) = 0 for x < 0. We consider tempering functions which can be written as where Q(·) is a probability measure. This class is a relevant one-dimensional subfamily of the general TS laws discussed in Grabchak (2020). It covers for p = 1 the set of completely monotone functions as in Rosiński (2007), for p = 2 the case of Rapidly Decreasing TS-OU discussed in Kim et al. (September 2010) and Bianchi et al. (2017), the Modified TS-OU studied in Kim et al. (2009), and the Bessel TS-OU discussed in Chung (2016).

Normal variance-mean mixtures
In order to introduce NTS laws we recall that a rv Y is said to be distributed according to a normal variance-mean mixture with mixture pdf f V (x) if it has the form where X ∼ N (0, 1) and V , independent from X , has pdf f V (x) defined on the positive real axis R + , θ , μ and σ > 0 Assuming for simplicity μ = 0, the law of Y coincides with the distribution of the position of a subordinated BM where the density of the subordinator at a fixed time is given by f V (x). A NTS distribution is then a normal mean-variance mixture distribution with a TS mixture law, accordingly, a time-changed a BM with drift θ and volatility σ with a TS subordinator, is a NTS process. Of course taking α = 1/2 and Q(ds) = δ β (ds) we have the well-known NIG process. Finally, it is possible to see that NTS distribution admits Lévy density ν N (x) which, according to Theorem 4.3 in Cont and Tankov, can be represented as follows

Symmetric NTS-OU processes and their simulation
In the section we focus on the transition law of OU processes whose stationary law is a symmetric NTS. According the aforementioned notation, such a process is labeled SNTS-OU. The symmetric NTS distribution assumes μ = θ = 0 in (13) and hereafter will be dubbed as SN T S(α, c, Q, p, σ ) to reflect its parameters. According to (4) a SNTS-OU process can be represented as and L X (·) is its BDLP. We remark that the case α = 1/2 and Q(ds) = δ β (ds), β > 0 has been already studied in Barndorff-Nielsen (1998) without the symmetric constrain.
On the other hand, Barndorff-Nielsen (1998) does not focus on Z X (·) but rather on L X (·). In contrast, we focus on the former process and show that this different approach, irrespective to the extension to a general SNTS-OU process, not only explicitly gives information on the transition law of X (·) but also yields a faster simulation algorithm of its skeleton. Moreover, finding the chf and the mgf of X (·) is straightforward.

Derivation from the stationary law
is a compound Poisson rv where P ω is an independent Poisson rv with parameter and J k , k > 0 are iid rv's with density namely, a mixture of a generalized gamma law GGa( p − α, p, V S) with random rate being the product of two independent rv's: S with cdf F S (ds) = s α Q(ds)/q α and V having pdf Proof Taking a = e −bt , the law of of the rv Z X (t) coincides with that of the a-remainderof a SN T S(α, c, Q, p, σ ) distribution which is infinitely divisible therefore, we can rely on the Lévy-Khinchin machinery. According to Proposition 3.1 in Sabino and Cufaro Petroni (2022) the Lévy density ν Z (x, t) of Z X (t) can be written as By change of variable ω u = y in the second term of the integral and then changing back y = u the above integral can be re-written as On the other hand, according to Theorem 4.3 in Cont and Tankov, we note that ν 2 (x, t) is the Lévy density of a normal mean-variance mixture law whose mixing distribution has Lévy density where we have used the change of variables z = u p and 3.434.1 in Gradshteyn and Ryzhik (2007) hence, we can write that concludes the proof.
Remark 1 From the standard properties of the Gaussian law, it also holds where X ∼ N (0, 1) and it turns out to be a more convenient and a computationally efficient way to simulate the skeleton of X (·).
Proposition 3.1 also provides an algorithm for the simulation of the skeleton of an SNTS-OU process on a time-grid also detailed in Algorithm 1. Of course, taking Q(ds) = δ β (ds), steps 9 simplifies and the TS law turns out to be a classical TS distribution. The generation of the rv's in Algorithm 1 is standard except that for the classical TS distributed rv M 1 , for V and S that depends on the specific tempering function. To this end, the sampling from a classical TS law has been widely studied by several authors (see for instance Devroye (2009) and Hofert (2012)), of course, taking α = 0.5 the TS law coincides with an IG law hence, one can rely on the many-two-one transformation method of Michael et al. (1976). Finally, V can be generated via the inverse method after observing that

Derivation from the BDLP
In any case, from equations (8) we can also characterize the BDLP of an SNTS-OU process. Such a process has already been studied in Barndorff-Nielsen (1998) including the asymmetric case under the assumption that α = 1/2, Q(ds) = δ β (ds) and p = 1. It turns out that the BDLP is the sum of three independent processes, one of them simplifies out for the symmetric case and more important, can only be represented in terms of its che, therefore it is not straightforward to simulate it and its cdf can only be obtained by

Algorithm 1
1: n ← P ωm ∼ P( ωm ), Generate a Poisson rv with ωm in (21)  8: Generate n generalized gamma rv's with the same scale 1 − α and random rates 11: X (t m ) ← a m X (t m−1 ) + z m . 14: end for numerical inversion. Of course, in this study we are considering the easier symmetric situation which is somehow justified by the above observations. The characterization of the BDLP of a SNTS-OU process is given by the following proposition, where we take p = 1 then the tempering function is completely monotone. We will comment on this last assumption afterward.
Proposition 3.2 Under the assumption of Proposition 3.1, but with p = 1 and then dropping the dependence on p from for simplicity, the BDLP L X (·) is the sum of two independent Lévy processes L 1 (·) and L 2 (·). L 1 (1) ∼ SN T S(α, 2 b α c, Q, σ ), hence it is a SNTS process, and namely, a compound Poisson process with intensity and the independent jump sizes have pdf  Cont and Tankov (2004), the Lévy densityν X (x) of a SN T S(α, c, Q, σ ) distribution can be represented as where we have used 6.561.16 page 676 in Gradshteyn and Ryzhik (2007), The pdf of the jumps is then that concludes the proof.
We remark that Proposition 3.2 and Theorem 3.1 in Barndorff-Nielsen (1998) are in agreement when one takes α = 1/2 and the symmetric distribution. Of course, our notation is different and it is a well known fact that the sum of standard normal rv's is chi-squared distributed which in turn is a (1/2, 1) law. As a consequence of Proposition 3.2 an SNTS-OU process can be alternatively represented as and accordingly its transition density can also be derived from (29). Once again for p = 1, and Q(ds) = δ β (ds) the law ofZ (t) has been studied in Sabino (2022), Proposition 3.1, and we can rely on the fact that where τ n , n ∈ N are the jumping times of N (·) and U n ∼ U([0, 1]), with H n , U n , τ n , n ∈ N all mutually independent. The last equality in distribution is a consequence of the observations in Lawrance (1980). All these results can be summarized in the following corollary

Corollary 3.3 Under the assumptions of Proposition 3.2 but with Q(ds) = δ β (ds), X(t) can be represented as the sum in distribution of four independent rv's
where X 0 = X (0) P-a.s.,

is a compound Poisson rv whereP ω is an independent Poisson rv with parameter
and J k , k > 0 are iid rv's with density This corollary then, entails an alternative procedure, summarized in Algorithm 2, for the simulation of an SNTS-OU process valid for p = 1 and Q(ds) = δ β (ds). Of course, it is not complicated to relax this last assumption, on the other hand, in our numerical illustrations we will only consider the case of exponentially tilted TS laws. Finally, Remark 1 is applicable to the corollary as well.
It is straightforward to notice that the numbers of steps of Algorithm 1 is much lower than Algorithm 2 that means that one can expect that the former one will be faster, as will be shown in the section with the numerical illustration. Moreover, as a side-product the two procedures provide easier ways to perform the parameter estimation from real data using Monte Carlo (MC) based techniques compared to those illustrated in Barndorff-Nielsen (1998) for the NIG process of OU type.

The chf and cgf of the transition law
The chf, lch and cgf of the transition law of a SNTS process are quantities which are very helpful for instance, for the pricing of financial derivatives or for the derivation of the risk-neutral conditions. Proposition 3.4 Under the assumptions of Proposition 3.1 but with p = 1, denoting ϕ i (u) and ψ i (u) the chf and the lch of M i , i = 1, 2 respectively, the chf ϕ X (u, t) and the lch ψ X (u, t) of X (t) in (17) are where ϕ J (u) is the chf of J .
Proof Equation (36) is a consequence of (15) and (18), whereas ψ 1 (u) is the lch of a TS subordinator, see Küchler and Tappe (2013). From Proposition 3.1 we know that M 2 is a compound Poisson rv therefore (38) expresses its lch in terms of the chf of the jumps sizes ϕ J (u) which in turn is Remark 3 As far as we are aware of, the chf of a T S(α, c, Q, p) distribution is known for p = 1 and p = 2 only. In this last case the chf can be written in terms of confluent hypergeometric functions (see Bianchi et al. (2017)). Moreover, for p = 1 it is not possible to write the Lévy density of the stationary law in terms of the modified Bessel function of the second type. These observations explain why we presented some propositions under the less general condition p = 1 only. On the other hand, we conjecture that for p > 1 Proposition 3.2 should be still valid after replacing the gamma laws with the generalized gamma laws and some other small changes.
Corollary 3.5 Taking p = 1 and Q(ds) = δ β (ds), the cgf m X (s, t) of X (t) in (17) where m Z X (s, t), m 1 (s) and m 2 (s) are the cgf of Z X (t), M 1 and M 2 , respectively.

Numerical illustrations
In this section, we analyze the performance and the effectiveness of the algorithms for the simulation of SNTS-OU processes and we present possible financial applications, namely derivative pricing. We also assume, p = 1 and Q(ds) = δ β (ds) and take CTS subordinators with an unbiased clock, E [S(t)] = t in which case they form a one-parameter family with Finally, all the numerical experiments in this paper have been conducted using Python with a 64-bit Intel Core i5-6300U CPU, 8GB.

Path simulations
The performance of the algorithms is measured in terms of the absolute percentage error relatively to the second and fourth cumulants denoted err % and defined as err % = |true value − estimated value| true value namely, the absolute percentage difference between the theoretical values of the second and fourth cumulants and their MC estimation. From Equation (10) one can calculate the cumulants c X ,k (x 0 , t), k = 1, 2, . . . of X (t) for X (0) = x 0 knowing the cumulants cX ,k of the stationary SNTS law, which in turn can be written as (see ) Of course, one could also legitimately start from (9) but that is less straightforward. Table 1 and 2 show the estimated cumulants taking (b, σ, ν) = (1, 0.201, 0.256) with two t's and then varying the number of MC repetition R. Tables 3 also reports the CPU times in seconds using Algorithm 1 and Algorithm 2, respectively. The choice t = 1/360 and t = 1/12 can be interpreted as daily and monthly timesteps, respectively. Apparently, the two different procedures return unbiased cumulants although a relatively high number of repetitions should be used to attain low percentage errors. The speed of convergence seems somehow slower for t = 1/360, our interpretation is that with this value of t the expected number of jumps of the compound Poisson component is low but still non-negligible and accordingly a higher number of repetitions is required to attain the same percentage error with t = 1/12. In contrast, when looking at the computational times in Table 3 the alternatives have a completely different behaviour. As expected, Algorithm 1 is faster because its implementation is slimmer with fewer steps than Algorithm 2. Also one can notice that when α = 0.5 the execution times are lower in both cases because the I G random variate can be drawn with the efficient many-two-one transformation method of Michael et al. (1976). In comparison, the difference in computational time is smaller for t = 1/360 mainly because the time is dominated by the generation of TS random variates whereas, for t = 1/12 the cost of simulating from the compound Poisson plays a more relevant role. Based on the above observations, one can conclude that Algorithm 1 is the more preferable solution. Finally, Fig. 1 plots the sample trajectories with Algorithm 2 and varying α.

Financial applications
We consider two financial applications: the pricing of a strip of daily European options in energy or commodity markets and the modeling of future markets. Table 1 Cumulants, c X ,2 is multiplied by 10 4 whereas, c X ,4 by 10 5 . (b, σ, ν) = (1, 0.201, 0.256), t = 1/360   A daily strip of M call options with maturity T and strike K is a contract with payoff where P(·) denotes the spot price process of a certain commodity, for instance gas of power. Such a contract normally encompasses monthly, quarterly and yearly maturities, but is not very liquid and is generally offered by brokers.
The spot price of power or gas and in general of commodities exhibit mean-reversion, seasonality and spikes, which are hard to be modeled in a pure Gaussian world. Different approaches have been investigated in order to somehow extend the classical Gaussian framework introduced in Lucia and Schwartz (Jan 2002) and Schwartz and Smith (2000). Among others, Cartea andFigueroa (2005), Meyer-Brandis andTankov (2008) and Sabino and Cufaro Petroni (2021a) have studied mean-reverting jump-diffusions to model sudden spikes, whereas . Cummins et al. (2017Cummins et al. ( , 2018, and Sabino ( , 2020 have considered VG, CGMY and OU-SNTS processes to price power or gas derivative contracts. We assume that the evolution of the spot dynamics is a onefactor model driven by a SNTS-OU process, this assumption is somehow similar to Benth and Šaltyté Benth (2004) and Benth et al. (2007) where the authors instead take a OU-NIG. Accordingly we have where h(t) is a deterministic function, F(0, t) is the forward curve derived from quoted products and reflects the seasonality, whereas X (·) is the SNTS-OU process in (4) under the assumption that p = 1 and Q(ds) = δ β (ds). For this financial application the chf and the cgf obtained in Sect. 3.3 play a key role for finding the risk-neutral conditions and the pricing relying on the FFT-based technique of Carr and Madan (1999). Indeed, as also observed in Lemma 3.1 in Hambly et al. (2009) the risk-neutral conditions are met when the deterministic function h(t) is given by where m X (s, t) is the cgf of X (t), therefore, assuming for simplicity X (0) = 0, we have with m i (·), i = 1, 2 are defined in Corollary 3.5.    Table 3 Computational times in seconds for R = 10 6 repetitions (b, σ, ν) = (1, 0.201, 0.256). A 1 and A 2 stand for Algorithm 1 and Algorithm 2, respectively whereas, the columns relative to A 2 are the multiplicative factors compared to the CPU's of Algorithm 1 As a consequence of Proposition 3.4, we can make use of the explicit form of the chf of the log P(t) = log F(0, t) + h(t) + X (t) to compute the price of such a strip of calls using the FFT-based technique of Carr and Madan (1999).
For this specific example we consider a set with realistic parameters, taken from Sabino (2022), (b, ν, σ ) = (10, 0.7, 0.2) and let α vary. In addition, in order to better highlight the dependency on the maturity and α, we take a flat forward curve F(0, t) = 20, indeed the seasonality changes the moneyness of the strip of call options and could shadow some features; therefore we also assume an at-the-money strike K = 20. Table 4 shows the numerical results relatively to different α's and different maturities spanning from one month to one year. Table 4 also compares the results obtained with the FFT method to those with the MC method with 10 5 simulations. As far as this last method is concerned, we also report the standard errors defined as the sample standard deviation divided by the square root of the relative number of simulations. As expected, we observe that the FFT and the MC approaches return comparable results, of course the FFT method is known to be faster and is then preferable.
As second application we consider forward contracts in energy markets, sometimes called swaps, whose value at time t, maturity T , t ≤ T ≤ T 1 < T 2 , and with delivery period [T 1 , T 2 ] is denoted by F(t, T 1 , T 2 ). Often the evolution of these contracts is modeled by exponential Lévy or additive processes as done for instance in Kiesel et al. (2009) or Goutte et al. (2014) with one or more driving factors. A discussion on the differences of the available approaches is beyond the scope of this study, what is nevertheless relevant is the fact that European call and put options written on these contracts exhibit a volatility smile and the so-called Samuelson effect which cannot be described by a pure Lévy or exponential Lévy mode. One has to resort then to additive processes. In the same vein of the setting of Piccirilli et al. (2021) but for simplicity with one-factor only, one can model F(t, T 1 , T 2 ) as where m(t, T 1 , T 2 ) is a deterministic functions that ensures that F(·, T 1 , T 2 ) is a martingale. Moreover, is another deterministic functions meant to capture the Samuelson effect (see also Jaeck and Lautier (2016)). Of course, such dynamics are related to the process studied in Section 3, because, after some algebra, it turns out hence the chf and in the particular, the simulation procedure of the skeleton of the additive process X (·, T 1 , T 2 ) can  be derived from those of Z X (·). Accordingly one can adopt FFT or MC-based techniques to price derivative contracts and in addition, can rely on simulation methods which are very useful for risk managers to derive information about the distribution of the traded portfolios and a view on their risk profile. To this end, it is fundamental to dispose of an exact and non-biased simulation scheme with sufficiently fast computational speed that can be run with different parameter settings and market conditions. Finally, Fig. 2 shows the evolution according to the model (4.2) of a swap contract with maturity in six months, T 1 =1/2, and delivery period over a quarter, T 2 − T 1 = 1/4, assuming (b, σ, ν) = (20, 0.201, 0.256), with initial price F(0, T 1 , T 2 ) = 20 and different values of α.

Conclusions
In this paper we have studied the properties of the OU process with a symmetric normal tempered stable stationary law, named symmetric NTS-OU process, with main focus on the transition law and the backward driving Lévy processes. Our findings complement the existing state of affairs because the characterization of a NTS-OU was still missing. On the other hand the full asymmetric case has not here been covered and will be the object of future studies. We have found that the transition distribution and the backward driving Lévy process can be represented in terms of the sum of independent random variables or processes, where two terms are always of normal tempered stable and of compound Poisson type. This result has two important implications: the design of fast simulation methods and the derivation of the characteristic function. We then present two applications to energy or commodity markets, namely the pricing of a strip of European call options and the modeling of forward markets.
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://creativecomm ons.org/licenses/by/4.0/. Schwartz, P., Smith, J.E.: Short-term variations and long-term dynamics in commodity prices. Manag. Sci. 46(7), 893-911 (2000) Zhang, B., Oosterlee, C.W.: Efficient pricing of European-style asian options under exponential Lévy processes based on fourier cosine expansions. SIAM J. Financial Math. 4(1), 399-426 (2013) Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.