Correlating Lévy processes with Self-Decomposability: Applications to Energy Markets
aa r X i v : . [ q -f i n . P R ] A p r Correlating L´evy processes with Self-Decomposability:Applications to Energy Markets ∗ Matteo Gardini † Piergiacomo Sabino ‡ Emanuela Sasso § April 9, 2020
Abstract
The aim of this paper is to build dependent stochastic processes using the notionof self-decomposability in order to model dependence across different markets andextend some recently proposed multivariate L´evy models based on subordination.Consequently, we study the properties of such processes, derive closed form expressionsfor characteristic function and linear correlation coefficient and develop Monte Carloschemes for their simulation. These results are instrumental to calibrate the modelson power and gas energy European markets and to price spread options written ondifferent underlying assets using Monte Carlo and Fourier techniques.
Keywords : Multivariate L´evy Processes, Self-Decomposability, Monte Carlo, FFT,Energy Markets, Spread Options.
During the last decades a lot of efforts have been done to go beyond the Black and Scholes[4] framework in Financial Modelling. The Black-Scholes (BS) formula is widely used bypractitioners but its limits are well-known: real markets exhibit non normally distributedlog-returns, jumps in prices dynamic, volatility clustering, presence of extreme events andso on. Over the years a lot of researchers - Merton [30], Madan and Seneta [28] and Heston[22] among others - have proposed more sophisticated models to overcome these issues.Nevertheless, their focus is mainly on the single asset modelling framework.In the multivariate setting things get harder and we are mainly faced to three differentissues: • How can we extend a univariate model to a multivariate setting preserving mathe-matical tractability? ∗ The views, opinions, positions or strategies expressed in this work are those of the authors and donot represent the views, opinions and strategies of, and should not be attributed to E.ON SE. † Department of Mathematics, University of Genoa, Via Dodecaneso 16146, Genoa, Italy, email [email protected] ‡ Quantitative Modelling E.ON SE Br¨usseler Platz 1, 45131 Essen, Germany, email [email protected] § Department of Mathematics, University of Genoa, Via Dodecaneso 16146, Genoa, Italy, [email protected] How can we calibrate this model? • Which techniques can be used for derivatives pricing?The BS model is easily extensible to a multivariate framework (see Shreve [39]): indeed,in a Gaussian world, all the dependence structure is totally defined by covariance matrix.Heath et al. [21] proposed a multivariate model based on the Gaussian hypothesis that iswidely used in interest rate markets and also adapted to energy markets (see Barth andBenth [3]).Another approach to model dependence is using copulas, L´evy copulas and L´evy seriesrepresentation as illustrated in Cont and Tankov [12], Cherubini et al. [11], Panov andSamarin [31] and Panov and Sirotkin [32]. Unfortunately, these approaches, especiallyL´evy copulas, are hard to hand mathematically and are often hard to calibrate.In this study we address the three issues above in the context of multi-dimensionalprocesses, that are at least marginally L´evy, using multivariate subordination. To thisend, several approaches are available in the literature, for instance Barndorff-Nielsen et al.[2] or Luciano and Schoutens [26] use a common subordinator. In particular, in a series ofpapers Semeraro [38], Luciano and Semeraro [27], Ballotta and Bonfiglioli [1], Buchmannet al. [5] and Buchmann et al. [6] have proposed models based on subordination to intro-duce dependence between L´evy process. The common idea of these papers is to definemultivariate processes that are the sum of an independent process and a common process.For example Ballotta and Bonfiglioli [1] define a multivariate process in the following way: Y ( t ) = ( Y ( t ) , · · · , Y n ( t )) T = ( X ( t ) + a Z ( t ) , · · · , X n ( t ) + a n Z ( t )) T where Z ( t ), X j ( t ), j = 1 , · · · , n are independent L´evy processes. In a financial market,one can see the common process Z ( t ) as a systemic risk, whereas the independent pro-cesses X j ( t ) can be considered as an idiosyncratic component. The model has a simpleeconomical interpretation and it is mathematically tractable.The technique proposed by Semeraro [38] and Luciano and Semeraro [27] is very similarbut it is applied on subordinators. Semeraro [38] defines the process Y ( t ) as: Y ( t ) = ( Y ( t ) , · · · , Y n ( t )) = ( W ( G ( t )) , · · · , W n ( G n ( t )))where W j ( t ) j = 1 , · · · , n are independent BM and G j ( t ) j = 1 , . . . , n are subordinators.Each of these subordinators G j ( t ) is built as: G j ( t ) = X j ( t ) + α j Z ( t ) j = 1 , · · · , n where Z ( t ), X j ( t ), j = 1 , · · · , n , are independent subordinators and α j ∈ R + , j =1 , . . . , n .Since the correlation obtained with this model might not be as high as the one observedin the market (see Wallmeier and Diethelm [40]), Luciano and Semeraro [27] use the sametechnique proposed by Semeraro [38] but applied to correlated BM ’s.The purpose of this study is to extend these models using self-decomposable ( sd )subordinators: the idea is to replace the common process Z ( t ) by depended processes H ( t ) , . . . , H n ( t ). These processes are built using the notion of sd that provides an easyway to correlate random variables ( rv ). This construction allows us to introduce a new2arket dynamic similar to that observed in Cufaro Petroni and Sabino [15]: the processes( H ( t ) , . . . , H n ( t )) can be interpreted as news in the market occurring with random delays therefore, this can help us modeling co-movements in the markets.Indeed, once a market receives a news (or shock) after a certain random delay weobserve a propagation onto another market. We refer to this new type of processes as synaptic processes as proposed by Cufaro Petroni and Sabino [16].We show that our approach preserves mathematical tractability, indeed, we can derivethe characteristic functions, the linear correlation coefficient and design algorithms to sim-ulate these processes in the case of Gamma subordinators. It is also worthwhile observingthat our model goes beyond the mathematical generalization of the original ones providedby Buchmann et al. [5] and Buchmann et al. [6]: they analyze the case where the subor-dinator is sd . As it will be clear from the sequel, the a -reminder part of the subordinator H j ( t ) , j > sd .We then apply our market models to the pricing of spread options that are com-mon derivative contracts in energy markets. Since spread options are often not liquidenough, all our models are calibrated using a two-step methodology: firstly, one calibratesthe marginal parameters on quoted vanilla products and then, one fits the correlationstructure using historical quotations. Moreover, for these models some standard pricingtechniques based on the Fourier transform for multivariate setting (see for instance Hurdand Zhou [23], Caldana and Fusai [7] and Pellegrino [33]), or Monte Carlo (MC) methodsare available.The article is organized as follow: in Section 2 we give the basic notions that we needin the sequel. In Sections 3 and 5 we detail how to extend the models of Semeraro [38],Luciano and Semeraro [27] and Ballotta and Bonfiglioli [1] using sd subordinators, whereasin Sections 4 and 6 we propose MC and Fourier algorithms for simulation and pricing. InSection 7 we discuss about calibration and, eventually, in Section 8 we apply these modelsto Power and Gas Forward markets and we price spread options. Section 9 concludes thepaper. All the proof are given in Appendix A. In this section we introduce the fundamental concepts we need in the sequel: sd laws andBrownian subordination. We look at sd as a natural way to generate correlated rv and weuse this notion to build dependent stochastic processes in continuous time. Once we areable to generate increasing dependent stochastic processes, we can use the subordinationtechnique to define dependent subordinated Brownian motions ( BM ). We refer to Contand Tankov [12], Sato [37] and Cufaro Petroni [13] for the details. Definition 2.1.
We say that the law of a random variable X with pdf f ( x ) is sd if forevery a ∈ (0 , we can always find two independent random variables Y (with the samelaw of X ) and Z a with pdf g a ( x ) such that: X d = aY + Z a . If X has chf φ ( u ) and Z a (hereafter called Z a the a -reminder) has chf χ a ( u ) then wehave that: φ ( u ) = φ ( au ) χ a ( u ) . < a < rv ’s, Y (with the same law of X ) and Z a (here called a -remainder , with pdf g a ( x ) and chf χ a ( u )) such that X d = aY + Z a . (1)It is easy to see that a plays the role of correlation coefficient between X and Y . Wecan use this fact to build correlated stochastic processes in continuous time. To do so, weneed a very important property of sd laws that is given by Sato [37, Proposition 15.5]. Proposition 2.1.
If a law is sd then is infinitely divisible (id) and for a given a ∈ (0 , the law of Z a is uniquely determined and id. Other important concepts are the notions of subordinators, that are almost surelynon-decreasing L´evy processes, and Brownian subordination (see Cont and Tankov [12]).The idea behind subordination, is to replace the deterministic time t with a stochasticrunning time. One can use a non-decreasing L´evy process, called subordinator, G ( t ) to time-change a L´evy process obtaining a new one (Cont and Tankov [12, Theorem 4.2]). Ifthe time-change is done on a BM this operation is then called Brownian subordination. Definition 2.2.
Consider a probability space (Ω , F , P ) , µ ∈ R and σ ∈ R + . Let W ( t ) bea BM and let G ( t ) be a subordinator. A subordinated BM X ( t ) with drift is defined as: X ( t ) = µG ( t ) + σW ( G ( t )) (2)If H is P - a.s . non-negative random variables with sd law we can build sd subordinatorsas follows Definition 2.3 (Self-decomposable subordinators) . Let ˜ H and ˜ H be P - a.s . non-negativerv with sd laws and define H ( t ) and Z a ( t ) as the L´evy processes such that ( H i (1)) d = (cid:16) ˜ H i (cid:17) , i = 1 , and Z a (1) d = ˜ Z a . sd subordinators are defined as: H ( t ) = aH ( t ) + Z a ( t ) (3)Note that the process H ( t ) defined in (3) is a L´evy process because it is a linear com-bination of two L´evy processes (Cont and Tankov [12, Theorem 4.1]). The previous con-struction can be extended to the case n > sd subordinators can be used to buildsubordinated BM . Define H ( t ) = ( H ( t ) , · · · , H n ( t )) , n ∈ N setting: H ( t ) H ( t ) = a H ( t ) + Z a ( t ) · · · H n ( t ) = a n − H n − ( t ) + Z a n − ( t )where a j , j = 1 , . . . , n − ∈ (0 ,
1) and then defining the subordinated BM ’s: X j ( t ) = µ j H j ( t ) + σ j W j ( H j ( t )) , j = 1 , · · · , n where ( H j ) nj =1 are r.v. with sd law, W j ( t ) , j = 1 , . . . , n are mutually independent BM ’s.In next sections we extend some recently proposed multivariate L´evy models using the sd subordinator H ( t ). Hereafter, for the sake of notation simplicity, we consider only thebivariate case, i.e. n = 2. 4 Model extensions with Self-Decomposability
In this section we extend the models presented by Semeraro [38], Luciano and Semeraro[27] and Ballotta and Bonfiglioli [1] using sd subordinators. The key idea is that a marketis supposed to be driven by independent and dependent drivers. Dependence in our articleis not introduced using a common driver but using dependent renewal processes based onthe concept of sd introduced in Section 2. Suppose that the couple ( H ( t ) , H ( t )) isdefined as in (3), then these process are linked by the parameter a and by the distributionof Z a ( t ). If one interprets H ( t ) , H ( t ) as stochastic times, this model is appropriate tointroduce delays in the propagation of an information across different dependent markets.This new models are a generalization of the original ones that can be retrieved with asuitable parameterization. In this section we extend Semeraro’s model using sd subordinators. Definition 3.1 ( sd -Semeraro Model) . Let I j ( t ) j = 1 , be independent subordinators,and H ( t ) , H ( t ) be sd subordinators defined in (3) , independent of I j ( t ) . Define thesubordinator G j ( t ) G j ( t ) = I j ( t ) + α j H j ( t ) , j = 1 , with α j ∈ R + . Let µ j ∈ R , σ j ∈ R + , W j ( t ) be standard independent BM’s and let G j ( t ) subordinators as is (4) . Define the subordinated BM with drift Y j ( t ) as: Y j ( t ) = µ j G j ( t ) + σ j W j ( G j ( t )) , j = 1 , . (5)The chf of the process definite in (5) can be easily computed. Proposition 3.1 (Characteristic Function) . The joint chf φ Y ( t ) ( u ) of the process Y ( t ) =( Y ( t ) , Y ( t )) at time t defined in (5) is given by: φ Y ( t ) ( u ) = φ I ( t ) (cid:18) u µ + i σ u (cid:19) φ I ( t ) (cid:18) u µ + i σ u (cid:19) φ Z a ( t ) (cid:18) u µ + i σ u (cid:19) φ H ( t ) (cid:18) α (cid:18) u µ + i σ u (cid:19) + aα (cid:18) u µ + i σ u (cid:19)(cid:19) (6) Note.
We can observe that if we take the limit for a → a → φ Y ( t ) ( u ) = φ I ( t ) (cid:18) u µ + i σ u (cid:19) φ I ( t ) (cid:18) u µ + i σ u (cid:19) φ H ( t ) (cid:18) α (cid:18) u µ + i σ u (cid:19) + α (cid:18) u µ + i σ u (cid:19)(cid:19) and this is exactly the chf for the original model proposed by Semeraro [38]. Proposition 3.2 (Correlation) . The correlation at time t ρ Y ( t ) ,Y ( t ) is given by: ρ Y ( t ) Y ( t ) = µ µ α α aV ar [ H ( t )] p V ar [ Y ( t )] V ar [ Y ( t )] (7)5 .1.1 2D - Variance-Gamma In this subsection we extend Semeraro’s model for the Variance Gamma process defined inMadan and Seneta [28] using sd -subordinators. We recall that a Gamma rv has a density( pdf ) f ( α, β ; x ) and chf given by: f ( α, β ; x ) = β α Γ ( α ) x α − e − βx x> ( x ) ,φ X ( u ) = (cid:18) − iuβ (cid:19) − α (8)with α, β ∈ R + . It is well-known that if X ∼ Γ ( α, β ) then cX ∼ Γ (cid:16) α, βc (cid:17) and if X ∼ Γ ( α , β ) and Y ∼ Γ ( α , β ) are independent, then X + Y ∼ Γ ( α + α , β ). Now set in (5): I j ∼ Γ (cid:18) A j , Bα j (cid:19) , H j ∼ Γ (
A, B ) , j = 1 , α j H j ∼ Γ (cid:16) A, Bα j (cid:17) we have G j ∼ Γ (cid:18) A j + A, Bα j (cid:19) , j = 1 , . Remembering that A j , A, B, α j ∈ R + we have the following conditions:1 A j + A = α j B , j = 1 , < α j ≤ BA , j = 1 , E [ G j ] = 1 and then E [ G j ( t )] = t . Note.
If we request condition (9), we have that:1 = α ( A + A ) B = α ( A + A ) B and so the parameter B is somehow redundant and we can assume B = 1.The same conclusion can be reached observing that, in Equation (7), we fit only thevariance of H ( t ): for this reason assuming B = 1 is not restrictive.From Propositions 3.1 and 3.2 we have that: Corollary 3.3.
The chf in 2D Variance-Gamma case is: φ H j ( t ) ( u ) = (cid:16) − i uB (cid:17) − tA j = 1 , φ I j ( t ) ( u ) = (cid:16) − α j i uB (cid:17) − tA j j = 1 , φ Z a ( t ) ( u ) = φ H ( t ) ( u ) φ H ( t ) ( au ) = (cid:18) B − iuB − iau (cid:19) − tA (11) and so chf φ Y ( t ) ( u ) in (6) can be computed. Corollary 3.4.
Linear correlation coefficient in 2D Variance-Gamma case is given by: ρ ( Y ( t ) ,Y ( t )) = µ µ α α aA p σ + µ α p σ + µ α .2 Extension of Semeraro-Luciano’s Model Here we extend the model proposed by Luciano and Semeraro [27] using sd subordinators. Definition 3.2 ( sd -Luciano and Semeraro’s model) . Let I j ( t ) , j = 1 , , subordinatorsand let H ( t ) and H ( t ) two sd subordinators independent from I j ( t ) . Define the followingprocess: Y ρ ( t ) = µ I ( t ) + σ W ( I ( t )) + α µ H ( t ) + √ α σ W ρ ( H ( t )) µ I ( t ) + σ W ( I ( t )) + α µ H ( t ) + √ α σ (cid:16) W ρ ( aH ( t )) + ˜ W ( Z a ( t )) (cid:17) ! (12) where W ( t ) and W ( t ) are independent BM ’s whereas E [ dW ρ ( t ) dW ρ ( t )] = ρdt and ˜ W ( t ) is independent from W ( t ) = ( W ( t ) , W ( t )) and W ρ ( t ) = ( W ρ ( t ) , W ρ ( t )) . Proposition 3.5 (Characteristic Function) . The joint chf φ Y ρ ( t ) ( u ) of the process Y ρ ( t ) =( Y ρ ( t ) , Y ρ ( t )) at time t defined in (12) is given by: φ Y ( t ) ρ ( u ) = φ I ( t ) (cid:18) u µ + i σ u (cid:19) φ I ( t ) (cid:18) u µ + i σ u (cid:19) φ H ( t ) (cid:18) i u α σ (1 − a ) + u T µ + i u T a Σ u (cid:19) φ Z a ( t ) (cid:18) u µ α + i u α σ (cid:19) where µ = [ α µ , aα µ ] and Σ = (cid:20) α σ √ α α σ σ ρ √ α α σ σ ρ α σ (cid:21) Following the technique proposed for the proof of Proposition 3.2 one can show thefollowing:
Proposition 3.6 (Correlation) . The correlation at time t , ρ Y ρ ( t ) ,Y ρ ( t ) is given by: ρ Y ρ ( t ) ,Y ρ ( t ) = a (cid:0) µ µ α α V ar [ H ( t )] + ρσ σ √ α α E [ H ( t )] (cid:1)p V ar [ Y ( t )] V ar [ Y ( t )] (13) In order to build a 2D-Variance Gamma process we choose I j ∼ Γ (cid:18) A j , Bα j (cid:19) , H j ∼ Γ (
A, B ) , j = 1 , j + α j H j ∼ Γ (cid:18) A j + A, Bα j (cid:19) , j = 1 , E [ G j ] = 1 and, consequently, E [ G j ( t )] = t for j = 1 ,
2. Following the same argument of Section 3.1.1, the expression of linear cor-relation coefficient and the chf for the 2D Variance Gamma case can be derived.
Corollary 3.7.
Linear correlation coefficient in 2D Variance-Gamma case is given by: ρ ( Y ρ ( t ) ,Y ρ ( t ) ) = a (cid:0) µ µ α α A + ρAσ σ √ α α (cid:1)p σ + µ α p σ + µ α The chf can be obtained combining Corollary 3.3 with Proposition 3.5.
In this section we extend the model presented by Ballotta and Bonfiglioli [1] to the casewith sd subordinators. Definition 3.3 ( sd -Ballotta and Bonfiglioli’s model) . Let H ( t ) and H ( t ) be sd subor-dinators as in (3) and define: Y ( t ) = ( Y ( t ) , Y ( t )) = ( X ( t ) + a R ( t ) , X ( t ) + a R ( t )) (14) where: • X j ( t ) is a subordinated Brownian motion with parameters ( β j , γ j , ν j ) , j = 1 , ,where β j ∈ R is the drift, γ j ∈ R + is the diffusion and ν j ∈ R + is the variance ofthe subordinator. We state the two independent subordinators of X j ( t ) with G j ( t ) .We have that: X j ( t ) = β j G j ( t ) + W j ( G j ( t )) , j = 1 , . • R ( t ) and R ( t ) are given by: R ( t ) = β R H ( t ) + γ R W ( H ( t )) R ( t ) = β R H ( t ) + γ R (cid:16) W ( aH ( t )) + ˜ W ( Z a ( t )) (cid:17) (15) where W ( t ) and ˜ W ( t ) are independent Brownian motions and β R j ∈ R and γ R j ∈ R + . Lemma 3.8.
The chf of the process defined in (15) at time t is given by: φ R ( t ) ( u ) = φ H ( t ) (cid:18) u β R + u β R a + i (cid:0) u γ R + 2 u u γ R γ R a + u aγ R (cid:1)(cid:19) φ Z a ( t ) (cid:18) u β R + i u γ R (cid:19) Now we can compute the chf of the process defined in (14) as follows.8 roposition 3.9 (Characteristic Function) . The chf of the process at time t defined in (14) is given by: φ Y ( t ) ( u , u ) = φ G ( t ) (cid:18) β u + i u γ (cid:19) φ G ( t ) (cid:18) β u + i u γ (cid:19) φ R ( t ) ( a ◦ u ) (16) where a = ( a , a ) and u = ( u , u ) and ◦ is the Hadamard product.Note. Observe that:lim a → β R ,β R → β Z γ R ,γ R → γ Z φ Y ( t ) ( u , u ) = φ G ( t ) (cid:18) β u + i u γ (cid:19) φ G ( t ) (cid:18) β u + i u γ (cid:19) φ Z ( t ) (cid:18) β Z ( a u + a u ) + i a u + a u ) γ Z (cid:19) and this is the chf obtained by Ballotta and Bonfiglioli [1]. Proposition 3.10.
The correlation coefficient at time t of the process Y ( t ) defined in (14) is given by: ρ Y ( t ) = a a ( β R β R aV ar [ H ( t )] + γ R γ R a E [ H ( t )]) p V ar [ Y ( t )] p V ar [ Y ( t )] (17) It’s easy to show that, if X j ( t ) and R j ( t ) , j = 1 ,
2, are subordinated BM ’s with subor-dinators from the same family, then Y j ( t ) is a subordinated process of the same type of X j ( t ) and R j ( t ) if the following Ballotta and Bonfiglioli [1] style convolution conditionshold: ν R := ν R = ν R (18)and (cid:26) α j µ j = ν R a j β R j j = 1 , α j σ j = ν R a j γ R j j = 1 , H ( t ) and H ( t ) have the same law and so they have thesame variance ν R . It is easy to check that if Equations (19) are satisfied then: µ j = β j + a j β R j , σ j = γ j + a j γ R j , α j = ν j ν R / ( ν j + ν R ) . We can construct a 2D - Variance-Gamma using Gamma subordinators as follows. • Let H ( t ) ∼ Γ (cid:16) tν R , ν R (cid:17) be a Gamma subordinator and set H ( t ) = aH ( t ) + Z a ( t ).9 Let R j ( t ) be a subordinated BM (with drift β R j and diffusion γ R j ) obtained usingthe Gamma subordinator H j ( t ) ∼ Γ (cid:16) tν R , ν R (cid:17) , j = 1 , • Let X j ( t ) be a subordinated BM (with drift β j and diffusion γ j ) obtained using aGamma subordinator G j ( t ) ∼ Γ (cid:16) tν j , ν j (cid:17) , j = 1 , • Set Y j ( t ) = X j ( t ) + a j R j ( t )With this construction we have that Y ( t ) ∼ V G ( µ j , σ j , α j ) , j = 1 ,
2, where µ j , σ j , α j respect convolution conditions (19).The joint chf φ Y ( t ) ( u , u ) can be easly derived using (16) and remembering theexpression of the chf of a Gamma rv : φ ( u ) = (cid:18) − iuβ (cid:19) − α Corollary 3.11.
The correlation of the 2D - Variance-Gamma process Y ( t ) at time t isgiven by: ρ Y ( t ) = a a a ( β R β R ν R + γ R γ R ) p σ + µ α p σ + µ α sd -Gamma subordinators In this section we present algorithms to simulate the stochastic processes introduced inprevious sections. The choice of Gamma process was driven by the fact that we know howto efficiency simulate the a -remainder Z a ( t ). We could also adopt an Inverse Gaussianlaw, that is sd as shown in Halgreen [19]: however, an efficient way to simulate its a -remainder is not available at the moment. A numerical MC scheme to simulate a VarianceGamma process with parameters ( µ, σ, ν ) is easy to implement and can be found in Contand Tankov [12, Algorithm 6.11].In order to simulate sd dependent Gamma processes H ( t ) and H ( t ) defined in (3),one must know how to simulate a r.v. Z a ( t ). Cufaro Petroni and Sabino [17] and [14] haveshown that the a -remainder of a Gamma distribution Γ ( α, λ ) can be simulated taking: Z a = S X j =1 X j where S ∼ B ( α, − a ) X j ∼ E ( λ/a ) X = 0 P − a.s. B ( α, − a ) denote a Polya or negative binomial distribution and E ( λ/a ) denotes anexponential distribution. Algorithm 1 shows how to simulate Z a ( t ) over a fixed time grid t , t , . . . , t N , ∆ t n = t n − t n − , n = 1 , . . . , N .Different types of dependent sd subordinators varying a are shown in Figure 1.Algorithms 2, 3, 4 show how to simulate processes defined, respectively, in Section 3.1,Section 3.2 and Section 3.3. 10 lgorithm 1 for n = 1 , . . . , N do b ← B ∼ B ( α, − a ) ⊲ Generate a Polya ( α, − a ) rv ∆ Z an ← Z ( n ) a ∼ E b (cid:0) λ / a (cid:1) ; ⊲ Generate an Erlang rv with rate λ/ a Z ( t n ) = Z ( t n − ) + ∆ Z an . end forAlgorithm 2 sd -Semeraro modelSimulation of a bivariate Variance Gamma ( Y ( t ) , Y ( t )) in equally-spaced time grid withstep ∆ t , t , · · · , t i , · · · , t n . Simulate n iid gamma rv ’s ∆ I j,i ∼ Γ (cid:16) ∆ t (cid:16) α j − A (cid:17) , α j (cid:17) , j = 1 , i = 1 , . . . , n . Simulate n iid gamma rv ’s ∆ H ,i ∼ Γ (∆ t · A, i = 1 , . . . , n . Simulate n iid ∆ Z a,i , i = 1 , . . . , n rv ’s as a -remainder of a Gamma distributionΓ (∆ t · A,
1) using Algorithm 1. Set ∆ H ,i = a ∆ H ,i + ∆ Z a,i , i = 1 , . . . , n . Set ∆ G j,i = ∆ I j,i + α j ∆ H j,i , j = 1 , i = 1 , . . . , n . Simulate n iid standard normal rv ’s, N i i = 1 , . . . , n Set ∆ Y j,i = σ j N i p ∆ G j,i + µ j ∆ G j,i , j = 1 , i = 1 , . . . , n . The discretized trajectories are Y j ( t i ) = P ik =1 ∆ Y j,k , j = 1 , Algorithm 3 sd -Semeraro-Luciano modelSimulation of a bivariate Variance Gamma ( Y ( t ) , Y ( t )) in equally-spaced time grid withstep ∆ t , t , · · · , t i , · · · , t n . Simulate n iid Gamma rv ’s ∆ I j,i ∼ Γ (cid:16) ∆ t (cid:16) α j − A (cid:17) , α j (cid:17) , j = 1 , i = 1 , . . . , n . Simulate n independent gamma rv ’s ∆ H ,i ∼ Γ (∆ t · A, i = 1 , . . . , n . Simulate n independent ∆ Z a,i , i = 1 , . . . , n rv ’s as a -remainder of a Gamma distribu-tion Γ (∆ t · A,
1) using Algorithm 1. Set ∆ H ,i = a ∆ H ,i + ∆ Z a,i , i = 1 , . . . , n Simulate n iid standard normal rv ’s N j,i , j = 1 , i = 1 , . . . , n Set ∆ Y I j,i = σ N j,i p ∆ I j,i + µ j ∆ I j,i . Simulate n jointly standard normal rv (cid:16) ∆ W ρ , , ∆ W ρ , (cid:17) , i = 1 , . . . , n with correlation √ a . Simulate n standard normal iid rv ’s W , , · · · , W ,n , i = 1 , . . . , n . Set: ∆ Y H ,i = √ α σ ∆ W ρ ,i p ∆ H ,i + α µ ∆ H ,i , i = 1 , . . . , n . Set: ∆ Y H ,i = √ α σ (cid:16) ∆ W ρ ,i p a ∆ H ,i + p ∆ Z a,i ∆ W ,i (cid:17) + α µ ∆ H ,i , i = 1 , . . . , n . Set Y H j ( t i ) = P ik =1 ∆ Y H j ,k , j = 1 , i = 1 , . . . , n . Set Y I j ( t i ) = P ik =1 ∆ Y I j ,k , j = 1 , i = 1 , . . . , n . Set Y j ( t i ) = Y I j ( t i ) + Y H j ( t i ), j = 1 , i = 1 , . . . , n .11
50 100 150 200 250 300 350 4000246810 a = 0.01 H (t) H (t) a = 0.5 H (t) H (t) a = 0.7 H (t) H (t) a = 0.99 H (t) H (t) Figure 1: Correlated subordinators H ( t ) and H ( t ) with different values of a . Similar to Cont and Tankov [12], in this section we consider financial markets driven byexponential L´evy models using the processes Y ( t ) defined in the previous sections: S j ( t ) = S j (0) e rt + ω j t + Y j ( t ) (20) S j ( t ) is the value of the underlying at time t , r is the risk-free rate and ω j is the driftcorrector. By introducing a new source of randomness, the market is incomplete and we donot have a unique equivalent martingale measure (Shreve [39, Chapter 5]). Non-arbitrageconditions can be obtained setting: ω j = − ϕ ( − i ) (21)If we focus on the Variance Gamma process with parameters ( µ j , σ j , ν j ) we obtain: ω j = 1 ν j log − σ j ν j − µ j ν j ! . The payoff of spread options depends on the difference of two or more underlying assets.For an in-depth discussion about spread option pricing in log-normal setting one can referto Kirk [24] and Carmona and Durrleman [8]. Cheang and Chiarella [10] extended resultobtained by Margrabe [29] to the case of exchange options where both stock price processescontain compound Poisson jump components whereas Pellegrino and Sabino [34] applieda three-dimensional Fourier cosine series expansion method in order to price and hedgemultiasset spread options. In the more general L´evy framework, Hurd and Zhou [23]proposed a Fourier technique based on the double inversion of the Fourier transform toprice spread options. In contrast, Caldana and Fusai [7] and Pellegrino [33] developed an12 lgorithm 4 sd -Ballotta-Bonfiglioli modelSimulation of a bivariate Variance Gamma ( Y ( t ) , Y ( t )) in equally-spaced time grid withstep ∆ t , t , · · · , t i , · · · , t n . Simulate a Variance Gamma process X j ( t i ) with parameters ( β j , γ j , ν j ) j = 1 , Simulate n iid Gamma rv ’s ∆ H ,i ∼ Γ (cid:16) ∆ tν R , ν R (cid:17) , i = 1 , . . . , n . Simulate n iid rv ’s ∆ Z a,i as a -remainder of a Gamma distribution Γ (cid:16) ∆ tν R , ν R (cid:17) , i =1 , . . . , n . Set ∆ H ,i = a ∆ H ,i + ∆ Z a,i , i = 1 , · · · , n . Simulate n standard normal iid rv ’s N j,i , j = 1 , i = 1 , · · · , n . Set ∆ R ,i = γ R N ,i p ∆ S ,i + β R ∆ S ,i , i = 1 , . . . , n . Simulate n iid standard normal rv ’s N , , · · · , N ,n correlated with N ,i i = 1 , · · · , n with correlation √ a . Set ∆ R ,i = γ R (cid:0) N ,i p a ∆ S ,i + p ∆ Z a,i N ,i (cid:1) + β R ∆ S ,i , i = 1 , · · · , n . Set R j ( t i ) = P ik =1 ∆ R j,k , j = 1 , i = 1 , . . . , n . Set Y j ( t i ) = X j ( t i ) + R j ( t i ), j = 1 , i = 1 , . . . , n .approximated formula based on a single Fourier inversion. All these techniques require anexplicit expression of the chf Φ T ( u ) = E h e i h u T , Y ( T ) i i of the joint log-process.As shown above, we derived closed form expressions for the chf ’s of all the introducedmodels. For this reason we can use any of these Fourier methods: we chose anyway theapproximated method proposed by Caldana and Fusai [7] and we refer to their papers fordetails.Let S ( t ) and S ( t ) be two underlying assets, then the price of an European spreadoption with maturity T is given by: C K (0) = e − rT E (cid:2) ( S ( T ) − S ( T ) − K ) + (cid:3) Setting u = ( u , u ) T ∈ R and Y ( t ) = (log S ( t ) , log S ( t )) T they have the followingproposition. Proposition 6.1.
The approximate spread option value C k,αK (0) is given in term of aFourier inversion formula as: C k,αK (0) = (cid:18) e − δk − rT π Z + ∞ e − iγk Ψ T ( γ ; δ, α ) dγ (cid:19) + , (22) where δ is the dumping factor and Ψ T ( γ ; δ, α ) = e i ( γ − iδ ) log(Φ T (0 ,iα )) i ( γ − iδ ) [Φ T (( γ − iδ ) − i, − α ( γ − iδ )) − Φ T ( γ − iδ, − α ( γ − iδ ) − i ) − K Φ T ( γ − iδ, − α ( γ − iδ ))] and α = F (0 , T ) F (0 , T ) + K , k = log ( F (0 , T ) + K ) , F (0 , T ) = Φ T (0 , − i )13 Calibration
We use a two-steps calibration procedure as in Ballotta and Bonfiglioli [1] and Lucianoand Semeraro [27]: firstly, we calibrate the marginal parameters on vanilla products thatare usually liquid in many markets. Secondly, since derivatives written on more productsare very illiquid, the dependence parameters are estimated fitting the correlation matrixon historical series of the underlying assets.
Suppose we can observe n vanilla contract prices denoted with C i , i = 1 , . . . , n . Denotewith θ the vector of unknown marginal parameters and with C θ i ( K, T ) model prices. Theoptimal parameters set θ ∗ can be found solving: θ ∗ = arg min θ n X i =1 (cid:16) C θ i ( K, T ) − C i (cid:17) . (23)Since calibration procedure requires the computation of a large number of derivative pricesfor different values of θ , we need a very fast option pricing method. To this end, meth-ods based on FFT techniques, as proposed in Carr and Madan [9] and Lewis [25], arerecommended to compute C θ i ( K, T ) i = 1 , . . . , n in each iteration of the minimizationalgorithm. In our numerical experiment both methods produce very close results, and wechoose the method proposed by Carr and Madan [9]. Once the margins have been fitted, we fit the dependence structure looking at the historicalrealizations of log-returns of two underlying assets: we denote the market correlation with ρ M . Two different approaches can be used: the first one consists of simply observing thecorrelation of log-returns in the market ρ M and of trying to find η ∗ such that: η ∗ = arg min η (cid:0) ρ Y ( t, η ) − ρ M (cid:1) . (24)where η is the set of dependence parameters and ρ Y ( t,η ) is the correlation of the model:of course, some parameters constraints can be added.Another possible approach consists of using the generalized method of moments ( GMM )as proposed by Hansen [20]. Since the chf φ Y ( t ) ( u , u ) for the process Y ( t ) is availablein closed form, it is possible to compute the moment generating function M ( u , u ) andthen computing the mixed moments: E [ Y n Y m ] = ∂M n + m ( u , u ) ∂u n ∂u m (cid:12)(cid:12)(cid:12)(cid:12) u =0 u =0 (25)Then the GMM is applied and vector η ∗ estimated.Both methods have advantages and drawbacks. The problem in (24) is undeterminedbecause we have only one equation and several parameters to fit: for this reason the setof estimated parameters is not unique. In (25) numerical computation of higher momentslacks of precision if the data-set is not large enough.14e have chosen the first presented method because the size of the data-set is not verylarge and the higher moments computation might be inaccurate. We check the performance of models derived in Section 3 and apply them to energymarkets. This section is split into two parts: in the first one we apply our models to theGerman and French power forward markets, whereas in the second part our models areapplied to German power forward market and to natural gas forward market.We chose those markets because, in the first case we deal with markets that are stronglycorrelated due to the configuration of European electricity network, whereas, in the secondcase, the correlation between markets is still positive, since natural gas can be used toproduce electricity. The correlation in this latter case is not as strong as in the previousone. This gives us the opportunity to test our models for different level of correlations.For the sake of concision we use the following notation: • ( SSD ): sd -Semeraro’s model presented in Section 3.1. • ( LSSD ): sd -Luciano and Semeraro’s model presented in Section 3.2. • ( BBSD ): sd -Ballotta and Bonfiglioli’s model presented in Section 3.3.In our experiments we adopt the calibration procedure of Section 7 to fit the marketparameters to price spread options on future prices, now denoted F i ( t ) , i = 1 , S i ( t ) , i = 1 ,
2, whose payoff is given by:Φ T = ( F ( T ) − F ( T ) − K ) + . It customary to reserve the name
Cross-Border or Spark-Spread option if the futures arerelative to power or gas markets, respectively. In all our experiments we use the MCtechnique with N sim = 10 simulations and the Fourier-based method in Section 6. In order to calibrate our model we need both derivatives contracts written on forward andhistorical time series of forward quotations. For this reason, the data-set we relied uponis composed as follow: • Forward quotations from 25 April 2017 to 12 November 2018 of Calendar 2019 powerforward. A Forward Calendar 2019 contract is a contract between two counterpartsto buy or sell a specific volume of energy in MWh at fixed price for all the hours of2019. Calendar power forward in German and France are stated respectively withDEBY and F7BY. • Call Options on power forward 2019 quotations for both countries with settlementdate 12 November 2018. We used strikes in a range of ±
10 [
EU R/M W h ] aroundthe settlement price of the Forward contract, i.e. we exclude deep ITM and OTMoptions. We assume risk-free rate r = 0 . • The historical correlation between markets is ρ mkt = 0 . ǫ i defined as: ǫ i = C θi ( K, T ) − C i C i . We can observe this error is really small, varying K : our model is able to replicate marketprices and therefore can be used for pricing purposes.If we look at the fitted correlation the situation is slightly different. The SSD modelpresented in Section 3.1 fits a correlation that is roughly zero. For this reason the modelis not recommendable for
Cross-Border option pricing because it can overestimate thederivative price as we can see from the upper picture in Figure 2. The
LSSD model ofSection 3.2 is better than the previous one and the fitted correlation is very close to theone observed in the market as we can see from Table 3. For this reason the
LSSD modelcan be used to price
Cross-Border options. The
BBSD model derived by in Section 3.3provides an even better fitting of market correlation. We conclude that the
BBSD modelis the best one for the valuation of
Cross-Border options. A comparison between modelscan be found in the upper part of Figure 2: the option prices provided by the
BBSD modelare the lowest ones due to the highest value of fitted correlation.One additional consideration is needed: we note that, for all models, the fitted valuefor the sd parameters a is very close to one. This is because the linear correlation isproportional to a and, in order to maximize the correlation, a reasonable choice is tolet a →
1. As mentioned before, if a → Cross Border options, there’s not an essential difference between original models and theextended ones.
Model µ µ σ σ α α SSD
LSSD
BBSD
Table 1: Fitted marginal parameters for German and French power markets.
Parameter Value A B a ρ mod Table 2:
SSD
Parameter Value A B ρ a ρ mod Table 3:
LSSD
Parameter Value Parameter Value β -0.00 β R β γ R γ γ R γ ν R ν a ν ρ mod β R Table 4:
BBSD Strike [EUR/MWh] -1%-0.5%0%0.5% R e l a t i v e E rr o r s (SSD)-FR(LSSD)-FR(BBSD)-FR(SSD)-DE(LSSD)-DE(BBSD)-DE
45 50 55 60 65 -1%0%1%2%3%
France and German Call Pricing Relative Errors (FFT)
Strike [EUR/MWh] P r i c e [ E UR / M W h ] Cross-Border Options (MC) (SSD)(LSSD)(BBSD)
Figure 2: Percentage errors and Cross Border option prices.
In this section we present the numerical results of our models applied to the German powerforward market and to the Natural Gas forward market (TTF). These two markets arepositively correlated but not as strongly as power markets are.As in the power case, data-set we relied upon is the following one • Forward quotations from 1 July 2019 to 09 September 2019 relative to the MonthJanuary 2019 for the Power Forward in Germany and the Gas TTF Forward. • Call Options on power forward January 2020 quotations for both Germany and TTFwith settlement date 9 September 2019. As done before, we use strikes prices K ina range of ±
10 [
EU R/M W h ] around the settlement price of the forward contract,i.e. we exclude deep ITM and OTM options. • We assume risk-free rate r = 0 . • The historical correlation between log-returns is ρ mkt = 0 . ǫ is very small. In Figure 3 the pictureat the top shows that the SSD model overprices the
Spark-Spread option due to the factthat captured correlation is very low.
LSSD and
BBSD models provide a lower price ofthe derivatives because they are able to catch the market correlation. Fitted parametersare shown in Table 5: we observe that the sd parameter a is no more as close to one as itwas in the forward power markets. odel µ µ σ σ α α SSD
LSSD
BBSD
Table 5: Fitted marginal parameters for power and gas forward markets.
Parameter Value A B a ρ mod Table 6:
SSD
Parameter Value A B ρ a ρ mod Table 7:
LSSD
Parameter Value Parameter Value β β R β γ R γ γ R γ ν R ν a ν ρ mod β R Table 8:
BBSD
16 18 20 22 24 26 28
Strike [EUR/MWh] -2%-1.5%-1%-0.5%0%0.5% R e l a t i v e E rr o r s (SSD)-TTF(LSSD)-TTF(BBSD)-TTF(SSD)-DE(LSSD)-DE(BBSD)-DE
44 46 48 50 52 54 56 58 60 62 64 -0.2%0%0.2%0.4%0.6%0.8%
TTF and German Call Pricing Relative Errors (FFT)
25 30 35 40 45
Strike [EUR/MWh] P r i c e [ E UR / M W h ] Spark-Spread Options (MC) (SSD)(LSSD)(BBSD)
Figure 3: Percentage errors and Spark-Spread option prices.18
Conclusions
Based on the concept of self-decomposability, in this paper we have presented a newmethod to build dependent stochastic processes that are, at least, marginally L´evy. Wehave developed the theoretical setting and we have shown how sd subordinators can bebuilt starting from sd laws which are also infinitely divisible. Such processes are extremelyuseful if one wants to model dependence across markets. Applying this technique, we haveextended some recent works based on multivariate subordinators presented by Semeraro[38], Luciano and Semeraro [27] and Ballotta and Bonfiglioli [1] and we have derived ex-plicit expressions for the chf and the correlation. These results are instrumental to designMonte Carlo schemes and Fourier techniques employed to calibrate the models to real datain energy markets and to price Cross Border and Spark Spread options. We have con-sidered German and French power and gas forward markets using a two steps calibrationtechnique, consisting in fitting firstly marginal parameters on quoted vanilla products andsecondly, the correlation on historical realizations. Our numerical experiments have shownthat our proposed models can catch even extreme values of correlation between assets.Our approaches, and the relative developed numerical techniques, have been applied toenergy markets with two underlying assets only. Nevertheless, our modeling frameworkis very general and can be extended to all those situations where an innovation occurredin one process has an impact, with a certain random delay, to another processes. Such aframework can be used, for example, in equity derivatives, with an arbitrary number ofstocks, or in credit risk to model a chain of defaults caused by a common market shockthat propagates across markets.On the other hand, from a more mathematical perspective some points are still openand will be the objective of future inquires. For instance, our models have L´evy marginsbut it is still unclear whether the couple is still a L´evy process.In addition, although most of our results are general, we focused on sd Gamma sub-ordinators. It will be worthwhile investigating the case of Inverse Gaussian processes,and therefore Normal Inverse Gaussian processes, in more detail, deriving for instance, anefficient Monte Carlo algorithm to simulate the relative a -remainder where some intuitionmay come from the results in Dassios et al. [18]. Finally, a topic deserving further inves-tigation is the time-reversal simulation of such processes in order to efficiently price othercontracts like swings and storages via backward simulation as detailed in Pellegrino andSabino [35] and Sabino [36]. A Proofs
A.1 Proof of Proposition 3.1 (See page 5)
Proof.
Substituting the expression of Y j ( t ), conditioning with respect G j ( t ) and since W j ( t ) are independent we get: φ Y ( t ) ( u ) = E h e i h u , Y ( t ) i i = E h e iu Y ( t )+ iu Y ( t ) i = E " e i (cid:18) u µ + i σ u (cid:19) G ( t ) e i (cid:18) u µ + i σ u (cid:19) G ( t ) G j ( t ) we have: φ Y ( t ) ( u ) = E " e i (cid:18) u µ + i σ u (cid:19) I ( t ) e i (cid:18) u µ + i σ u (cid:19) I ( t ) e i (cid:18) u µ + i σ u (cid:19) α Z a ( t ) e i (cid:18)(cid:18) u µ + i σ u (cid:19) α + (cid:18) u µ + i σ u (cid:19) α a (cid:19) H ( t ) and, observing that I j ( t ), H ( t ) and Z a ( t ), are mutually independent the thesis follows. A.2 Proof of Proposition 3.2 (See page 5)
Proof.
We have to compute: cov ( Y ( t ) , Y ( t )) = E [ Y ( t ) Y ( t )] − E [ Y ( t )] E [ Y ( t )]Substituting the expressions of Y j ( t ) and G j ( t ) and observing that E [ H ( t ) H ( t )] = aV ar [ H ( t )]the thesis follows from straightforward computations. A.3 Proof of Proposition 3.5 (See page 7)
Proof.
Rewrite Y ρ ( t ) as: Y ρ ( t ) = Y I ( t ) + Y H ( t ) where: Y I ( t ) = (cid:18) µ I ( t ) + σ W ( I ( t )) µ I ( t ) + σ W ( I ( t )) (cid:19) and: Y H ( t ) = α µ H ( t ) + √ α σ W ρ ( H ( t )) α µ H ( t ) + √ α σ (cid:16) W ρ ( aH ( t )) + ˜ W ( Z a ( t )) (cid:17) ! The characteristic function is given by: φ Y ( t ) ρ ( u ) = E h e i h u , Y ρ ( t ) i i = E h e i h u , Y I ( t )+ Y H ( t ) i i = E h e i h u , Y I ( t ) i i E h e i h u , Y H ( t ) i i (26)We now compute the two last term separately. Substituting the expression of Y I , condi-tioning respect I j ( t ) , j = 1 , W ( t ) and W ( t ) are idependent wehave: 20 h e h u , Y I ( t ) i i = E h e i ( u µ + i u σ ) I ( t ) i E h e i ( u µ + i u σ ) I ( t ) i = φ I ( t ) (cid:18) u µ + i σ u (cid:19) φ I ( t ) (cid:18) u µ + i σ u (cid:19) (27)Following the same approach we can compute the second term, obtaining: E h e h u , Y H ( t ) i i = E h E h e iu α µ H ( t )+ iu √ α σ W ρ ( H ( t ))+ iu α µ aH ( t )+ iu √ α σ W ρ ( aH ( t )) | H ( t ) i E h e iu α µ Z a ( t )+ iu √ α σ ˜ W ( Z a ( t )) | Z a ( t ) ii Now we compute the inner expected values separately. The second inner expected valueis: E h e iu α µ Z a ( t )+ iu √ α σ ˜ W ( Z a ( t )) | Z a ( t ) i = e i ( u α + i u α σ ) Z a ( t ) For the second therm we have that, since H ( t ) is known: E h e iu α µ H ( t )+ iu √ α σ W ρ ( H ( t ))+ iu α µ aH ( t )+ iu √ α σ W ρ ( aH ( t )) | H ( t ) i = e iu α µ H ( t )+ iu α µ aH ( t ) E h e iu √ α σ W ρ ( H ( t ))+ iu √ α σ W ρ ( aH ( t )) | H ( t ) i The only unknown terms is the expected value. We have that: E h e iu √ α σ W ρ ( H ( t ))+ iu √ α σ W ρ ( aH ( t )) | H ( t ) i = e − u α σ (1 − a ) H ( t ) e − a u T a Σ u H ( t ) where Σ = (cid:20) α σ √ α α σ σ ρ √ α α σ σ ρ α σ (cid:21) and u = [ u , u ]. Setting µ = [ α µ , aα µ ] we can conclude that: E h e h u , Y H ( t ) i i = φ Z a ( t ) (cid:18) u α + i u α σ (cid:19) φ H ( t ) (cid:18) u T µ + i u α σ (1 − a ) + i a u T a Σ u (cid:19) (28)Using (27) and (28) in (26) we have the thesis. A.4 Proof of Lemma 3.8 (See page 8)
Proof.
Replacing the definition of R ( t ) and R ( t ) we get: φ R ( t ) ( u ) = E h e iu R ( t )+ iu R ( t ) i = E h e iu β R H ( t )+ iu aβ R H ( t )+ iu β R Z a ( t ) E h e iu γ R W ( H ( t ))+ iu γ R ( W ( aH ( t ))+ ˜ W ( Z a ( t )) ) | H ( t ) , Z a ( t ) ii
21e compute now the inner expected value: E h e iu γ R W ( H ( t ))+ iu γ R ( W ( aH ( t ))+ ˜ W ( Z a ( t )) ) | H ( t ) , Z a ( t ) i = E h e iu γ R W ( H ( t ))+ iu γ R W ( aH ( t )) | H ( t ) i E h e iu γ R ˜ W ( Z a ( t )) | Z a ( t ) i The second computation of the second expected value is immediate. E h e iu γ R ˜ W ( Z a ( t )) | Z a ( t ) i = e − u γ R Z a ( t ) For the first term we have: E h e iu γ R W ( H ( t ))+ iu γ R W ( aH ( t )) | H ( t ) i = e − (cid:16) u γ R +2 u u γ R γ R a + au γ R (cid:17) H ( t ) Observing that H ( t ) and Z a ( t ) are idependent the thesis follows. A.5 Proof of Proposition 3.9 (See page 9)
Proof.
Replacing the expression of Y j j = 1 , E h e h u , Y ( t ) i i = E h e iu X ( t ) i E h e iu X ( t ) i φ R ( t ) ( a ◦ u )Observe that, conditioning to G j ( t ), we have that: E h e iu j X j ( t ) i = E h e i ( u j β j + i u j γ j ) G j ( t ) i = φ G j ( t ) (cid:18) u j β j + i u j γ j (cid:19) This observation jointly with Lemma 3.8 complete the proof.
A.6 Proof of Proposition 3.10 (See page 9)
Proof.
Computing the covariance between Y ( t ) and Y ( t ) we have that: cov ( Y ( t ) , Y ( t )) = a a cov ( R ( t ) , R ( t )) (29)But, by direct computations, one can show that: cov ( R ( t ) , R ( t )) = β R β R aV ar [ H ( t )] + γ R γ R a E [ H ( t )] (30)where we used the following property: E [ W ( H ( t )) W ( aH ( t ))] = a E [ H ( t )]Using (29) and (30) we have the thesis. 22 eferences [1] L. Ballotta and E. Bonfiglioli. Multivariate Asset Models Using L´evy Processes andApplications. The European Journal of Finance , 13(22):1320–1350, 2013.[2] O.E. Barndorff-Nielsen, J. Pedersen, and K. Sato. Multivariate Subordination, Self-Decomposability and Stability.
Advances in Applied Probability , 33(1):160–187, 2001.[3] A. Barth and F.E. Benth. The Forward Dynamics in Energy Markets – Infinite-Dimensional Modelling and Simulation.
Stochastics , 86(6):932–966, 2014.[4] F. Black and M. Scholes. The Pricing of Options and Corporate Liabilities.
Journalof Political Economy , 81(3):637–654, 1973.[5] B. Buchmann, B. Kaehler, R. Maller, and A. Szimayer. Multivariate Subordina-tion Using Generalised Gamma Convolutions with Applications to Variance GammaProcesses and Option Pricing.
Stochastic Processes and their Applications , 127(7):2208–2242, 2017.[6] B. Buchmann, K. Lu, and D. Madan. Self-Decomposability of Variance GeneralisedGamma Convolutions. arXiv:1712.03640 [math.PR], 2019.[7] R. Caldana and G. Fusai. A General Closed-Form Spread Option Pricing Formula.
Journal of Banking & Finance , 12(37):4863–4906, 2016.[8] R. Carmona and V. Durrleman. Pricing and Hedging Spread Options.
SIAM Review ,45:627–685, 2003.[9] P. Carr and D.B. Madan. Option Valuation Using the Fast Fourier Transform.
Journal OF Computational Finance , 2:61–73, 1999.[10] G.H.L. Cheang and C. Chiarella. Exchange Options under Jump-Diffusion Dynamics.
Applied Mathematical Finance , 18(3):245–276, 2011.[11] U. Cherubini, E. Luciano, and Vecchiato V.
Copula Methods in Finance . WileyFinance, 2013.[12] R. Cont and P. Tankov.
Financial Modelling with Jump Processes . Chapman andHall, 2003.[13] N. Cufaro Petroni. Self-decomposability and Self-similarity: a Concise Primer.
Phys-ica A, Statistical Mechanics and its Applications , 387(7-9):1875–1894, 2008.[14] N. Cufaro Petroni and P. Sabino. Fast Pricing of Energy Derivatives with Mean-reverting Jump-diffusion Processes. available at: https://arxiv.org/abs/1908.03137.[15] N. Cufaro Petroni and P. Sabino. Coupling Poisson Processes by Self-decomposability.
Mediterranean Journal of Mathematics , 14(2):69, 2017.[16] N. Cufaro Petroni and P. Sabino. Pricing Exchange Options with Correlated JumpDiffusion Processes.
Quantitative Finance , 0(0):1–13, 2018.2317] N. Cufaro Petroni and P. Sabino. Gamma Related Ornstein–Uhlenbeck Processesand their Simulation. available at: https://arxiv.org/abs/2003.08810, 2020.[18] A. Dassios, Y. Qu, and H. Zhao. Exact Simulation for a Class of Tempered Stableand Related Distributions.
ACM Trans. Model. Comput. Simul. , 28(3), July 2018.[19] C. Halgreen. Self-decomposability of the Generalized Inverse Gaussian and HyperbolicDistributions.
Zeitschrift f¨ur Wahrscheinlichkeitstheorie und Verwandte Gebiete , 47(1):1432–2064, 1979.[20] L.P. Hansen. Large Sample Properties of Generalized Method of Moments Estimators.
Econometrica , 50:1029–1054, 1982.[21] D. Heath, R. Jarrow, and A. Morton. Bond Pricing and the Term Structure of InterestRates: A New Methodology for Contingent Claims Valuation.
Econometrica , 60(1):77–105, 1992.[22] S. L. Heston. A Closed-Form Solution for Options with Stochastic Volatility withApplications to Bond and Currency Options.
The Review of Financial Studies , 6(2):327–343, 1993.[23] T.R. Hurd and Z. Zhou. A Fourier Transform Method for Spread Option Pricing.https://arxiv.org/pdf/0902.3643.pdf, 2009.[24] E. Kirk. Correlations in the Energy Markets.
In Managing Energy Price Risk, RiskPublications and Enron , pages 71–78, 1995.[25] A. Lewis. A Simple Option Formula for General Jump-Diffusion and Other Expo-nential L´evy Processes. available at http://optioncity.net/pubs/ExpLevy.pdf, 2001.[26] E. Luciano and W. Schoutens. A Multivariate Jump-driven Finan-cial Asset Model.
Quantitative Finance , 6(5):385–402, 2006. URL https://doi.org/10.1080/14697680600806275 .[27] E. Luciano and P. Semeraro. Multivariate Time Changes for L´evy Asset Models:Characterization and Calibration.
Journal of Computational and Applied Mathemat-ics , 233(1):1937–1953, 2010.[28] D. B. Madan and E. Seneta. The Variance Gamma (V.G.) Model for Share MarketReturns.
The Journal of Business , 63(4):511–524, 1990.[29] W. Margrabe. The Value of an Option to Exchange One Asset for Another.
Journalof Finance , 33(1):177–186, 1978.[30] R.C. Merton. Options Pricing when Underlying Shocks are Discontinuous.
Journalof Financial Economics , 3:125–144, 1976.[31] V. Panov and E. Samarin. Multivariate Asset-Pricing Model Based on SubordinatedStable Processes.
Applied Stochastic Models in Business and Industry , 35(4):1060–1076, 2019. 2432] V. Panov and I. Sirotkin. Series Representations for Bivariate Time-Changed L´evyModels.
Methodology and Computing in Applied Probability , 19:97–119, 2017.[33] T. Pellegrino. A General Closed Form Approximation Pricing Formula for Basketand Multi-Asset Spread Options.
The Journal of Mathematical Finance , 6(5):944–974, 2016.[34] T. Pellegrino and P. Sabino. Pricing and Hedging Multiasset Spread Options Usinga Three-Dimensional Fourier Cosine Series Expansion Method.
Journal of EnergyMarkets , 7(2):71–92, 2014.[35] T. Pellegrino and P. Sabino. Enhancing Least Squares Monte Carlo with DiffusionBridges: an Application to Energy Facilities.
Quantitative Finance , 15(5):761–772,2015.[36] P. Sabino. Forward or Backward Simulations? A Comparative Study.
QuantitativeFinance , 2020. In press.[37] K. Sato.
L´evy Processes and Infinitely Divisible Distributions . Cambridge U.P., Cam-bridge, 1999.[38] P. Semeraro. A Multivariate Variance Gamma Model For Financial Applications.
International Journal of Theoretical and Applied Finance , 11(1):1–18, 2008.[39] S.E. Shreve.