On the control of the difference between two Brownian motions: an application to energy markets modeling

We derive a model based on the structure of dependence between a Brownian motion and its reflection according to a barrier. The structure of dependence presents two states of correlation: one of comonotonicity with a positive correlation and one of countermonotonicity with a negative correlation. This model of dependence between two Brownian motions $B^1$ and $B^2$ allows for the value of $\mathbb{P}\left(B^1_t - B^2_t \geq x\right)$ to be higher than $\frac{1}{2}$ when $x$ is close to 0, which is not the case when the dependence is modeled by a constant correlation. It can be used for risk management and option pricing in commodity energy markets. In particular, it allows to capture the asymmetry in the distribution of the difference between electricity prices and its combustible prices.

1. Introduction 1.1. Motivation. One of the major issues in commodity energy markets is the pricing and hedging of multi-assets options, in particular the spread options. For instance, if we denote by X t the price of electricity at time t, by Y t the price of coal at time t, and by H the heat rate (conversion factor) between the two, the income of the coal plant can be modeled by (X t − HY t − K) + , with x + = max (x, 0) and K representing a fixed cost. To evaluate the value of this coal plant, one needs to model jointly the price of electricity and the price of coal. Because the coal is a fuel for electricity, the two prices can not be considered independent and dependence between the two needs to be modeled. For more information on spread options, the reader can refer to [6].
To model energy commodities forward prices, two different approaches exist: the first one consists in the modeling of the spot price and is equivalent to the Vasicek modeling of interest rates [20], the second one consists in the modeling of the forward curve and is similar to the a Heath-Jarrow-Morton approach [12]. In the first approach, one way to model dependence is the use of structural models [1; 2; 5]. In structural models, electricity is a function of the residual demand and of the fuels used to produce it. Some constraints are imposed in order for the electricity price to be higher than the minimum price of its combustibles with a high probability. An other way to model dependence is the use of co-integration between the different commodities spot prices [17]. However, structural models and co-integration models are very computational costly and are not adapted for practitioners. They prefer to use the second type of models, using a forward curve.
We denote by f i (t, T ) the forward price of commodity i at time t with maturity T , that is of the delivery of commodity i during one unit of time. The most common model for f i (t, T ) is the two-factor model, see [3] for instance. The forward price of commodity i, i = 1, .., n is modeled by the following stochastic differential equation: (1) df i (t, T ) = f i (t, T ) σ i s e −α i (T −t) dB s,i t + σ i l dB l,i with B s,i and B l,i , i = 1, .., n, 2n brownian motions. The dependence between the Brownian motions is usually modeled by a constant correlation matrix. In the following, we are interested only in factorial models with two commodities, electricity and one of its fuel. Marginals model (if we consider only one commodity) are really efficient and allow us to price efficiently options based on one underlying. However, dependence modeling is not satisfying because it does not capture the asymmetry in the distribution of the difference between the forward price of electricity and the one of its fuel. Furthermore, the probability for the price of electricity to be lower than the price of its fuel is closed to 1 2 which is not consistent with the reality. Indeed, the fuel is used to produce the electricity.
Modeling the dependence between the forward prices is equivalent to the modeling of the dependence between the Brownian motions. We consider only two Brownian motions. To capture asymmetry, it is needed to consider an other approach than the constant correlation model. A common approach to construct a pair of Brownian motions is the use of stochastic correlation. Stochastic correlation models are a generalization in a multivariate framework of stochastic volatility models, such as the Heston model [13] where the volatility is modeled by a Cox-Ingersoll-Ross process. The matrix of volatility-correlation is stochastic and can be modeled for instance by a Wishart processes [11]. In a stochastic correlation framework, as the difference between the two Brownian motions does not follow a normal law, it is possible to capture asymmetry. However, if the stochastic correlation (ρ s ) s≥0 is independent from the two Brownian motions, we have for x ≥ 0: with Φ the normal cumulative distribution function. Stochastic correlation does not allow to have higher value then 1 2 for P B 1 t − B 2 t ≥ x . An other way to construct a pair of Brownian motions is the use of a local correlation. The concept of local correlation is directly derived from the one of local volatility. In a Black and Scholes framework, the volatility is constant with the maturity and strikes which is not coherent with the implied volatilities from call and put option prices. Dupire introduces the local volatility in order to have a price model which is compatible with the volatility smiles and which is a complete market model [9]. Langnau introduces local correlation model which is the generalization of local volatility for a multi-dimensional framework [16]. A less common approach is the use of copulae. Copulae are used to model the dependence between random variables and have many applications in finance [7]. Indeed, Sklar's theorem [19] states that modeling the distribution of a couple of random variables (X, Y ) is equivalent than modeling the law of X, the law of Y and a copula function C corresponding to the dependence between the two. However, use of copulae is more complicated in a continuous time framework, that is when processes are involved. In [4] and [14], a partial derivative equation is derived linking the copula between two Brownian motions and their local correlation function based on the Kolmogorov forward equation. Constraints on the copula to be admissible for Brownian motions are very restrictive, especially if one want to find asymmetric copulae admissible for Brownian motions. Deschatre [8] derives families of copula that are admissible for Brownian motions and asymmetric. Furthermore, he studies the range of the function C → P C B 1 t − B 2 t ≥ η for η > 0 and t > 0 where B 1 and B 2 are two Brownian motions and P C denotes the measure of probability when C is the copula of B 1 , B 2 . Some Markovian constraints are imposed on the copula C. The range of this function is equal to 0, 2Φ −η 2 √ t and the supremum is achieved with the copula between the Brownian motion and its reflection according to the barrier η 2 . However, those results are not adapted to a modeling framework because of the degenerescence of the model: the Brownian motions are either correlated to 1 or to -1 depending on the value of B 1 .

