Correlating L\'evy processes with Self-Decomposability: Applications to Energy Markets

The aim of this paper is to build dependent stochastic processes using the notion of self-decomposability in order to model dependence across different markets and extend some recently proposed multivariate L\'evy models based on subordination. Consequently, we study the properties of such processes, derive closed form expressions for characteristic function and linear correlation coefficient and develop Monte Carlo schemes for their simulation. These results are instrumental to calibrate the models on power and gas energy European markets and to price spread options written on different underlying assets using Monte Carlo and Fourier techniques.


Introduction
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 by practitioners but its limits are well-known: real markets exhibit non normally distributed log-returns, jumps in prices dynamic, volatility clustering, presence of extreme events and so 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 different issues: • How can we extend a univariate model to a multivariate setting preserving mathematical tractability?
• 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 is widely used in interest rate markets and also adapted to energy markets (see Barth and Benth [3]). Another approach to model dependence is using copulas, Lévy copulas and Lévy series representation as illustrated in Cont and Tankov [12], Cherubini et al. [11], Panov and Samarin [31] and Panov and Sirotkin [32]. Unfortunately, these approaches, especially Lévy 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-dimensional processes, that are at least marginally Lévy, using multivariate subordination. To this end, 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 of papers Semeraro [38], Luciano and Semeraro [27], Ballotta and Bonfiglioli [1], Buchmann et al. [5] and Buchmann et al. [6] have proposed models based on subordination to introduce dependence between Lévy process. The common idea of these papers is to define multivariate 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 1 (t) , · · · , Y n (t)) T = (X 1 (t) + a 1 Z (t) , · · · , X n (t) + a n Z (t)) T where Z (t), X j (t), j = 1, · · · , n are independent Lévy processes. In a financial market, one can see the common process Z (t) as a systemic risk, whereas the independent processes X j (t) can be considered as an idiosyncratic component. The model has a simple economical interpretation and it is mathematically tractable. The technique proposed by Semeraro [38] and Luciano and Semeraro [27] is very similar but it is applied on subordinators. Semeraro [38] defines the process Y (t) as: , · · · , Y n (t)) = (W 1 (G 1 (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 observed in the market (see Wallmeier and Diethelm [40]), Luciano and Semeraro [27] use the same technique 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 1 (t) , . . . , H n (t). These processes are built using the notion of sd that provides an easy way to correlate random variables (rv ). This construction allows us to introduce a new market dynamic similar to that observed in Cufaro Petroni and Sabino [15]: the processes (H 1 (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 we observe 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 derive the characteristic functions, the linear correlation coefficient and design algorithms to simulate these processes in the case of Gamma subordinators. It is also worthwhile observing that our model goes beyond the mathematical generalization of the original ones provided by Buchmann et al. [5] and Buchmann et al. [6]: they analyze the case where the subordinator is sd. As it will be clear from the sequel, the a-reminder part of the subordinator H j (t) , j > 1 process is infinitely divisible but not sd.
We then apply our market models to the pricing of spread options that are common derivative contracts in energy markets. Since spread options are often not liquid enough, all our models are calibrated using a two-step methodology: firstly, one calibrates the marginal parameters on quoted vanilla products and then, one fits the correlation structure using historical quotations. Moreover, for these models some standard pricing techniques based on the Fourier transform for multivariate setting (see for instance Hurd and Zhou [23], Caldana and Fusai [7] and Pellegrino [33]), or Monte Carlo (MC) methods are available.
The article is organized as follow: in Section 2 we give the basic notions that we need in 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, whereas in Sections 4 and 6 we propose MC and Fourier algorithms for simulation and pricing. In Section 7 we discuss about calibration and, eventually, in Section 8 we apply these models to Power and Gas Forward markets and we price spread options. Section 9 concludes the paper. All the proof are given in Appendix A.

Preliminaries
In this section we introduce the fundamental concepts we need in the sequel: sd laws and Brownian subordination. We look at sd as a natural way to generate correlated rv and we use this notion to build dependent stochastic processes in continuous time. Once we are able to generate increasing dependent stochastic processes, we can use the subordination technique to define dependent subordinated Brownian motions (BM ). We refer to Cont and 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 for every a ∈ (0, 1) we can always find two independent random variables Y (with the same law of X) and Z a with pdf g a (x) such that: If X has chf φ (u) and Z a (hereafter called Z a the a-reminder) has chf χ a (u) then we have that: Looking at the definition 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 (here called a-remainder, with pdf g a (x) and chf χ a (u)) such that It is easy to see that a plays the role of correlation coefficient between X and Y . We can use this fact to build correlated stochastic processes in continuous time. To do so, we need 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, 1) the law of Z a is uniquely determined and id.
Other important concepts are the notions of subordinators, that are almost surely non-decreasing Lévy processes, and Brownian subordination (see Cont and Tankov [12]). The idea behind subordination, is to replace the deterministic time t with a stochastic running time. One can use a non-decreasing Lévy process, called subordinator, G (t) to time-change a Lévy process obtaining a new one (Cont and Tankov [12,Theorem 4.2]). If the 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) be a BM and let G (t) be a subordinator. A subordinated BM X (t) with drift is defined as: If H is P-a.s. non-negative random variables with sd law we can build sd subordinators as follows Definition 2.3 (Self-decomposable subordinators). LetH 1 andH 2 be P-a.s. non-negative rv with sd laws and define H 1 (t) and Z a (t) as the Lévy processes such that (H i (1)) d = H i , i = 1, 2 and Z a (1) d =Z a . sd subordinators are defined as: Note that the process H 2 (t) defined in (3) is a Lévy process because it is a linear combination of two Lévy processes (Cont and Tankov [12,Theorem 4.1]). The previous construction can be extended to the case n > 2 and sd subordinators can be used to build subordinated BM. Define H (t) = (H 1 (t) , · · · , H n (t)) , n ∈ N setting: H n (t) = a n−1 H n−1 (t) + Z a n−1 (t) where a j , j = 1, . . . , n − 1 ∈ (0, 1) and then defining the subordinated BM 's: where (H j ) n j=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évy models using the sd subordinator H (t). Hereafter, for the sake of notation simplicity, we consider only the bivariate case, i.e. n = 2.

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 market is supposed to be driven by independent and dependent drivers. Dependence in our article is not introduced using a common driver but using dependent renewal processes based on the concept of sd introduced in Section 2. Suppose that the couple (H 1 (t) , H 2 (t)) is defined as in (3), then these process are linked by the parameter a and by the distribution of Z a (t). If one interprets H 1 (t) , H 2 (t) as stochastic times, this model is appropriate to introduce 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 a suitable parameterization.

Extension of Semeraro's Model
In this section we extend Semeraro's model using sd subordinators.
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: The chf of the process definite in (5) can be easily computed.
In this subsection we extend Semeraro's model for the Variance Gamma process defined in Madan and Seneta [28] using sd-subordinators. We recall that a Gamma rv has a density (pdf ) f (α, β; x) and chf given by: (5): Remembering that A j , A, B, α j ∈ R + we have the following conditions: Given the condition (9) and (10) we have that E [G j ] = 1 and then E [G j (t)] = t.
Note. If we request condition (9), we have that: 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 the variance of H 1 (t): for this reason assuming B = 1 is not restrictive. From Propositions 3.1 and 3.2 we have that: 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:

Extension of Semeraro-Luciano's Model
Here we extend the model proposed by Luciano and Semeraro [27] using sd subordinators.
Following the technique proposed for the proof of Proposition 3.2 one can show the following:

2D -Variance-Gamma
In order to build a 2D-Variance Gamma process we choose we have that: and, imposing conditions (9) and (10), we have get 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 correlation 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: The chf can be obtained combining Corollary 3.3 with Proposition 3.5.

Extension of Ballotta-Bonfiglioli's Model
In this section we extend the model presented by Ballotta and Bonfiglioli [1] to the case with sd subordinators.

Definition 3.3 (sd-Ballotta and Bonfiglioli's model).
Let H 1 (t) and H 2 (t) be sd subordinators as in (3) and define: where: • X j (t) is a subordinated Brownian motion with parameters (β j , γ j , ν j ) , j = 1, 2, where β j ∈ R is the drift, γ j ∈ R + is the diffusion and ν j ∈ R + is the variance of the subordinator. We state the two independent subordinators of X j (t) with G j (t).
We have that: • R 1 (t) and R 2 (t) are given by: where W (t) andW (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: Now we can compute the chf of the process defined in (14) as follows.
where a = (a 1 , a 2 ) and u = (u 1 , u 2 ) and • is the Hadamard product.
Note. Observe that: 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:

Convolution Conditions
It's easy to show that, if X j (t) and R j (t) , j = 1, 2, are subordinated BM 's with subordinators 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 conditions hold: and Relation (18) holds because H 1 (t) and H 2 (t) have the same law and so they have the same variance ν R . It is easy to check that if Equations (19) are satisfied then:

2D -Variance-Gamma
We can construct a 2D -Variance-Gamma using Gamma subordinators as follows.
• Let R j (t) be a subordinated BM (with drift β R j and diffusion γ R j ) obtained using • Let X j (t) be a subordinated BM (with drift β j and diffusion γ j ) obtained using a 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 1 , u 2 ) can be easly derived using (16) and remembering the expression of the chf of a Gamma rv : Corollary 3.11. The correlation of the 2D -Variance-Gamma process Y (t) at time t is given by:

Simulation Algorithms for sd -Gamma subordinators
In this section we present algorithms to simulate the stochastic processes introduced in previous sections. The choice of Gamma process was driven by the fact that we know how to efficiency simulate the a-remainder Z a (t). We could also adopt an Inverse Gaussian law, that is sd as shown in Halgreen [19]: however, an efficient way to simulate its aremainder is not available at the moment. A numerical MC scheme to simulate a Variance Gamma process with parameters (µ, σ, ν) is easy to implement and can be found in Cont and Tankov [12,Algorithm 6.11].
In order to simulate sd dependent Gamma processes H 1 (t) and H 2 (t) defined in (3), one must know how to simulate a r.v. Z a (t). Cufaro Petroni and Sabino [17] and [14] have shown that the a-remainder of a Gamma distribution Γ (α, λ) can be simulated taking: B (α, 1 − a) denote a Polya or negative binomial distribution and E (λ/a) denotes an exponential distribution. Algorithm 1 shows how to simulate Z a (t) over a fixed time grid t 0 , t 1 , . . . , t N , ∆t n = t n − t n−1 , 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. ∆Z a n ← Z (n) ⊲ Generate an Erlang rv with rate λ/ a 4: Z (t n ) = Z (t n−1 ) + ∆Z a n . 5: end for Algorithm 2 sd-Semeraro model Simulation of a bivariate Variance Gamma (Y 1 (t) , Y 2 (t)) in equally-spaced time grid with step ∆t, t 1 , · · · , t i , · · · , t n .

Financial Models
Similar to Cont and Tankov [12], in this section we consider financial markets driven by exponential Lévy models using the processes Y (t) defined in the previous sections: S j (t) is the value of the underlying at time t, r is the risk-free rate and ω j is the drift corrector. By introducing a new source of randomness, the market is incomplete and we do not have a unique equivalent martingale measure (Shreve [39,Chapter 5]). Non-arbitrage conditions can be obtained setting: If we focus on the Variance Gamma process with parameters (µ j , σ j , ν j ) we obtain:

Fourier Method for Spread Options
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 refer to Kirk [24] and Carmona and Durrleman [8]. Cheang and Chiarella [10] extended result obtained by Margrabe [29] to the case of exchange options where both stock price processes contain compound Poisson jump components whereas Pellegrino and Sabino [34] applied a three-dimensional Fourier cosine series expansion method in order to price and hedge multiasset spread options. In the more general Lévy framework, Hurd and Zhou [23] proposed a Fourier technique based on the double inversion of the Fourier transform to price spread options. In contrast, Caldana and Fusai [7] and Pellegrino [33] developed an Algorithm 4 sd-Ballotta-Bonfiglioli model Simulation of a bivariate Variance Gamma (Y 1 (t) , Y 2 (t)) in equally-spaced time grid with step ∆t, t 1 , · · · , t i , · · · , t n .
approximated formula based on a single Fourier inversion. All these techniques require an explicit expression of the chf Φ T (u) = E e i u T ,Y (T ) of the joint log-process.
As shown above, we derived closed form expressions for the chf 's of all the introduced models. For this reason we can use any of these Fourier methods: we chose anyway the approximated method proposed by Caldana and Fusai [7] and we refer to their papers for details.
Let S 1 (t) and S 2 (t) be two underlying assets, then the price of an European spread option with maturity T is given by: Setting u = (u 1 , u 2 ) T ∈ R 2 and Y (t) = (log S 1 (t) , log S 2 (t)) T they have the following proposition.
Proposition 6.1. The approximate spread option value C k,α K (0) is given in term of a Fourier inversion formula as: where δ is the dumping factor and

Calibration
We use a two-steps calibration procedure as in Ballotta and Bonfiglioli [1] and Luciano and Semeraro [27]: firstly, we calibrate the marginal parameters on vanilla products that are usually liquid in many markets. Secondly, since derivatives written on more products are very illiquid, the dependence parameters are estimated fitting the correlation matrix on historical series of the underlying assets.

Margins calibration
Suppose we can observe n vanilla contract prices denoted with C i , i = 1, . . . , n. Denote with θ the vector of unknown marginal parameters and with C θ i (K, T ) model prices. The optimal parameters set θ * can be found solving: Since calibration procedure requires the computation of a large number of derivative prices for different values of θ, we need a very fast option pricing method. To this end, methods based on FFT techniques, as proposed in Carr and Madan [9] and Lewis [25], are recommended to compute C θ i (K, T ) i = 1, . . . , n in each iteration of the minimization algorithm. In our numerical experiment both methods produce very close results, and we choose the method proposed by Carr and Madan [9].

Dependence structure calibration
Once the margins have been fitted, we fit the dependence structure looking at the historical realizations 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 the correlation of log-returns in the market ρ M and of trying to find η * such that: 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 1 , u 2 ) for the process Y (t) is available in closed form, it is possible to compute the moment generating function M (u 1 , u 2 ) and then computing the mixed moments: Then the GMM is applied and vector η * estimated. Both methods have advantages and drawbacks. The problem in (24) is undetermined because we have only one equation and several parameters to fit: for this reason the set of estimated parameters is not unique. In (25) numerical computation of higher moments lacks of precision if the data-set is not large enough.
We have chosen the first presented method because the size of the data-set is not very large and the higher moments computation might be inaccurate.

Numerical Results
We check the performance of models derived in Section 3 and apply them to energy markets. This section is split into two parts: in the first one we apply our models to the German and French power forward markets, whereas in the second part our models are applied 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 strongly correlated due to the configuration of European electricity network, whereas, in the second case, the correlation between markets is still positive, since natural gas can be used to produce electricity. The correlation in this latter case is not as strong as in the previous one. 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.
In our experiments we adopt the calibration procedure of Section 7 to fit the market parameters to price spread options on future prices, now denoted F i (t) , i = 1, 2 instead of S i (t) , i = 1, 2, whose payoff is given by: It customary to reserve the name Cross-Border or Spark-Spread option if the futures are relative to power or gas markets, respectively. In all our experiments we use the MC technique with N sim = 10 6 simulations and the Fourier-based method in Section 6.

Application to German and French Power Markets
In order to calibrate our model we need both derivatives contracts written on forward and historical time series of forward quotations. For this reason, the data-set 1 we relied upon is composed as follow: • We assume risk-free rate r = 0.015.
• The historical correlation between markets is ρ mkt = 0.94.
From Table 1 we see that the three models provide the same set of marginal parameters. In the lower box of Figure 2 we report the percentage error ǫ i defined as: We can observe this error is really small, varying K: our model is able to replicate market prices and therefore can be used for pricing purposes.
If we look at the fitted correlation the situation is slightly different. The SSD model presented in Section 3.1 fits a correlation that is roughly zero. For this reason the model is not recommendable for Cross-Border option pricing because it can overestimate the derivative price as we can see from the upper picture in Figure 2. The LSSD model of Section 3.2 is better than the previous one and the fitted correlation is very close to the one observed in the market as we can see from Table 3. For this reason the LSSD model can be used to price Cross-Border options. The BBSD model derived by in Section 3.3 provides an even better fitting of market correlation. We conclude that the BBSD model is the best one for the valuation of Cross-Border options. A comparison between models can be found in the upper part of Figure 2: the option prices provided by the BBSD model are the lowest ones due to the highest value of fitted correlation.
One additional consideration is needed: we note that, for all models, the fitted value for the sd parameters a is very close to one. This is because the linear correlation is proportional to a and, in order to maximize the correlation, a reasonable choice is to let a → 1. As mentioned before, if a → 1 we obtain the original models of Semeraro [38], Luciano and Semeraro [27] and Ballotta and Bonfiglioli [1]. For this reason, for Cross Border options, there's not an essential difference between original models and the extended ones. Model

Application to Power German and TTF Gas Future Market
In this section we present the numerical results of our models applied to the German power forward market and to the Natural Gas forward market (TTF). These two markets are positively correlated but not as strongly as power markets are. As in the power case, data-set 2 we relied upon is the following one • Call Options on power forward January 2020 quotations for both Germany and TTF with settlement date 9 September 2019. As done before, we use strikes prices K in a 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.015.
• The historical correlation between log-returns is ρ mkt = 0.54.
In the picture at the bottom of Figure 3 we can see that all models provide a good fitting of quoted market options because the error ǫ is very small. In Figure 3 the picture at the top shows that the SSD model overprices the Spark-Spread option due to the fact that captured correlation is very low. LSSD and BBSD models provide a lower price of the derivatives because they are able to catch the market correlation. Fitted parameters are shown in Table 5: we observe that the sd parameter a is no more as close to one as it was in the forward power markets. Model

Conclusions
Based on the concept of self-decomposability, in this paper we have presented a new method to build dependent stochastic processes that are, at least, marginally Lévy. We have developed the theoretical setting and we have shown how sd subordinators can be built starting from sd laws which are also infinitely divisible. Such processes are extremely useful if one wants to model dependence across markets. Applying this technique, we have extended some recent works based on multivariate subordinators presented by Semeraro [38], Luciano and Semeraro [27] and Ballotta and Bonfiglioli [1] and we have derived explicit expressions for the chf and the correlation. These results are instrumental to design Monte Carlo schemes and Fourier techniques employed to calibrate the models to real data in energy markets and to price Cross Border and Spark Spread options. We have considered German and French power and gas forward markets using a two steps calibration technique, consisting in fitting firstly marginal parameters on quoted vanilla products and secondly, the correlation on historical realizations. Our numerical experiments have shown that our proposed models can catch even extreme values of correlation between assets. Our approaches, and the relative developed numerical techniques, have been applied to energy markets with two underlying assets only. Nevertheless, our modeling framework is very general and can be extended to all those situations where an innovation occurred in one process has an impact, with a certain random delay, to another processes. Such a framework can be used, for example, in equity derivatives, with an arbitrary number of stocks, or in credit risk to model a chain of defaults caused by a common market shock that propagates across markets. On the other hand, from a more mathematical perspective some points are still open and will be the objective of future inquires. For instance, our models have Lévy margins but it is still unclear whether the couple is still a Lévy process.
In addition, although most of our results are general, we focused on sd Gamma subordinators. It will be worthwhile investigating the case of Inverse Gaussian processes, and therefore Normal Inverse Gaussian processes, in more detail, deriving for instance, an efficient Monte Carlo algorithm to simulate the relative a-remainder where some intuition may come from the results in Dassios et al. [18]. Finally, a topic deserving further investigation is the time-reversal simulation of such processes in order to efficiently price other contracts like swings and storages via backward simulation as detailed in Pellegrino and Sabino [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 e i u,Y (t) = E e iu 1 Y 1 (t)+iu 2 Y 2 (t) =E e i u 1 µ 1 +i and, observing that I j (t), H 1 (t) and Z a (t), are mutually independent the thesis follows.