Exact Simulation of Variance Gamma related OU processes: Application to the Pricing of Energy Derivatives
aa r X i v : . [ q -f i n . C P ] A p r Exact Simulation of VarianceGamma related OU processes:Application to the Pricing ofEnergy Derivatives. ∗ Piergiacomo
Sabino † Quantitative Methods, E.ON SEBr¨usseler Platz 1, 45131 Essen, Germany
Abstract
In this study we define a three-step procedure to relate the self-decomposabilityof the stationary law of a generalized Ornstein-Uhlenbeck process to the law ofthe increments of such processes.Based on this procedure and the results of Qu et al. [36], we derive the exactsimulation, without numerical inversion, of the skeleton of a Variance Gamma,and of a symmetric Variance Gamma driven Ornstein-Uhlenbeck process. Ex-tensive numerical experiments are reported to demonstrate the accuracy andefficiency of our algorithms.These results are instrumental to simulate the spot price dynamics in energymarkets and to price Asian options and gas storages by Monte Carlo simulationsin a framework similar to the one discussed in Cummins et al. [18, 19].
Keywords
Monte Carlo, Exact simulation, Non-Gaussian Ornstein-Uhlenbeck(OU) processes, OU-Variance-Gamma processes, Energy Markets, Energy Deriva-tives.
The modeling based on non-Gaussian Ornstein-Uhlenbeck (OU) processes has re-ceived a considerable attention in the recent literature in an attempt to accommodatefeatures such as jumps, heavy tails and asymmetry which are well evident in real phe-nomena. For instance, with regards to financial and econometric applications, energymarkets, and commodity markets in general, exhibit mean-reversion, seasonality andsudden spikes; mean-reversion in particular, cannot be captured by ordinary L´evyprocesses. ∗ The views, opinions, positions or strategies expressed in this work are those of the author anddo not represent the views, opinions and strategies of, and should not be attributed to E.ON SE. † [email protected] he availability of simulation techniques of easy implementation is important foranalysis, validation and estimation purposes. Indeed, direct likelihood analysis isoften impracticable for these models, whereas, Monte Carlo (MC) based techniquesand generalized method of moments (GMM) approaches can be a viable route toestimate the model parameters.As observed in a series of papers by Barndorff-Nielsen [3], Barndorff-Nielsen etal. [2], Barndorff-Nielsen and Shephard [4, 5], the concept of self-decomposability (seeSato [38] and Cufaro Petroni [14]) plays an essential role in the theory of generalizedOU-processes. In this paper we define a simple three-steps procedure to determinethe characteristic function, and therefore the cumulant function, of a OU process andits relation to the characteristic function of what we name the a-reminder of a self-decomposable law. This machinery demonstrates to be a powerful tool to determinesimulation algorithms for generalized OU-processes as well as to simplify alreadyexisting proofs (see for instance Qu et al. [36] and Bianchi et al. [10]).Relying on the results of Qu et al. [36], the main contribution of this article isthe development of exact simulation schemes to generate the skeleton of VarianceGamma (VG) driven OU processes (OU-VG) discussed in Cummins et al. [18, 19].The extensive simulation experiments show that our algorithms are efficient andaccurate therefore, suitable for concrete applications.To this end, the modeling of energy markets with non-Gaussian OU processes hasbeen discussed, among others, in Benth et al. [7], Meyer-Brandis and Tankov [34] andrecently in Benth and Pircalabu [8] in the context of modeling wind power futures.Compared to mean-reverting jump-diffusion models (see for instance, Cartea andFigueroa [12] and Kjaer [28]) these models exhibit the competitive advantage ofhaving less parameters.We illustrate the applicability of our schemes in the pricing Asian options andgas storages by MC simulation using market dynamics similar to those discussedin Cummins et al. [18, 19]. Once more, our algorithms demonstrate to be efficientand reasonably fast to compute the fair values of such energy derivatives. AlthoughMC methods are not as fast as other numerical techniques as FFT and quantizationmethods (see for instance Jaimungal and Surkov [26] and Bardeau et al. [1]), theynevertheless, give the possibility to compute different quantiles of the price distribu-tion and are independent on contract payoffs.The remainder of the article is organized as follows. Section 2 recalls the prop-erties of generalized OU processes and introduces the conceptual procedure whichwe will use in order to develop the simulation schemes. In Section 3 we derive thecharacteristic function of the law of the increments of OU-VG processes and developsimulation schemes for the skeleton of such processes. In this section, we also demon-strate the effectiveness of our algorithms through extensive numerical experiments.Section 4 illustrates some financial applications: we consider the pricing of Asianoptions by MC simulations in a 2-factor market driven by the sum of a standard VGprocess and a OU-VG process then, we consider the pricing of gas storages using aone-factor spot dynamics similar to the setting discussed in Cummins et al. [18, 19].Finally, Section 5 concludes the paper with an overview of future inquiries and furtherpossible applications. 2 Preliminaries
Following Barndorff-Nielsen and Shephard [4], we consider a L´evy process Z ( t ) andthe generalized OU process defined by the SDE dX ( t ) = − kX ( t ) dt + dZ ( t ) X (0) = X P - a.s. k > . (1)with solution X ( t ) = X (0) e − kt + Z t e k ( t − v ) dZ ( v ) . (2)Here Z ( t ) is called the Backward Driving L´evy Process ( BDLP ), and we will adoptthe following notation: if D is the stationary law of X ( t ), we will say that X ( t )is a D - OU process; if on the other hand, Z (1) is distributed according to the id (infinitely divisible) law e D , then we will say that X ( t ) is an OU - e D process. Nowa well known result (see for instance Cont and Tankov [13] or Sato [38]) is that, agiven one-dimensional distribution D always is the stationary law of a suitable OU - e D process if and only if D is self-decomposable.We recall that a law with probability density ( pdf ) f ( x ) and characteristic function( chf ) ϕ ( u ) is said to be self-decomposable ( sd ) (see Sato [38] or Cufaro Petroni [14])when for every 0 < a < pdf g a ( x ) and chf χ a ( u ) suchthat ϕ ( u ) = ϕ ( au ) χ a ( u ) (3)We will accordingly say that a random variable ( rv ) X with pdf f ( x ) and chf ϕ ( u ) is sd when its law is sd : looking at the definition, this means that for every 0 < a < independent rv ’s, Y (with the same law of X ) and Z a (herecalled a -remainder ), with pdf g a ( x ) and chf χ a ( u ) such that X d = aY + Z a P - a.s. (4)As observed in Barndorff [2], X ( t ) is stationary if and only if the chf of the sd stationary law ϕ X ( u ) is of the form ϕ X ( u ) = ϕ X ( u e − kt ) χ ( u, t ), where χ ( u, t ) denotesthe chf of the second term of Equation (2). Defining the cumulant function of a rv Y as κ Y ( u ) = log E (cid:2) e uY (cid:3) , it turns out that there is precise relation between thecumulant function of the stationary distribution ¯ κ X ( u ), that of Z (1), denoted κ Z ( u ),and that of the second term of Equation (2), denoted ̺ X ( u, t ) (see also Taufer andLeonenko [41] and Schoutens [39]).¯ κ X ( u ) = Z + ∞ κ Z ( u e − kt ) ds (5) ̺ X ( u, t ) = ¯ κ X ( u ) − ¯ κ X ( u e − kt ) (6)The last equation means that the law of the second term of Equation (2) coincideswith that of the a -remainder of the law of the stationary distribution if one takes a = e − kt . A similar observation was also mentioned in Gaver and Lewis [22], Lawrence [30]and later in Wolfe [42] in the context of first order auto-regressive processes X n = ρX n − + ǫ n , in case 0 ≤ ρ <
1, that are the discrete-time equivalent of OU processes.This facts give a useful machinery to determine the chf or the cumulant of OUprocesses which of course, can be used to find simulation algorithms.3
Find the cumulant function of the stationary distribution given the BDLP. • Find the cumulant function of the a -remainder of the marginal distribution. • Set a = e − kt .On the other hand, based on the observations above, the sequential generation of theskeleton of X ( t ) on a time grid t , . . . t M consists in finding a simulation algorithmfor the a -remainder of the stationary law assuming at each step a i = e − k ( t i − t i − ) , i =1 , . . . , M . Hereafter, without loss of generality, we will assume an equally-spaced timegrid with ∆ t = t i − t i − , ∀ i = 1 , . . . , M .Finally, because the cumulant function κ X ( u, t ) = log E (cid:2) e uX ( t ) (cid:3) can also be writ-ten in terms of the cumulant function κ Z ( u ) of Z (1) as (see Cont and Tankov [13]Lemma 15.1) κ X ( u, t ) = uX (0) e − kt + Z t κ Z (cid:0) ue − k ( t − v ) (cid:1) dv = uX (0) e − kt + ̺ X ( u, t ) , (7)one can relate the cumulants κ X,n of X ( t ) to the cumulants κ Z,n of Z (1) E [ X ( t )] = X (0) e − kt + κ Z, k (cid:0) − e − kt (cid:1) (8) κ X,n = κ Z,n n k (cid:0) − e − nkt (cid:1) , n ≥ The VG process, introduced in Madan and Seneta [33], can be seen as a Brown-ian Motion (BM) where the clock ticks with a random time described by a gammasubordinator G ( t ).We recall that the gamma law Γ( α, β ) is famously sd (see Grigelionis[24]) and hasthe following pdf and chf f α,β ( x ) = β α x α − Γ( α ) e − βx x ≥ φ Γ ( u ) = (cid:18) ββ − iu (cid:19) α where Γ( · ) is the Euler gamma function, α > β > shape and rate parameters, respectively.The gamma process G ( t, α, β ) is a continuous-time process with stationary, inde-pendent gamma increments such that for any h > G ( t + h, α, β ) − G ( t, α, β ) ∼ Γ( αh, β ) . (10)therefore, it is a subordinator (see Sato [38]). In order to guarantee that the stochasticclock G ( t ) is an unbiased reflection of calendar time (see Joshi [27]) we need to set4 [ G ( t )] = t . The law of the increments of G ( t ) now depends on one parameter only ν = α = β . Denoting now G ( t, ν ) the gamma process with the above parametersrestriction, the VG process is defined as follows: V ( t ) = θG ( t, ν ) + σW ( G ( t, ν )) , (11)with characteristic exponent ( che ) ψ V G ( u ): ψ V G ( u ) = log E (cid:2) e iuV (1) (cid:3) = − ν log (cid:18) − iuθν + u σ ν (cid:19) . (12)for σ > θ ∈ R constants.Since the VG is a process of finite variation, it can be written as difference of twoincreasing gamma processes V ( t ) = γ p ( t, µ p , ν p ) − γ n ( t, µ n , ν n )with µ p , ν p , µ n , ν u , ν satisfying the following equations µ p = 12 r θ + 2 σ ν + θ ,µ n = 12 r θ + 2 σ ν − θ ,ν p = µ p ν,ν n = µ n ν. The che can then be rewritten as ψ V G ( u ) = − ν log (cid:18) − iu ν p µ p (cid:19) − ν log (cid:18) iu ν n µ n (cid:19) = ψ Γ p ( u ) + ψ Γ n ( − u ) , (13)where ψ Γ p ( u ) and ψ Γ n ( u ) are the che ’s of a Γ( ν , µ p ν p ) and a Γ( ν , µ n ν n ) law, respectively(therefore of the difference of two independent gamma-distributed rv ’s). When θ = 0- in case of a symmetric VG (SVG) - it simplifies to ψ SV G ( u ) = − ν log (cid:18) u σ ν (cid:19) = ψ Γ ( u ) − ψ Γ ( − u ) . (14)where now ψ Γ ( u ) denotes the che of a Γ( v , σ ν ) law. Using the machinery illustrated inSection 2, we can calculate the cumulant function ¯ κ X ( u ) of the stationary distributionof the VG driven OU process X ( t ) = X (0) e − kt + Z t e − k ( t − v ) dV ( v ) (15)and the cumulant function ̺ X ( u, t ). Proposition 3.1.
The cumulant function ¯ κ X ( u ) and ̺ X ( u, t ) of a OU-VG processare given by ¯ κ X ( u ) = 1 kν (cid:18) Li (cid:18) u µ p ν p (cid:19) + Li (cid:18) − u µ n ν n (cid:19)(cid:19) . (16) ̺ X ( u, t ) = 1 kν (cid:18) Li (cid:18) u µ p ν p (cid:19) − Li (cid:18) u µ p e − kt ν p (cid:19) + Li (cid:18) − u µ n ν n (cid:19) + Li (cid:18) − u µ n e − kt ν n (cid:19)(cid:19) (17)5 roof. Denoting κ V G ( u ) = ψ V G ( − iu ), we have¯ κ X ( u ) = Z + ∞ κ V G ( ue − ks ) ds = Z + ∞ κ Γ p ( ue − ks ) ds + Z + ∞ κ Γ n ( − ue − ks ) ds where κ Γ p ( u ) = ψ Γ p ( − iu ) and κ Γ n ( u ) = ψ Γ n ( − iu ). With the change of variable x = e − ks we have¯ κ X ( u ) = − kν Z log(1 − u x ν p µ p ) x dx + Z log(1 + u x ν n µ n ) x dx ! Of course, the last two terms are the cumulant functions of a stationary OU − Γprocesses and of its negative counterpart (see Qu et al. [36] and Table 2 in Barndorff-Nielsen and Shephard [5]) whichx can be written in terms of the dilogarithmicSpencer’s function Li ( z ) = − R z − y ) y dy, z ∈ C (see Gradshteyn and Ryzhik [23])¯ κ X ( u ) = 1 kν (cid:18) Li (cid:18) u µ p ν p (cid:19) + Li (cid:18) − u µ n ν n (cid:19)(cid:19) . Hence ̺ X ( u, t ) = − kν Z e − kt log (cid:16) − u x ν p µ p (cid:17) x dx + Z e − kt log (cid:16) u x ν n µ n (cid:17) x dx = 1 kν (cid:18) Li (cid:18) u µ p ν p (cid:19) − Li (cid:18) u µ p e − kt ν p (cid:19) + Li (cid:18) − u µ n ν n (cid:19) − Li (cid:18) − u µ n e − kt ν n (cid:19)(cid:19) that concludes the proof.By simply setting θ = 0 we retrieve the cumulant function relative to the OU-SVGprocess of Cummins et al. [18].¯ κ X ( u ) = 12 kν Li (cid:18) u σ ν (cid:19) , (18)and ̺ X ( u, t ) = 12 kν (cid:18) Li (cid:18) u σ ν (cid:19) − Li (cid:18) u σ ν e − kt (cid:19)(cid:19) (19) From the results of the previous section we can conclude that the simulation of aOU-VG process consists in the repetition of the simulation a OU-Γ process two timesand then take the difference. To this end, Qu et al. [36] found that a rv Y withcumulant function ̺ Y ( u, ∆ t ) = − αk Z e − k ∆ t log (cid:16) u xβ (cid:17) x dx can be decomposed into the sum of a gamma-distributed rv Y ∼ Γ( α ∆ t, β e k ∆ t )and a compound Poisson process Y = P Mm =1 J m with intensity λ = αk ∆ t and6xponentially-distributed jumps J m with random rate β e k ∆ t √ U and U ∼ U ([0 , t , t , . . . , t M with step ∆ t consists in nothing less than simulating a Y -like rv two times at each step and then taking the difference as illustrated in Algorithm 1.On the other hand, the simulation procedure of a symmetric VG can be simplified Algorithm 1 for m = 1 , . . . , M do Generate G p ∼ Γ( ∆ tν , µ p ν p e k ∆ t ) Generate G n ∼ Γ( ∆ tν , µ n ν n e k ∆ t ). Generate r ∼ P ( k ∆ t ν ) ⊲ Poisson rv with intensity k ∆ t ν Generate s ∼ P ( k ∆ t ν ) Generate r iid uniform rv ’s u = ( u , . . . , u r ) ∼ U ([0 , r ). Generate s iid uniform rv ’s v = ( v , . . . , v s ) ∼ U ([0 , s ). β p,i ← µ p ν p e k ∆ t √ u i , i = 1 , . . . , r . β n,j ← µ n ν n e k ∆ t √ v j , j = 1 , . . . , s . Generate r iid J p,i ∼ E ( β p,i ) , i = 1 , . . . , r , ⊲ Exponential rv ’s with rate β p,i Generate s iid J n,j ∼ E ( β n,j ) , j = 1 , . . . , s , ⊲ Exponential rv ’s with rate β n,j X ( t m ) ← X ( t m − ) e − k ∆ t + G p − G n + P ri =1 J p,i − P sj =1 J n,j . end for observing that for θ = 0, µ p ν p = µ n ν n = σ q ν and that the difference C = C − C of two iid compound Poisson processes C = P N n U n and C = P N n D n with intensity λ and with exponentially distributed jumps has the same law of P Nn ( U n − D n ) with in-tensity 2 λ . It is well known that the difference of exponentially distributed rv ’s withthe same scale parameter µ is distributed according to a central Laplace law La ( µ ).Such a rv can be efficiently generated using the inverse transformation method (seeDevroye [20]) observing that the inverse of the cumulative distribution is F − L ( y ) = − µ sign( y − .
5) ln(1 − | y − . | ) . Based on these observations, the steps of the sequential simulation of the skeleton ofa symmetric OU-SVG process are summarized in Algorithm 2
Algorithm 2 for m = 1 , . . . , M do Generate G p ∼ Γ( ∆ tν , σ q ν e k ∆ t ) Generate G n ∼ Γ( ∆ tν , σ q ν e k ∆ t ). Generate r ∼ P ( k ∆ t ν ) ⊲ Poisson rv with intensity k ∆ t ν Generate r iid uniform rv ’s u = ( u , . . . , u r ) ∼ U ([0 , r ). µ i ← σ p ν e − k ∆ t √ u i , i = 1 , . . . , r . Generate r iid J i ∼ La ( µ i ) , i = 1 , . . . , r , ⊲ Laplace rv ’s with parameter µ i X ( t m ) ← X ( t m − ) e − k ∆ t + G p − G n + P ri =1 J i . end for .2 Numerical Experiments In this section, we illustrate the performance and effectiveness of our algorithmsthrough extensive numerical experiments. All the simulation experiments in thepresent paper have been conducted using
MATLAB R2019a with a 64-bit Intel Corei5-6300U CPU, 8GB . As an additional validation, the comparisons of the simulationcomputational times have also been performed with R and Python leading to the sameconclusions.The numerical validation and tests for our algorithms are based on the com-parison to the true expected value, variance, skewness and kurtosis of the OU-VGprocess. Because of Equations (8)and (9) and the relation between the cumulants,these quantities relative to X ( t + ∆ t ) = X ( t ) e − k ∆ t + R ∆ t e − k ( t − v ) V ( v ) are E [ X ( t + ∆ t )] = a X ( t ) + (1 − a ) θk V [ X ( t + ∆ t )] = (1 − a ) σ + θ ν k Skew [ X ( t + ∆ t )] = 2 √ k − a )(1 − a ) / θ ν + 3 σ θν ( σ + θ ν ) / Kurt [ X ( t + ∆ t )] = k (1 + a )(1 − a ) × σ ν + 12 σ θ ν + 6 θ ν ( σ + θ ν ) + 3where a = e − k ∆ t whereas, the symmetric case is simply obtained with θ = 0.Table 1 reports the CPU times in seconds and compares the MC estimated valuesof the true E [ X ( T )], V [ X ( T )], Skew [ X ( T )] and Kurt [ X ( T )]. The values at thetop of the table are obtained with a single time step T = ∆ t = 1 / /
5. Varyingthe number of simulations N S , we can conclude that our algorithm is efficient andconvergent, although it seems that at least N S = 10 simulations is required to achievea good estimate. However, although the algorithm provides an exact simulation ofa OU-VG process, the generation of an entire trajectory, especially over a time gridwith several points is not extremely fast compared to the simulation of other OUprocesses (see for instance Cufaro Petroni and Sabino [16, 17]). Figure 1a shows asample trajectory using a time grid of 365 points that is a quite common choice infinancial applications relative to the pricing of a one year contract.We conclude this section illustrating the results of the numerical experimentsrelative to a OU-SVG process. E [ X ( T )] = X (0) e − kT and the skewness is zerotherefore, in Table 2 we show the CPU times in seconds and the MC estimated valuesof the true V [ X ( T )] and Kurt [ X ( T )] only. The simulation has been conductedusing the same time grid as the previous case whereas, the process parameters arethose presented in Cummins et al. [18]; Figure 1b shows a sample path of such aprocesses over a time grid of 365 points. Once again, the simulation algorithm isconvergent and captures the true values of the benchmarks quite well. Moreover,Algorithm 2 is faster than Algorithm 1 of almost a factor 2 and provides an efficientsolution when the drift θ can be neglected. The relative codes are available at https://github.com/piergiacomo75/OUVarianceGamma M = 365, ∆ t = 1 / Number of steps -0.1-0.0500.050.10.150.20.250.30.35
Trajectory of a OU-VG process (a) OU-VG X (0) = 0, θ = 0 . k = 0 . ν = 0 . σ = 0 . Number of steps -0.06-0.05-0.04-0.03-0.02-0.0100.010.020.030.04
Trajectory of a OU-Symmetric-VG process (b) OU-SVG with X (0) = 0, k = 0 . ν =0 . σ = 0 . E [ X ( T )] = 0 . V [ X ( T )] = 0 . Skew [ X ( T )] = 0 . Kurt [ X ( T )] = 4 . T = 1 / , ∆ t = 1 / N S CPU MC error % MC error % MC error % MC error %2500 0 .
05 0 . .
14 0 . .
68 0 .
537 1 .
52 4 .
786 2 . .
20 0 . .
17 0 . .
57 0 .
502 5 .
08 4 .
637 1 . .
77 0 . .
31 0 . .
07 0 .
529 0 .
07 4 .
703 0 . .
06 0 . .
82 0 . .
48 0 .
535 1 .
17 4 .
666 0 . .
32 0 . .
19 0 . .
02 0 .
522 1 .
39 4 .
646 0 . .
39 0 . .
18 0 . .
11 0 .
527 0 .
38 4 .
678 0 . E [ X ( T )] = 0 . V [ X ( T )] = 0 . Skew [ X ( T )] = 0 . Kurt [ X ( T )] = 3 . T = 1 , ∆ t = 1 / .
28 0 . .
31 0 . .
97 0 .
214 9 .
81 3 .
220 3 . .
02 0 . .
18 0 . .
83 0 .
248 4 .
31 3 .
332 0 . .
02 0 . .
28 0 . .
15 0 .
236 0 .
70 3 .
315 0 . .
03 0 . .
09 0 . .
10 0 .
230 3 .
42 3 .
316 0 . .
46 0 . .
20 0 . .
02 0 .
239 0 .
55 3 .
341 0 . .
51 0 . .
06 0 . .
10 0 .
238 0 .
30 3 .
338 0 . Table 1:
CPU times in seconds and comparison among the true E [ X ( T )], V [ X ( T )], Skew [ X ( T )]and Kurt [ X ( T )] of a OU-VG process with ( k, θ, ν, σ, X (0)) = (0 . , . , . , . ,
0) and their relativeestimated values with N S MC scenarios. T = 1 / , ∆ t = 1 / T = 1 , ∆ t = 1 / V [ X ( T )] = 0 . Kurt [ X ( T )] = 6 . V [ X ( T )] = 0 . Kurt [ X ( T )] = 3 . N S CPU MC error % MC error % CPU MC error % MC error %2500 0 .
04 0 . .
96 6 .
89 0 .
67 0 .
20 0 . .
69 3 .
341 11 . .
12 0 . .
73 6 .
66 2 .
62 0 .
64 0 . .
07 3 .
748 0 . .
53 0 . .
66 6 .
96 1 .
79 2 .
50 0 . .
29 3 .
853 1 . .
04 0 . .
73 6 .
91 1 .
03 10 .
16 0 . .
71 3 .
717 1 . .
06 0 . .
18 6 .
78 0 .
88 39 .
98 0 . .
38 3 .
793 0 . .
10 0 . .
32 6 .
85 0 .
16 159 .
84 0 . .
04 3 .
778 0 . Table 2:
CPU times in seconds and comparison among the true V [ X ( T )] and Kurt [ X ( T )]of a OU-Symmetric VG process with ( k, ν, σ, X (0)) = (0 . , . , . ,
0) and their relativeestimated values with N S MC scenarios
Pricing derivative contracts or energy facilities is often accomplished using MC meth-ods; for this purpose is therefore, necessary to rely on efficient and eventually, fastpath-generation techniques. On the other hand, the day-ahead (also called spot) priceof power or gas and in general of commodities exhibit mean-reversion, seasonality9nd spikes, this last feature is particularly difficult to be captured in a pure Gaussianworld. Different approaches have been investigated in order to somehow extend theclassical Gaussian framework introduced in Lucia and Schwarz [32] and Schwartz andSmith [40]. Among others, Cartea and Figueroa [12], Kjaer [28], Meyer-Brandis andP. Tankov [34] have studied mean-reverting jump-diffusions to model sudden spikes,whereas, Benth et al. [7] and Benth and A. Pircalabu [8] have considered differentnon-Gaussian OU processes in order to price power or wind derivative contracts. Re-cently, Cummins et al. [18, 19] have addresses the pricing of gas storages via FFT ina market driven by a OU-SVG process.In the following subsections we illustrate the effectiveness of the algorithms pre-sented in subsection 3.1 when applied to the pricing of Asian options and of gasstorages using market models similar to those considered in Cummins et al. [18, 19].Of course, MC methods are known to be sometimes slower than FFT and othertechniques, nevertheless, they provide a view on the distribution of the potentialcash-flows of derivative contracts giving a precious information to risk managers orto trading units. The calibration and in general, the parameter estimation of OU-VGprocesses is not the aim of the paper. However, as observed in Wolfe [42], discrete first-order autoregressive processes are embedded into continuous OU processes, thereforeone could use the generalized method of moments (GMM) to derive Yule-Walker-likeequations and estimate the model parameters from historical data.
In this section we assume that the spot price of a gas market is driven by the following2-factors process S ( t ) = F (0 , t ) e h ( t )+ X ( t )+ X ( t ) = F (0 , t ) e h ( t )+ H ( t ) (20)where h ( t ) is a deterministic function, F (0 , t ) is the forward curve and X ( t ) is a OU-VG process with parameters ( k, θ , ν , σ ). In contrast to Cummins et al. [19], we adda second independent VG process X ( t ) with parameters ( θ , ν , σ ) to capture thelong-term behavior. Using the risk-neutral arguments of the Lemma 3.1 in Hamblyet al. [25], the deterministic function h ( t ) consistent with forward curve is h ( t ) = − κ H (1 , t ) . (21)where κ H ( u, t ) is the cumulant function of the process H ( t ) at time t , then becauseof Equations (12) and (17) h ( t ) = − kν (cid:18) Li (cid:18) µ p ν p (cid:19) − Li (cid:18) µ p e − kt ν p (cid:19) + Li (cid:18) − µ n ν n (cid:19) − Li (cid:18) − µ n e − kt ν n (cid:19)(cid:19) − tν log (cid:18) − θν − σ ν (cid:19) , (22)where µ p , ν p , µ n , ν n are relative to the parameters of X ( t ).Finally, we recall that the payoff at maturity T of an Asian option with Europeanstyle and strike price K is A ( T ) = d X i =1 ω i S ( t i ) − K ! + . Trajectory of the 2-Factors Model (a) 2-factors Model
Trajectory of the One Factor Model (b) One factor Model
In our numerical experiments we assume an at-the-money Asian option K = F (0 ,
0) =15 having one year maturity, T = 1, with equal weights ω i = 1 /d and with a flatforward curve. Although, as already mentioned, we do not focus on the parame-ters estimation, the values in Table 3a could be considered realistic because they arebased on the estimations presented in Gardini et al. [21] for X ( t ) and are similar tothose shown in Cummins et al. [18] for X ( t ); Figure 2a shows one sample price-pathgenerated with Algorithm 1 using the market model of Equation (20).Table 3: Asian option in a 2-Factors market dynamicsParameter Value F (0 ,
0) 15 K T d κ . θ . ν . σ . θ . ν . σ . (a) Parameters N S CPU price stdev error %-error1000 7 .
59 1 .
219 2 .
24 0 .
071 5 . .
87 1 .
229 2 .
34 0 .
023 1 . .
03 1 .
227 2 .
26 0 .
016 1 . .
10 1 .
231 2 .
26 0 .
010 0 . .
03 1 .
236 2 .
27 0 . . (b) Results Table 3b shows the estimated prices obtained by MC varying the number ofsimulations N S , along with the overall computational times (CPU) in seconds. Thecolumns stdev and error report the standard deviations of the MC estimator and theerrors around the estimated option prices defined as the standard deviation dividedby √ N S . The results illustrate that our simulation scheme is accurate indeed, theprice of the option converges rapidly and the errors are very small; it seems that 10000simulations are good enough to have a reliable Asian option price. On the other hand,although here combined with the simulation of a standard VG process, the simulation11f a OU-VG process is not extremely fast compared to the one of other generalizedOU-processes (see for instance Cufaro Petroni and Sabino [16, 17]). Nevertheless, itprovides additional information, such as quantiles of the price distribution, that canbe employed to derive risk premia and support decision making and in particular, cangive an insight whether the calibrated parameters imply realistic price trajectories. In contrast to the previous subsection, in the following we consider a one factor gasmarket similar to the model discussed in Cummins et al. [18] S ( t ) = F (0 , t ) e h ( t )+ X ( t ) , (23)where X ( t ) is OU-SVG with parameters ( k, ν, σ ) and once more, because of risk-neutral arguments h ( t ) = − kν (cid:18) Li (cid:18) σ ν (cid:19) − Li (cid:18) σ ν e − kt (cid:19)(cid:19) . (24)We then adopt this market dynamics for the pricing of a fast-churn gas storage.To this end, denote by C ( t ) the volume of a gas storage at time t with C min ≤ C ( t ) ≤ C max . The holder of such an energy asset is faced with a timing problem that consistsin deciding when to inject, to withdraw or to do-nothing.Denoting J ( t, x, c ) the value of a gas storage at time t given S ( t ) = x , C ( t ) = c ,one can write: J ( t, x, c ) = sup u ∈U E (cid:20)Z Tt φ u ( S ( s )) ds + q ( S ( T ) , C ( T )) (cid:12)(cid:12)(cid:12)(cid:12) S ( t ) = x, C ( t ) = c (cid:21) , (25)where U denotes the set of the admissible strategies, u ( t ) ∈ {− , , } is the regimeat time t such that φ − ( S ( t )) = − S ( t ) − K in a in , injection φ ( S ( t )) = − K N , do nothing φ ( S ( t )) = S ( t ) − K out a w withdrawal , (26) a in and a w are the injection and withdrawal rates, K in , K out and K N , respectively,represent the costs of injection, do-nothing and withdrawal, and q takes into accountthe possibility of final penalties. Based on the Bellman recurrence equation (seeBertsekas [9]), one can perform the following backward recursion for i = 1 , . . . , d : J ( t i , x, c ) = sup k ∈{− , , } { φ k S ( t i ) + E [ J ( t i +1 , S ( t i +1 ) , ˜ c k ) | S ( t i ) = x, C ( t i ) = c ] } , i = 1 , . . . , d, (27)where ˜ c − = min( c + a in , C max )˜ c = c ˜ c = min( c − a w , C min ) . (28)A standard approach to price gas storages is a modified version of the Least-SquaresMonte Carlo (LSMC), introduced in Longstaff-Schwartz [31], detailed in Boogert and12e Jong [11]. With this approach, the backward recursion is obtained by defining afinite volume grid of G steps for the admissible capacities c of the plant and then applythe LSMC methodology to the continuation value per volume step. In alternative,one may solve the recursion by adapting the method proposed by Ben-Ameur et al. [6]or might use the quantization method as explained in Bardou et al. [1] or even FFTand Fourier techniques described for instance in Jaimungal and Surkov [26].Finally, we consider a one-year fast-churn storage with the parameters shown inTable 4a such that 20 days are required to fill or empty the storage. Our frameworkis similar to the one presented in Cummins et. al. [18] and we select the parametersof the OU-SVG shown therein. Once again, we assume a flat forward curve thatdoes not change the validity of our experiment because we price a fast-churn storagewhose extrinsic value is dominated by the short-term decisions rather than by theseasonality of the forward curve. In this case, Figure 2b shows one sample price-pathgenerated with Algorithm 2 using the market model of Equation (23).The columns of Table 4b have the same meaning of those presented in the anal-ysis of the Asian option. In addition, the column CPU ∗ reports the computationaltimes relative to the path simulation only, using Algorithm 2. Indeed, the computa-tional cost of the standard LSMC method can be split into a path generation stepand into a stochastic optimization step, with the computational cost of the latterone being independent on the price dynamics and being the dominant factor to theoverall computational time (CPU in Table 4b). In this example, the CPU ∗ ’s areapproximately, one forth of the overall time, and in particular, are almost a half ofthe overall time required to price an Asian option with the same number of simula-tions. Once again, the results show that our methodology is accurate and reasonablyfast. However, although still acceptable, the CPU ∗ ’s are higher that those requiredto simulate other generalized OU processes. For instance, one of the reasons whyCummins et al. [18] investigated the use of OU-SVG was to reduce the number ofmodel parameters compared to the mean-reverting jump-diffusion model discussed inKjaer [28]. In this last case however, Cufaro Petroni and Sabino [16] have presented afast simulation algorithm for mean-reverting jump-diffusion dynamics including thecase of the mean-reverting equivalent of the Kou [29] model. These computationaltimes are faster that those presented in this study because they almost completelycut off the time needed for the path simulation that therefore, becomes negligiblecompared to that of the optimization step.It is worthwhile noticing that the performance of standard LSMC could be im-proved relying on the backward simulation of the price dynamics as explained inPellegrino and Sabino [35] and Sabino [37]. This means that one has to design abackward simulation algorithm for the OU-VG process that will be the objective ofa future research. The results of Table 4b show that one should rely on at least N S = 10000 sample paths to get an acceptable price. Of course, this also dependson the granularity of the volume grid: we have chosen 100 equally-spaced steps. Thecomputational performance of the LSMC based on Algorithm 2 is inferior to that ofother numerical techniques such as FFT based approaches. Nevertheless, it has theadvantage to be applicable to any payoff function, in contrast to FFT techniques thathave to be adapted to each particular contract.13able 4: Gas Storage in a 1-Factors market dynamicsParameter Value F (0 ,
0) 15 T d κ . ν . σ . C (0) 0 C ( T ) 0 a in a w C max (a) Parameters N S CPU CPU ∗ price stdev error %-error1000 12 .
92 4 .
74 8 .
79 0 .
33 0 .
011 0 . .
84 47 .
08 5 .
20 0 .
49 0 .
005 0 . .
39 92 .
16 5 .
24 0 .
34 0 .
002 0 . .
98 139 .
01 5 .
10 0 .
37 0 .
002 0 . .
09 465 .
41 5 .
13 0 .
36 0 .
001 0 . (b) Storage Results In this paper we have introduced a three-steps procedure to determine the law of theincrement of generalized OU processes relying on the role of self-decomposability inthe theory of such processes. Based on this machinery and the results of Qu et al. [36],we have developed efficient and accurate algorithms for the exact simulation of theOU-VG and OU-SVG processes discussed in Cummins et al. [18, 19] and Cufaro etal. [15]. The algorithms are accurate, efficient, and have been numerically tested, withthe associated performance reported in detail. In addition, our three-steps proceduresimplifies some of the proofs presented in the cited papers and could be employed tofind simulation algorithms of other generalized OU processes that will the object offuture inquires.These results are instrumental to design algorithms to price derivative contractsin energy markets by MC simulation. To this end, we have considered the case of anAsian option in a market driven by the sum of a standard VG and a OU-VG processand the case of a fast-churn gas storage in a market driven by a OU-SVG processusing the LSMC method of Boogert and de Jong [11]. Although, MC methods areslower than other numerical solutions, our algorithms give the possibility to com-pute the entire price distribution of derivative contracts. Although the parametersestimation is not the focus of our study, MC based techniques can also be a viableroute to estimate the model parameters because the likelihood method is impracti-cable. Moreover, they can provide a graphical evidence that calibrated parameterscorrespond to realistic sample paths.It is also worth noticing that all the algorithms that we have discussed are basedon the sequential generation of processes. Therefore, a last topic deserving furtherinvestigation is the possibility to simulate OU-VG processes backward in time ex-tending the results of Pellegrino and Sabino [35] and Sabino [37].14 eferences [1] O. Bardou, S. Bouthemy, and G. Pag´es. Optimal quantization for the pricing ofswing options.
Applied Mathematical Finance , 16(2):183–217, 2009.[2] O. E. Barndorff-Nielsen, J. L. Jensen, and M. Sørensen. Some StationaryProcesses in Discrete and Continuous Time.
Advances in Applied Probability ,30(4):9891007, 1998.[3] O.E. Barndorff-Nielsen. Processes of Normal Inverse Gaussian Type.
Financeand Stochastics , 2(1):41–68, 1998.[4] O.E. Barndorff-Nielsen and N. Shephard. Non-Gaussian Ornstein-Uhlenbeck-based Models and some of their Uses in Financial Economics.
Journal of theRoyal Statistical Society: Series B , 63(2):167–241, 2001.[5] O.E. Barndorff-Nielsen and N. Shephard. Integrated OU Processes and non-Gaussian OU-based Stochastic Volatility Models.
Scandinavian Journal ofStatistics , 30(30):277–295, 2003.[6] H. Ben-Ameur, M. Breton, L. Karoui, and P. L’Ecuyer. A Dynamic Program-ming Approach for Pricing Options Embedded in Bonds.
Journal of EconomicDynamics and Control , 31(7):2212–2233, July 2007.[7] F.E. Benth, J. Kallsen, and T. Meyer-Brandis. A non-gaussian ornstein-uhlenbeck process for electricity spot price modeling and derivatives pricing.
Applied Mathematical Finance , 14(2):153–169, 2007.[8] F.E. Benth and A. Pircalabu. A non-Gaussian Ornstein-Uhlenbeck Model forPricing Wind Power Futures.
Applied Mathematical Finance , 25(1), 2018.[9] D. P. Bertsekas.
Dynamic Programming and Optimal Control, Volume I . AthenaScientific, Belmont, Mass., third edition, 2005.[10] M.L. Bianchi, S.T. Rachev, and F.J. Fabozzi. Tempered Stable Ornstein-Uhlenbeck Processes: A Practical View.
Communications in Statistics - Simu-lation and Computation , 46(1):423–445, 2017.[11] A. Boogert and C. de Jong. Gas Storage Valuation Using a Monte Carlo Method.
Journal of Derivatives , 15:81–91, 2008.[12] A. Cartea and M. Figueroa. Pricing in Electricity Markets: a Mean RevertingJump Diffusion Model with Seasonality.
Applied Mathematical Finance, No. 4,December 2005 , 12(4):313–335, 2005.[13] R. Cont and P. Tankov.
Financial Modelling with Jump Processes . Chapmanand Hall, 2004.[14] N. Cufaro Petroni. Self-decomposability and Self-similarity: a Concise Primer.
Physica A, Statistical Mechanics and its Applications , 387(7-9):1875–1894, 2008.1515] N. Cufaro Petroni, S. De Martino, S. De Siena, and F. Illuminati. StationaryDistribution of non-Gaussian Ornstein-Uhlenbeck Processes for Beam Halos. In
Proceedings of the Conference CYCLOTRONS 07 , page 424, 2008. GiardiniNaxos, Italy, 1-5 October 2007; INFN Catania.[16] N. Cufaro Petroni and P. Sabino. Fast Pricing of Energy Deriva-tives with Mean-reverting Jump-diffusion Processes. Available at:https://arxiv.org/abs/1908.03137.[17] N. Cufaro Petroni and P. Sabino. Gamma Related OrnsteinUhlenbeck Processesand their Simulation. available at: https://arxiv.org/abs/2003.08810.[18] M. Cummins, G. Kiely, and B. Murphy. Gas Storage Valuation under L´evyProcesses using Fast Fourier Transform.
Journal of Energy Markets , 4:43–86,2017.[19] M. Cummins, G. Kiely, and B. Murphy. Gas Storage Valuation under MultifactorL´evy Processes.
Journal of Banking and Finance , 95:167–184, 2018.[20] L. Devroye.
Non-Uniform Random Variate Generation . Springer-Verlag, 1986.[21] M. Gardini, P. Sabino, and E. Sasso. Correlating L´evy Processes withSelf-decomposability: Applications to Energy Markets. Available at:https://arxiv.org/abs/2004.04048.[22] D. P. Gaver and P. A. W. Lewis. First-order Autoregressive Gamma Sequencesand Point Processes.
Advances in Applied Probability , 12(3):727745, 1980.[23] I. S. Gradshteyn and I. M. Ryzhik.
Table of Integrals, Series, and Products .Elsevier/Academic Press, Amsterdam, seventh edition, 2007.[24] B. Grigelionis. On the Self-Decomposability of Euler’s Gamma Function.
Lithua-nian Mathematical Journal , 43(3):295–305, 2003.[25] B. Hambly, S. Howison, and T. Kluge. Information-Based Models for Financeand Insurance.
Quantitative Finance , 9(8):937–949, 2009.[26] S. Jaimungal and V. Surkov. L´evy Based Cross-commodity Models and Deriva-tive Valuation.
SIAM Journal of Financial Mathematics , 2:464–487, 2011.[27] M. Joshi.
The Concept and Practice of Mathematical Finance . Cambridge Uni-versity Press, 2003.[28] M. Kjaer. Pricing of Swing Options in a Mean Reverting Model with Jumps.
Applied Mathematical Finance , 15(5-6):479–502, 2008.[29] S. G. Kou. A Jump-Diffusion Model for Option Pricing.
Manage. Sci. ,48(8):1086–1101, August 2002.[30] A.J Lawrance. Some Autoregressive Models for Point Processes. In P. Bartfaiand J. Tomko, editors,
Point Processes and Queueing Problems (Colloquia Math-ematica Societatis J´anos Bolyai 24) , volume 24, pages 257–275. North Holland,1980. 1631] F. A. Longstaff and E.S. Schwartz. Valuing American Options by Simulation:a Simple Least-Squares Approach.
Review of Financial Studies , 14(1):113–147,2001.[32] J.J. Lucia and E.S. Schwartz. Electricity Prices and Power Derivatives: Evidencefrom the Nordic Power Exchange.
Review of Derivatives Research , 5(1):5–50,2002.[33] D. B. Madan and E. Seneta. The Variance Gamma (V.G.) Model for ShareMarket Returns.
The Journal of Business , 63(4):511–24, 1990.[34] T. Meyer-Brandis and P. Tankov. Multi-factor Jump-diffusion Models ofElectricity Prices.
International Journal of Theoretical and Applied Finance ,11(5):503–528, 2008.[35] T. Pellegrino and P. Sabino. Enhancing Least Squares Monte Carlo with Dif-fusion Bridges: an Application to Energy Facilities.
Quantitative Finance ,15(5):761–772, 2015.[36] Y. Qu, A. Dassios, and H. Zhao. Exact Simulation of Gamma-driven Ornstein-Uhlenbeck Processes with Finite and Infinite Activity Jumps.
Journal of theOperational Research Society , 0(0):1–14, 2019.[37] P. Sabino. Forward or Backward Simulations? A Comparative Study.
Quanti-tative Finance , 2020. Forthcoming.[38] K. Sato.
L´evy Processes and Infinitely Divisible Distributions . Cambridge U.P.,Cambridge, 1999.[39] W. Schoutens.
L´evy Processes in Finance: Pricing Financial Derivatives . JohnWiley and Sons Inc, 2003.[40] P. Schwartz and J.E. Smith. Short-term Variations and Long-term Dynamics inCommodity Prices.
Management Science , 46(7):893–911, 2000.[41] E. Taufer and N. Leonenko. Simulation of L´evy-driven OrnsteinUhlenbeck Pro-cesses with Given Marginal Distribution.