1.2.
Objectives and results. The main objective of this paper is to construct a model of dependence for solutions f 1 (t, T ) and f 2 (t, T ) of the stochastic differential equations (1). This model of dependence must create asymmetry in the difference between the two processes. In particular, we want to have high value for P f 1 (t, T ) − f 2 (t, T ) ≥ x with x close to 0. The dependence between the two processes is determined by the dependence between the Brownian motions. We reduce our problem to the construction of two Brownian motions B 1 = B 1 t t≥0 and B 2 = B 2 t t≥0 presenting asymmetry in their dependence and with values for P B 1 t − B 2 t ≥ x higher than 1 2 when x is close to 0. Our model is based on the work of Deschatre [8]. The value of P B 1 t − B 2 t ≥ η for η > 0 is maximized when B 2 is the reflection of B 1 according to the barrier η 2 . The copula between those two Brownian motions presents two states of dependence: one of comonotonicity corresponding to a correlation of 1 and one of countermonotonicity corresponding to a correlation of -1. We release these two states of dependence by allowing lower correlations in absolute value. This gives the copula of Proposition 1. This copula is asymmetric and Proposition 2 gives the survival function of difference between the two Brownian motions coupled with this copula: This model of dependence gives higher values for P B 1 t − B 2 t ≥ x than the constant correlation case and than 1 2 when x close to 0 and for ρ high enough. We generalize this model by allowing several reflections: it is the multi-barrier correlation model. We define two barriers ν and η with ν < η. We consider two independent Brownian motions X and B Y , and we construct the Brownian motion Y n that is correlated toX n : withX n the Brownian motion equal to −X at the beginning and reflecting when X −Y n hits a twostate barrier equal to η before the first reflection and switching from η to ν or from ν to η at each reflection. For a given x ∈ [η, ν] and t > 0, Corollary 1 states that the sequence P (X t − Y n t ≥ x) is increasing with n. Furthermore, the number of reflections in [0, t] N t is finite almost surely, see Proposition 3 (iii). We then consider the process Y t = Y Nt t which is a Brownian motion, see Proposition 3 (iv), that corresponds to the case n → ∞. Proposition 4 gives the survival function of X t − Y t , which is higher than in the constant correlation case and higher than 1 2 when x ∈ [η, ν] and ρ is high enough. This model can be transposed to a local correlation model: This system of stochastic differential has a strong solution (X, Y ), see Proposition 5. This model seems to be equivalent to the multi-barrier model when the two barriers have close values and ρ 2 = −ρ 1 = ρ. The solution has the advantage to be Markovian.
The multi-barrier correlation model is applied to the factorial model (1) in order to model jointly forward prices of electricity and forward prices of coal. Empirical results show that the model works well for products with a long delivery maturity (3 Month Ahead and 6 Month Ahead): the difference between the two products has an asymmetric distribution and the probability for the electricity product to be higher than the coal one is high. However, it is not the case for products with a short delivery maturity, such as the spot. This can be explained by a difference of volatility too high between the electricity spot price and the coal spot price. Indeed, the electricity and coal volatilities of the long term factors that drives the prices of long maturity products are close to each other whereas they are very different for the short term factors. An other limitations of our model is that it is highly sensitive to initial conditions, that is the initial prices of electricity and coal products. We also estimate prices of European spread options in our model with Monte Carlo. Results are the same in the local correlation model. 1.3. Structure of the paper. In Section 2, we provide a first model to construct two Brownian motions with a two-state correlation structure based on the dependence between a Brownian motion and its reflection. We give a closed formula for the survival function of the difference between the two Brownian motions: the distribution of the difference is asymmetric and can take higher values than in the constant correlation case. In Section 3, we improve the model of Section 2 by allowing several reflections to construct a multi-barrier correlation model. We give results about the survival function between the two Brownian motions and show that it takes higher values than the one in the model of Section 2. We also derive a local correlation model which gives the same results than the multi-barrier correlation model. Section 3 is our major contributions. In Section 4, we apply our results to the modeling of the forward prices of two commodities which are electricity and coal and to the pricing of spread options. Proofs are given in Section 5.

A two-state correlation copula
In this section, we derive a copula based on the Brownian motion and its reflection according to a barrier. As seen in [8], this copula contains two states depending on the value of the difference between the two Brownian motions: one of comonotonicity and one of countermonotonicity, that is correlation equal to 1 and -1. This copula maximizes P B 1 t − B 2 t ≥ η when the barrier is equal to η 2 , see [8,Proposition 3]. However, the dependence between the two Brownian motions when it is modeled by these copulae is degenerated in the sense that the difference between the two Brownian motions becomes constant in an infinite horizon. In this section, we construct a copula which does not present this degeneracy but which allows higher values for P B 1 t − B 2 t ≥ η than in the Gaussian copula case. The idea is to relax the correlation: instead of having states of correlation with correlations equals to 1 and -1, we have states of correlation with correlations equals to ρ and −ρ, |ρ| < 1.

2.1.
Model. Let us consider a filtered probability space (Ω, F, (F t ) t≥0 , P) with (F t ) t≥0 satisfying the usual hypothesis (right continuity and completion) and Thus,B k is a F Brownian motion according to the reflection principle (see [15, Theorem 3.1.1.2, p. 137]). Let ρ ∈ (0, 1) and Z a Brownian motion independent from B 1 . We consider the stochastic process . According to Skar's theorem [19], if X and Y are two random variables with continuous distribution function F X and F Y , there exists an unique copula C such that P (X ≤ x, Y ≤ y) = C F X (x) , F Y (y) . We will call C the copula of (X, Y ).
In the following, we will denote by Φ the cumulative distribution function of a standard normal random variable and by Φ ρ the cumulative distribution function of a bivariate gaussian vector of two standard normal random variables correlated with correlation ρ, ρ ∈ (−1, 1). Proposition 1 gives the copula between B 1 and B 2 .
is the copula between B 1 t and B 2 t at time t which are defined in the model of Section 2.1. This copula is clearly asymmetric in the sense that The copula contains two states of correlation: one of positive dependence (ρ > 0) and one of negative dependence (ρ < 0). Figure 1 gives the copula of Proposition 1 with ρ = 0.95 and in the degenerated case ρ = 1, h = 2 and t = 1.  2.3. Distribution of the difference between the two Brownian motions. Proposition 2 gives the survival function of B 1 t − B 2 t . Proposition 2. Let t > 0, h > 0, ρ ∈ (0, 1) and x ∈ R. Let B 1 and B 2 the two Brownian motions defined in the model of Section 2.1. We have:  This model allows us to have higher values than in the Gaussian copula case for P B 1 t − B 2 t ≥ z when z is close to 0. However, it presents some limitations in terms of modeling: Let us suppose that the barrier has already been crossed at time s, i.e.
Thus, the correlation between B 1 and B 2 at times t ≥ s is equal to ρ and does not change. We are in the same case than in the Gaussian copula case after time s, and then we do not optimize P

Multi-barrier correlation model
In this section, we improve the model of Section 2 by allowing several reflections. In the model of Section 2, once the reflection has happened, the two Brownian motions stay correlated with correlation ρ even if the difference between the two becomes low. We want to have two Brownian motions X and Y with the following correlation structure: if the value of X − Y is under a certain level that we denote by ν, X and Y have a negative correlation −ρ and if it is over an other level denoted by η, their correlation is positive and equal to ρ. One way to obtain this structure is to start with two Brownian motions having a negative correlation. When the difference between them reaches the barrier η, Y reflects and the correlation becomes positive. If the correlation is positive (resp. negative) and X − Y reaches ν (resp. η), Y reflects and the correlation becomes negative (resp. positive). The number of reflection that can happen is a parameter of our model denoted by n. Y is then correlated to a reflection of X reflecting each time the difference between the two reaches one of the two barriers. Figure 3 gives an illustration of our model. In Section 3.3, we develop a local correlation model based on the same principle. The local correlation model seems to be equivalent to the multi-barrier correlation model when the two barriers are close. Furthermore, in the local correlation model, the couple (X, Y ) is Markovian. 3.1. Model. Let B X and B Y be two independent Brownian motions defined on a common filtered probability space (Ω, F, (F t ) t≥0 , P) with (F t ) t≥0 satisfying the usual properties. We will denote indifferently B X by X.
where R(B, τ ) is the reflection Brownian motion of B with the reflection happening at time τ and τ a stopping time, i.e. R(B, Let N t = ∞ n=1 1 τn≤t be the number of reflections that happened before time t and Y t = Y Nt t . Y N is well defined because N t < ∞ almost surely according to Proposition 3 (iii). Proposition 3 gives results about the model.
k≥0 is a sequence of (F t ) t≥0 Brownian motions and (τ k ) k≥0 is a sequence of (F t ) t≥0 stopping times.
(iv) Y is a Brownian motion. Figure 4 is the empirical copula of (X t , Y n t ) for different n at time t = 1. The copula is asymmetric and we observe two states of correlation, as for the model of Section 2.

3.2.
Results on the distribution of the difference between the two Brownian motions. Proposition 4 gives an analytic formula for the survival function for X t − Y n t and X t − Y t .
Proposition 4. Let t > 0 and x ∈ R. Let (p n (t, x)) n≥0 the sequence defined by: where (u n ) n≥0 is the sequence defined by: and . is the floor function.
We have: and For x ∈ [ν, η], the survival function takes higher values than in the constant correlation case and than 1 2 . Furthermore, it is possible to increase the value of P (X t − Y n t ≥ x) by increasing the number of reflections with this model, which is why the case n = ∞ is considered.
Results of Proposition 1 are illustrated in Figure 5a. The case n = 0 corresponds to the Gaussian case. We can see that in [ν, η], the survival function is increasing with n. In Figure 5a, the curves for n = 5, n = 10 and n = 50 are the same. At time t = 1, the probability to cross more than 5 barrier is very weak then the Brownian reflection reflects less than 5 times with a high probability. The convergence in n at small time is fast. In Figure 5b, we can observe the difference between the cases n = 5, n = 10 and n = 50 at time t = 20. The survival function continues to grow. The survival function does not present the problem of the one the model of Section 2: its value stays high when t = 20 which is caused by the several reflections.
The results are confirmed with Figure 6. The higher the number of reflections is, more X − Y n is concentrated in the region [ν, η]. However, in the positive part of the plan, X − Y n take lower values than in the Gaussian case n = 0. One explanation comes from the martingality of and is higher than in the case of a constant correlation between the two Brownian motions. The probability mass in the positive part of the real line increases, but the expectation on all the real line stay the same: values taken by the random variables become lower in the positive part of the real line and becomes higher in the negative.We also remark that the symmetry present in the case n = 0 disappears when n is higher.  Let B X and B Y be two independent Brownian motions defined on a filtered probability space (Ω, F, (F t ) t≥0 , P).
Let us assume thatρ is Lipschitz.
Let us consider the following system of stochastic differential equations: Proposition 5 gives results about the solution of (5).
Proposition 5. The system of stochastic differential equations (5) has an unique strong solution (X, Y ) with X and Y two Brownian motions. Furthermore, (X, Y ) is Markovian.
Contrary to the multi-barrier correlation model, the local correlation model has the advantage to give a Markovian solution, which has some importance in practice. However, less analytical results are available for this model. In the following, we gives empirical results about it. Results are close to the ones of the multi-barrier correlation model.   As the local correlation function is asymmetric, i.e. ρ (x, y) = ρ (x, y) , x, y ∈ R, the copula of the solution of (5) is expected to be asymmetric. Figure 7 represents the copula of (X t , Y t ) at time t = 1. It is similar to the one of the multi-barrier correlation model. Figure 8 represents the survival function of the X t − Y t in the local correlation model at time t = 1 and t = 20 with parameters ν = 0, η = 0.5, ρ min = −0.9 and ρ max = 0.9. The local correlation function is chosen linear between ν and η. As for the multi-barrier correlation model, the distribution of X t − Y t is asymmetric. The survival function seems equivalent to the one of the multi-barrier correlation model. Between ν and η, the survival function is over 1 2 (Gaussian copula case). The survival function increases at the right of ν between time t = 1 and t = 20.

An application for joint modeling of commodity prices on energy market
In this section, we use the multi-barrier correlation model for the joint modeling of the forward prices of two commodities, electricity and coal. Coal is a fuel used to produce electricity which implies an asymmetry in the distribution of the difference between the two prices ; it is more likely that price of coal is lower than price of electricity (in the same unit). Modeling the dependence with a Gaussian copula is then not adapted. An advantage of our model is that it contains asymmetry in the distribution of the difference between the two prices. Furthermore, it allows not to change the marginal models.

4.1.
Model. Let us consider a two-factor model for both electricity and coal. For more information on the two-factor model, we refer to the study of Benth and Koekebakker [3].
Let f E (t, T ) (resp. f C (t, T )) the forward price of the electricity (resp. coal) at time t with maturity T , that is of the delivery of electricity (resp. coal) at maturity T during one day. Stochastic differential equation (6) gives dynamic of these products.
where B E,s , B E,l , B C,s , B C,l are standard Brownian motions defined on a common probability space (Ω, F, P).
In the dynamic of each commodity, there is one factor corresponding to the short term factor with a volatility σ i s e −α i s (T −t) , i = E, C . This short term factor is used to model the Samuelson effect [18], that is the decrease of volatility with time to maturity. The other factor is the long term factor with a constant volatility σ i l , i = E, C. Products traded on the market have a delivery period, except for the spot. We denote by f i (t, T, θ) , i = E, C the price of the product at time t that delivers i at time T during a period θ. By absence of arbitrage opportunities, we have In the following, we will only consider n Month Ahead (nMAH), n ≥ 1, which are products with a delivery period of one month and a delivery date which is the 1 st of the n th following month from today. Equation (7) gives the solutions of (6).
The spot price of electricity is given by S E t = f E (t, t) and the one of coal by S C t = f C (t, t). Then we have We model the dependence as follow: • B E,s and B E,l are independent, • B C,s and B C,l are independent, • B E,s and B C,s are independent, • B E,l and B C,l are constructed following the multi-barrier correlation model defined in Section 3. Usually, a constant correlation matrix is used to model the dependence between the 4 Brownian motions.

Parameters.
We consider the parameters of the marginal laws given in Table 1. Units are taken according to the year. We use the forward prices of electricity and of coal during 2014 in France to estimate these parameters. The method used for estimation is the first one of [10].
Parameters for the multi-barrier correlation model used to model the dependence between B E,l and B C,l are chosen arbitrarily ; we choose ν = 0, η = 0.5, ρ = 0.9, n = ∞.
In the benchmark model where dependence between B E,l and B C,l is modeled by a constant correlation, the correlation is equal to 0.275. The other correlation are equals to 0.
We assume that f E (0, T ) − Hf C (0, T ) = 0 and f E (0, T ) = 100 for all T (which does not represent the reality because we do not take into account the seasonality of the prices of electricity and coal). H is a conversion factor between the unit of electricity prices and the unit of coal prices and is called the heat rate.

4.3.
Numerical results. We are interested in the difference between f E (t, T ) and Hf C (t, T ). We only are interested in the multi-barrier correlation model ; results are the same for the local correlation model. Figure 10 represents the survival function of the difference between spot, 1MAH, 3MAH, and 6MAH prices. In the multi-barrier correlation model, the probability for the difference between the two spot prices to be non negative is close to 50%, which is the same value than in the benchmark model. However, we have good results if we consider long term products as 1MAH, 3MAH and 6MAH: we have probabilities closed to 60% for the 1MAH, and 70% for the 3MAH and  6MAH in the multi-barrier correlation model whereas we have probabilities closed to 50% in the benchmark model. The probability increases with the time to maturity. In the case of spot prices, the volatilities of the prices of the commodities is dominated by the short term factor, which we do not control ; in the other cases, these volatilities are small and the long term factor which we control dominates. This explains that we do not increase a lot the probability for the difference between the spot prices to be non negative. We also observed that in the multi-barrier correlation model, the survival function decreases faster than in the benchmark model and probability of being superior to 20 is closed to 0, which is not the case in the benchmark model. Figure 11 represents one trajectory of the different products. In the case of the spot prices, since electricity has a high volatility, it is difficult to control the difference between the two processes. For the other products, as the short term volatility decreases, we see that there is a control between the two processes. Remark 1. Using a multi-barrier correlation model to model the dependence between B E,s and B C,s does not improve the results for the different survival functions. That is why we consider them independent.
Results are sensitive to initial conditions. If we choose f E (0, T ) = 100 and Hf C (0, T ) = 120 for instance, f E (0, T ) − Hf C (0, T ) = −20 and we will have a distribution that is concentrated around -20, because the difference between the price is a martingale. The probability to be greater than -20 is higher in the multi-barrier correlation model than in the benchmark model but the probability to be positive is lower than in the benchmark model: it is closed to 0 in the multibarrier correlation model whereas it is closed to 10% in the benchmark model. Figure 12 represents the survival function of the difference between prices of electricity and coal for different products with ν = 0 and η = 0.5. As we choose a barrier near 0, the survival function will be maximized around -20.
One way to improve the value of the survival function around 0 is to choose a higher η. The idea in our model is that we want B E,l to go over B C,l + η, using correlation of -1 when the two  if we neglect the short term factors. We have . In the case with the same initial conditions, the right hand side term is equal to 0 and we choose a barrier of η. Heuristically, we then choose a barrier of ≈ 170.5 and ν = 170. Figure 13 gives the survival function of the different products in the multi-barrier correlation model with barriers ν = 170 and η = 170.5.
We can see that around 0, the values of the survival function are much better than in the benchmark model: around 20% in the multi-barrier correlation model and around 10% in the benchmark model. However, the values are still low. Indeed, even in the maximal case where the second Brownian motion is the reflection of the first one and the volatilities are equals, the probability for the difference between the Brownian motions to be positive knowing that one starts at −x, x > 0 and the other at 0 is equal to 2Φ −x 2 √ t which decreases with x.

4.4.
Pricing of European spread options. In this Section, we compare prices of European spread options in the factorial model (6) with different structures of dependence: the multi-barrier correlation model (m-b) with correlation equals to 0.3, 0.6, 0.9 and the benchmark model (constant correlation) with correlation equals to 0 and 0.275. Benchmark model with correlation equal to 0 is the same model than multi-barrier correlation model with correlation equal to 0. We price options with payoff (X t − HY t ) + , where X t is an electricity product, Y t a coal product and H is the conversion factor between electricity and coal. X t and Y t are Spot, 1MAH, 3MAH and 6MAH. Parameters used are those of Table 1. The price of the option is equal to E (X t − Y t ) + . We use Monte Carlo to estimate this expectation with a number of simulations equal to 10000. To simulate the processes, we use a step time of 1 hour. Table 2 gives 95% confidence intervals for the price of spread options with maturity 1 year when X 0 = HY 0 = 100. For the multi-barrier correlation model, we choose ν = 0 and η = 0.5. In the multi-barrier correlation model, the value of the option decreases with the correlation parameters. Indeed, when the correlation parameters increases, the probability to be over 0 is higher, but the values taken by the difference X t − HY t are smaller and smaller. The increase in the probability do not compensate the decrease in the values that can be taken and the expectation, i.e. the value of the option decreases. Value of the option in the benchmark model with correlation equal to 0.275 is close to the one in the multi-barrier correlation model with correlation equal to 0.6. We also observe that the value of the option decreases with the product maturity, in all the models. Table 3 gives 95% confidence intervals for the price of spread options with maturity 1 year when  Table 2. Values of European Spread options (X t − HY t ) + between electricity and coal products in the benchmark model and in the multi-barrier correlation model with parameters ν = 0, η = 0.5 with X 0 = HY 0 = 100. X 0 = 100 and HY 0 = 120. For the multi-barrier correlation model, we choose ν = 170 and η = 170.5. Contrarily to results of Table 2, the value of the option increases with the correlation parameter in the multi-barrier correlation model. Furthermore, the value of the option in the multi-barrier case is greater than the one of the benchmark model, for all the given correlations. In the constant correlation case, the probability to be greater than 0 is very low. The increase of probability in the multi-barrier correlation model is enough for the option value to be higher.  Table 3. Values of European Spread options (X t − HY t ) + between electricity and coal products in the benchmark model and in the multi-barrier correlation model with parameters ν = 170, η = 170.5 with X 0 = 100 and HY 0 = 120.

Proofs
5.1. Preliminary results. We start with well known results that will be useful for the proofs of propositions. Lemma 1. Let B = (B t ) t≥0 be a standard Brownian motion on a filtered probability space Ω, F, (F t ) t≥0 , P . We have: Proof The reader is referred to [15,  Lemma 2. Let B 1 = B 1 t t≥0 and B 2 = B 2 t t≥0 be two independent standard Brownian motion defined on a common filtered probability space Ω, F, (F t ) t≥0 , P with (F t ) t≥0 having all the good properties. Let h ≥ 0 and τ h = inf{t ≥ 0 : B 2 t = h}. We have: The same argument can be used to prove that Then we have We can conclude using Lemma 1.
Lemma 3. Let a, b and x ∈ R. We have: The reader is referred to [8, Proof of Lemma 19, Section 5.3].

5.2.
Proof of Proposition 2. Let B 1 and Z two independent Brownian motion. We consider B 2 = ρB h + 1 − ρ 2 Z withB h the reflection of B according to the barrier h. We have: Thus, P B 1 t − B 2 t ≥ x is the sum of the three following terms: Since B 1 and Z are independent, (i) is equal to the sum of the three following terms: with the use of Lemma 1. According to Lemma 3 (i), (9) is equal to Using Lemma 3 (ii), we find that the first term of (i) (12) is equal to In the same way, the second term of (i) (10) is equal to: The last one (11) is equal to The proof of (iii) can be done. {N t = n} = {τ n ≤ t, τ n+1 > t} and then we have According to Lemma 4, P (τ n ≤ t) = 2 Then nP (τ n ≤ t) = O n→∞ 1 n 2 and E (N t ) < ∞ by comparison theorem of positive series, implying N t < ∞ almost surely. (iv) If n ≥ N t , the number of reflections of X − Y n between time 0 and time t is equal to N t and Y Nt t = Y n t = almost surely.

5.4.
Proof of Proposition 4 and Corollary 1. We start with Lemma 5.
Lemma 5. For t > 0, x ∈ R, Proof We have: However, according to Equation (24), τ n ∼ τ = inf{t ≥ 0 : B t = u n } where B t is a standard Brownian motion. Then we have, using Lemma 2, The proof is the same for P (X t − Y n t ≤ x, t ≥ τ n+1 ).