L A Abbas-Turki

D Lamberton

European Options Sensitivity with Respect to the Correlation for Multidimensional Heston Models

This paper is devoted to the sensitivity study of the European option prices according to the correlation parameters when dealing with the multi-asset Heston model. When the Feller condition is not fulfilled, the CIR flow regularity is needed to prove the differentiability of the price according to the correlation. In the bidimensional case when the Feller condition is satisfied, the regularity of the volatility according to the correlation allows us to establish an asymptotic expression of the derivative of the price with respect to the correlation. This approximation provides the monotony for the exchange options then heuristically for spread option prices at short maturities. We also obtain this monotony for some restrictive choices of the products

where η i is the volatility of the volatility and ρ i is the asset/volatility correlation coefficient. Then, we explain how to extend the overall study to options written on more than two assets and on models that are derived from Heston model, like the double Heston model. We conclude by a large number of simulations that comfort the theoretical results.

Introduction

For a convex payoff, the authors of [START_REF] Karoui | Robustness of the Black & Scholes formula[END_REF] prove the monotony of the price of a European contract according to the volatility of the Black & Scholes (B&S) model. In the same fashion, let us prove the monotony according to the correlation parameter for the bidimensional B&S model

dS 1 t = S 1 t σ 1 dW 1 t , S 1 0 = x 1 0 , dS 2 t = S 2 t σ 2 (ρdW 1 t + √ 1 -ρ 2 dW 2 t ), S 2 0 = x 2 0 .
(

) 1
Let f be the convex payoff

f (s) = (a 1 s 1 + a 2 s 2 ± K) + = max(a 1 s 1 + a 2 s 2 ± K, 0), (a 1 , a 2 ) ∈ R 2 (2)
and F (t, x) the price of the studied contract, given under the risk-neutral probability by

F (t, x) = E ( f (S T ) S t = x ) = E t,x (f (S T )) .
Associated to the model [START_REF] Karoui | Robustness of the Black & Scholes formula[END_REF] and to the convex payoff [START_REF] Dimitroff | A parsimonious multi-asset Heston model: Calibration and derivatives pricing[END_REF], the price function

F (t, x) ∈ C 1,2 ( (0, T ) × R 2 + )
. This can be justified by the fact that the asset vector S T has a log-normal distribution which is sufficient to perform the wanted differentiations. Besides, F (t, x) satisfies the Black & Scholes PDE ∂F ∂t (t, x) + 1 2

( σ 2 1 ∂ 2 F ∂x 2 1 (t, x)x 1 1 + σ 2 2 ∂ 2 F ∂x 2 2 (t, x) ) x 2 2 + ρσ 1 σ 2 ∂ 2 F ∂x 1 ∂x 2 (t, x)x 1 x 2 = 0, F (T, x) = f (x),
We suppose now that the misspecified asset vector has the dynamic (1) but with a misspecified correlation ρ ̸ = ρ, that is to say

dS 1 t = S 1 t σ 1 dW 1 t , S 1 0 = x 1 0 , dS 2 t = S 2 t σ 2 (ρdW 1 t + √ 1 -ρ 2 dW 2 t ), S 2 0 = x 2 0 .
We have already seen that F (t, x) ∈ C 1,2 ( (0, T ) × R 2 +

) and using Ito calculus

F (T, S T ) = F (0, S 0 ) + 2 ∑ i=1 ∫ T 0 ∂F ∂x i (t, S t )dS i t + ∫ T 0 ∂F ∂t (t, S t )dt + 1 2 ∫ T 0 ( σ 2 1 ∂ 2 F ∂x 2 1 (t, S t ) [ S 1 t ] 2 + σ 2 2 ∂ 2 F ∂x 2 2 (t, S t ) [ S 2 t ] 2 + 2ρσ 1 σ 2 ∂ 2 F ∂x 1 ∂x 2 (t, S t )S 1 t S 2 t ) dt,
Combining the previous equality with the Black & Scholes PDE we get

E(F (T, S T )) = F (0, S 0 ) + (ρ -ρ)E {∫ T 0 σ 1 σ 2 ∂ 2 F ∂x 1 ∂x 2 (t, S t )S 1 t S 2 t dt } . ( 3
)
To compute the cross derivative, we consider the derivatives of S T with respect to S t = x (t < T ),

∂ x 1 S 1 T = S 1 T S 1 t , ∂ x 2 S 2 T = S 2 T S 2 t , ∂ x 2 S 1 T = ∂ x 1 S 2 T = 0.
When the payoff (2) is used then

∂ 2 F ∂x 1 ∂x 2 (t, x) = a 1 a 2 E t,x ( S 1 T x 1 S 2 T x 2 ε(a 1 S 1 T + a 2 S 2 T ± K) ) , ( 4
)
where ε is the Dirac distribution that can be justified thanks to the log-normal distribution g of

( S 1 T x 1 , S 2 T x 2
) :

∂ 2 F ∂x 1 ∂x 2 (t, x) = - a 2 x 1 ∫ R 2 + 1 a 1 x 1 v 1 +a 2 x 2 v 2 ≥±K ∂ v 1 [v 1 v 2 g(v 1 , v 2 )] dv 1 dv 2 = - a 1 x 2 ∫ R 2 + 1 a 1 x 1 v 1 +a 2 x 2 v 2 ≥±K ∂ v 2 [v 1 v 2 g(v 1 , v 2 )] dv 1 dv 2 .
From equality (4), a 1 a 2

∂ 2 F ∂x 1 ∂x 2 (t, x
) is clearly positive and the price is monotonous with respect to ρ. The direction of the latter monotony depends on the sign of the product a 1 a 2 . As an analogue of the implied volatility, thanks to the uniqueness of ρ one can define it as the implied correlation obtained from the market calibration of two assets that has the bidimensional B&S dynamics. As we will see in section 3, this notion of implied correlation is difficult to prove theoretically when using more complex models, like the Heston model.

In this paper, the assumed bidimensional version of the Heston model presumes the following dynamic for the couples asset/volatility (S i T , ν i T ) i=1,2

S 1 T = x 1 exp (∫ T t √ ν 1 s dW 1 s - 1 2
∫ T t ν 1 s ds ) , ( 5
)
S 2 T = x 2 exp (∫ T t √ ν 2 s ( ρdW 1 s + √ 1 -ρ 2 dW 2 s ) - 1 2
∫ T t ν 2 s ds ) , ( 6
)
ν 1 T = y 1 + κ 1 ∫ T t (θ 1 -ν 1 s )ds + η 1 ∫ T t √ ν 1 s dB 1 s , B 1 s = ρ 1 W 1 s + √ 1 -ρ 2 1 W 1 s , (7)
ν 2 T = y 2 + κ 2 ∫ T t (θ 2 -ν 2 s )ds + η 2 ∫ T t √ ν 2 s dB 2 s , B 2 s = ρ 2 ( ρW 1 s + √ 1 -ρ 2 W 2 s ) + √ 1 -ρ 2 2 W 2 s , (8)
where (W 1 , W 2 , W 1 , W 2 ) is a four-dimensional Brownian motion (these four Brownian motions are independent). We point out that the model specified by the previous SDEs does not include all the bidimensional Heston models. Indeed, the choice of this correlation structure is justified from a practitioner's point of view because it allows to calibrate simply each asset to the onedimensional put and call options, then add a correlation parameter ρ that can be calibrated from a spread option. Thus, the overall model will reproduce the prices of vanilla options and spread options. Although this model was already considered by various authors (see for example [START_REF] Dimitroff | A parsimonious multi-asset Heston model: Calibration and derivatives pricing[END_REF]) and widely used by practitioners, one of its drawbacks comes from constraining the correlation, between the Brownian motions of the two volatilities, to be equal to ρρ 1 ρ 2 .

Using the results of Bessel flow regularity in [START_REF] Vostrikova | On regularity properties of Bessel flow[END_REF], we study in section 2 the regularity of the CIR flow related to the SDEs [START_REF] Keller | Moment explosions and long-term behavior of affine stochastic volatility models[END_REF] and [START_REF] Del Bano Rollin | A new look at the Heston characteristic function[END_REF] then the volatility regularity with respect to the correlation of the Brownian motions. In section 3, we prove the differentiability of the price according to the correlation when the Feller condition is not fulfilled and we study some restrictive cases for which the price is monotonous with respect to the correlation. The derivative of ν 2 according to ρ is needed to establish in section 4 an asymptotic expression of the derivative of the price that works well for maturities T ≤ 0.3. In sections 3 and 4, we present also the basic ideas that allow to generalize our results to the multi-asset Heston and to models that are derived from Heston model, like the double Heston model. Thanks to a parallel implementation on the GPU Nvidia 480GTX, section 5 shows several tests of the error of our asymptotic approximation and it provides various Monte Carlo simulations that illustrate the monotony.

CIR flow regularity and volatility regularity with respect to the correlation

For a fixed t ≥ 0 and for s ≥ t, ν 1 and ν 2 share the same common CIR SDE given by

dν s = κ(θ -ν s )ds + η √ ν s ( rdW 1 s + √ 1 -r 2 dW 2 s ) , ν t = y, r ∈ [-1, 1], (9)
where here the Brownian motions W 1 and W 2 are independent but are not the same as the ones used in the previous section. However, it is quite clear that studying the flow of ν in ( 9) is equivalent to studying the flow of ν 1 and ν 2 in ( 7) and [START_REF] Del Bano Rollin | A new look at the Heston characteristic function[END_REF]. Moreover, the differentiability results of ν 2 with respect to ρ are similar to the differentiability results of ν with respect to r.

In this section, we use either the Feller condition

(A0) y > 0, 2κθ ≥ η 2 ,
or the following weaker assumption

(A1) y > 0, 4κθ > η 2 .
Introducing the process (0, ∞) ∋ y → τ 0 (y) defined by τ 0 (y) = inf {s ≥ t, ν s (y) = 0} , [START_REF] Alfonsi | High order discretization schemes for the CIR process: Application to affine term structure and heston models[END_REF] we refer for example to [START_REF] Lamberton | Introduction to Stochastic Calculus Applied to Finance[END_REF] for the proof of the finiteness of τ 0 (y) once (A0) is not satisfied, which means for a fixed y > 0 we have P (τ 0 (y) < ∞) = 1. The result of this section is summarized in the following theorem Theorem 2.1 Let ν be a CIR process driven by the SDE [START_REF]The Volatility Surface[END_REF]. Under the assumption (A0), both applications (0, ∞) ∋ y → ν s and (-

0 ≤ t ≤ s < τ 0 (y) then νs √ ν s = 1 √ ν t exp ( - κ(s -t) 2 - 1 2
( κθ - η 2 4 ) ∫ s t du ν u ) . ( 11
)
for these same assumptions and taking t = 0, ∂ r ν satisfies the following SDE

∂ r ν s = -κ ∫ s 0 ∂ r ν u du + η ∫ s 0 ∂ r ν u 2 √ ν u ( rdW 1 u + √ 1 -r 2 dW 2 u ) + η ∫ s 0 √ ν u ( dW 1 u - r √ 1 -r 2 dW 2 u ) , ∂ r ν 0 = 0, ( 12
)
that can be solved by a variation of constants method, to obtain

∂ r ν s = νs ( ηr ∫ s 0 √ ν u νu ( dW 1 u - r √ 1 -r 2 dW 2 u )) , ( 13
)
where in the latter equality, ν is the flow derivative at t = 0 (replace t by 0 in [START_REF] Glasserman | Monte Carlo Methods in Financial Engineering[END_REF]).

Note that (12) is only valid before time τ 0 (y). Therefore, in order to prove the differentiability of the price with respect to the correlation under (A1), we need additional work. This will be the main goal of section 3, in which we use the infinitesimal generator and the regularity of the flow. Unfortunately, the latter trick does not allow us to establish an asymptotic approximation and the only thing that we were able to do is to show that the asymptotic approximation, established when (A0) is fulfilled, works well numerically even for cases when only (A1) is satisfied.

Proof of Theorem 2.1:

We subdivide this proof into three steps

Step1: Proving the regularity of the flow. The solution of ( 9) is locally differentiable with respect to y, this means that we can differentiate with respect to y up to the time τ 0 (y) which is the upper limit of τ n 1/n (y) = inf{s ≥ t : ν n s ≤ 1/n}, such that ν n s is the solution of the truncated SDE associated to (9) with ν 1 t = y (we refer to [START_REF] Kunita | Stochastic flows and stochastic differential equations[END_REF] for more details). For s ∈ [t, τ 0 (y)[, we get νs = 1 -κ

∫ s t νu du + η ∫ s t νu 2 √ ν u ( rdW 1 u + √ 1 -r 2 dW 2 u ) . ( 14
)
By a change of variable using the logarithmic function, we obtain the solution of (14) for s < τ 0 (y) νs = exp

( -κ(s -t) + η ∫ s t rdW 1 u + √ 1 -r 2 dW 2 u 2 √ ν u - 1 2 η 2 ∫ s t du 4ν i u ) . ( 15
)
Moreover, by another change of variable X s = ln ( e κ(s-t) ν s ) , this time on the ν SDE, using Ito calculus we get exp ( η

∫ s t rdW 1 u + √ 1 -r 2 dW 2 u 2 √ ν u ) = √ ν u y exp ( κ(s -t) 2 - 1 2 ( κθ - η 2 2 ) ∫ s t du ν u ) ,
which combined with (15) provides [START_REF] Glasserman | Monte Carlo Methods in Financial Engineering[END_REF] for s < τ 0 (y).

According to [START_REF] Vostrikova | On regularity properties of Bessel flow[END_REF] (Proof of theorem 1.3), when δ ∈]1, 2[ the bessel flow (0, ∞) ∋ x → Θ(x, s), that satisfies (16) driven by the Brownian motion β

Θ(x, s) = x + β s + δ -1 2
∫ s t du Θ(x, u) , Θ(x, t) = x, ( 16
)
has a modification that admits a continuous derivative in probability sense that vanishes when s ≥ τ 0 (y). Consequently, one can use the same modification for the CIR flow because they are both related by the following equalities

ν s = exp [-κ(s -t)] Θ 2 ( √ y, η 2 4κ (e κ(s-t) -1) ) , δ = 4κθ η 2 . ( 17
)
To prove (17), we use Ito calculus on Z s = exp [κ(s -t)] ν s and we employ the time change

l s = 1 κ ln ( 4κs η 2
) + t to get

Z ls = y + 4κθ η 2 (s -t) + 2 ∫ s t √ Z lu dβ u .
Finally, defining Θ( √ y, s) = √ Z ls , we obtain (16) with x = √ y.

Step2: Proving the continuity of ν with respect to r. We define two Brownian motions

B s = rW 1 s + √ 1 -r 2 W 2 s and B s = rW 1 s + √ 1 -r 2 W 2
s thanks to which we set

dν s = κ(θ -ν s )ds + η √ ν s dB s , ν 0 = y, dν s = κ(θ -ν s )ds + η √ ν s dB s , ν 0 = y (18)
and we will prove that lim r→r ν s = ν s a.s. Let a n be a positive decreasing sequence defined by a n = a n-1 e -n , that satisfies

∫ a n-1 an dx x = n. (19) Afterwards, we set ϕ n ∈ C ∞ c (R) a mollifier function with support equal to [a n , a n-1 ] such that 0 ≤ ϕ n (x) ≤ ( * )
2 nx and (19) allows to have ∫ R ϕ n (x)dx = 1. Thanks to ϕ n , we define an approximation

ψ n (x) = ∫ |x| 0 dy ∫ y 0 ϕ n (z)dz (20) of the function absolute value | • |. Indeed |x| = ∫ |x| 0 dy ( ∫ y 0 ϕ n (z)dz + ∫ ∞ y ϕ n (z)dz ) = ψ n (x) + ∫ |x| 0 dy ∫ ∞ y ϕ n (z)dz thus, |x| ≥ ψ n (x) and because ∫ ∞ y ϕ n (z) ≤ 1 [0,a n-1 ] (y), then |x|-a n-1 ≤ ( * * ) ψ n (x).
In addition, for |x| ≥ a n (otherwise the first and the second derivative are equal to zero),

ψ ′ n (x) = Sgn(x) ∫ |x| 0 ϕ n (z)dz with |ψ ′ n (x)| ≤ ( * * * ) 1 [0,a n-1 ] (|x|) & ψ ′′ n (x) = ϕ n (|x|) ∈ [ 0, 2 n|x| ] .
Applying Ito calculus to ψ n (∆ s ), with ∆ s = ν s -ν s , we obtain

ψ n (∆ s ) = ∫ s 0 ψ ′ n (∆ u )d∆ u + 1 2 ∫ s 0 ψ ′′ n (∆ u )d ⟨∆⟩ u = -κ ∫ s 0 ψ ′ n (∆ u )∆ u du + L s + 1 2 ∫ s 0 ψ ′′ n (∆ u )d ⟨∆⟩ u ,
where

L s = ∫ s 0 ψ ′ n (∆ u )d(M u + N u
) is a square integrable martingale, because of the inequality ( * * * ) and that both M and N are two square integrable martingales (we refer the reader to [START_REF] Dufresne | The integrated square-root process[END_REF]) with

M s = η ∫ s 0 √ ν u d(B u -B u ), N s = η ∫ s 0 ( √ ν u - √ ν u )dB u . Employing Doob's L 2 -inequality on L * s = sup 0≤u≤s L s and ( * * * ) provide E ((L * s ) 2 ) ≤ 4E (L 2 s ) = 4E (∫ s 0 (ψ ′ n ) 2 (∆ u )d ⟨M + N ⟩ u ) ≤ 4E (∫ s 0 1 [0,a n-1 ] (|∆ u |)d ⟨M + N ⟩ u ) . ( 21
) Besides d ⟨∆, ∆⟩ u = d ⟨M + N, M + N ⟩ u ≤ 2d ⟨M, M ⟩ u + 2d ⟨N, N ⟩ u = 2η 2 ( ν u d ⟨ B -B ⟩ u + ( √ ν u - √ ν u ) 2 du ) . ( 22
)
Using both inequalities ( * * ) and ( * * * )

|∆ s | ≤ a n-1 + κ ∫ s 0 |∆ u |du + L * s + 1 2 ∫ s 0 ψ ′′ n (∆ u )d ⟨M + N, M + N ⟩ u .
Denoting the supremum X s = sup 0≤u≤s |∆ u | and using the inequality (a

+ b + c + d) 2 ≤ 4(a 2 + b 2 + c 2 + d 2 ) we get X 2 s ≤ 4a 2 n-1 + 4 ( κ ∫ s 0 X u du ) 2 + 4(L * s ) 2 + (∫ s 0 ψ ′′ n (∆ u )d ⟨M + N, M + N ⟩ u ) 2 ,
afterwards, we take the expectation and we use (21) and Cauchy-Schwarz inequality for the first integral term (s ≤ T )

E [X 2 s ] ≤ 4a 2 n-1 + 4T κ 2 E [∫ s 0 X 2 u du ] + 16E [∫ s 0 1 [0,a n-1 ] (|∆ u |)d ⟨M + N ⟩ u ] + E [ (∫ s 0 ψ ′′ n (∆ u )d ⟨M + N ⟩ u ) 2 ]
then ( 22), ( * ) and (

√ ν u - √ ν u ) 2 ≤ |ν u -ν u | provide E [X 2 s ] ≤ 4a 2 n-1 + 4κ 2 T E [∫ s 0 X 2 u du ] + 16η 4 E [∫ s 0 ν u d ⟨ B -B ⟩ u + ∫ s 0 1 [0,a n-1 ] (|∆ u |)|∆ u |du ] + 16η 4 n 2 E ( s + ∫ s 0 1 [an,a n-1 ] (|∆ u |) ν u |∆ u | d ⟨ B -B ⟩ u ) 2 ,
by continuing the computations, we obtain

E [X 2 s ] ≤ 4a 2 n-1 (1 + 4sη 4 ) + 4κ 2 T E [∫ s 0 X 2 u du ] + 16η 4 E [∫ s 0 ν u d ⟨ B -B ⟩ u ] + 16η 4 n 2 E ( s + ∫ s 0 ν u a n d ⟨ B -B ⟩ u ) 2 .
Let us take a sequence B = B k of Brownian motions that converges a.s. to B, once we apply Gronwall lemma

E [ X 2 s ] ≤ α(n, k)e 4κ 2 T 2 (23) with α(n, k) = 4a 2 n-1 (1 + 4sη 4 ) + 16η 4 s/n 2 + (16η 4 /n 2 )E (∫ s 0 ν u a n d ⟨ B k -B ⟩ u ) 2 + 16η 4 E (∫ s 0 ν u d ⟨ B k -B ⟩ u ) .
To conclude, we take n such that 4a

k such that 16η 4 n 2 E (∫ s 0 ν u a n d ⟨ B k -B ⟩ u ) 2 + 16η 4 E (∫ s 0 ν u d ⟨ B k -B ⟩ u )
< ϵ/2. The latter fact is possible because ν u admits moments of all orders (see the reference [START_REF] Dufresne | The integrated square-root process[END_REF]). Finally, we complete the proof of the continuity by using Fatou lemma on the left side of inequality (23).

Step3: Proving the differentiability of ν with respect to r.

Taking ϵ ∈]0, y[, we define τ ϵ (y), τ ϵ (y) and τ ϵ (y) as

τ ϵ (y) = inf{s > 0 : ν s (y) = ϵ}, τ ϵ (y) = inf{s > 0 : ν s (y) = ϵ}, τ ϵ (y) = τ ϵ (y) ∧ τ ϵ (y), (24)
We introduce the stochastic processes ∆ s = (ν s -ν s )/(r -r) and Λ s that satisfy the following SDEs

d∆ s = -κ∆ s ds + η √ ν s - √ ν s r -r dB s + η √ ν s d ( B s -B s r -r ) = -κ∆ s ds + η ∆ s √ ν s + √ ν s dB s + η √ ν s d ( B s -B s r -r ) , ∆ 0 = 0, dΛ s = -κΛ s ds + η Λ s √ ν s + √ ν s dB s , Λ 0 = 1,
which provides, by a variation of constants technique ∆ s∧ τϵ(y) = C s∧ τϵ(y) Λ s∧ τϵ(y) with

Λ s∧ τϵ(y) = exp ( -κ(s ∧ τ ϵ (y)) + η ∫ s∧ τϵ(y) 0 dB u √ ν u + √ ν u - η 2 2 ∫ s∧ τϵ(y) 0 du ( √ ν u + √ ν u ) 2 ) (25)
and

C s∧ τϵ(y) = ∫ s∧ τϵ(y) 0 η √ ν u Λ u d ( B u -B u r -r ) - ∫ s∧ τϵ(y) 0 η 2 √ ν u Λ u ( √ ν u + √ ν u ) d ⟨ B, B -B r -r ⟩ u ( 26
)
where

B u -B u r -r = W 1 u + √ 1 -r 2 - √ 1 -r 2 r -r W 2 u and ⟨ B, B -B r -r ⟩ u = ( r + √ 1 -r 2 √ 1 -r 2 - √ 1 -r 2 r -r ) u (27)
Now, we are going to take the limit as r → r in equality (25). For this task, we use the continuity result established in Step2, the lower bounds {ν u } u<s∧τ ϵ(y) ≥ ϵ, {ν u } u<s∧τϵ(y) ≥ ϵ and applying the dominated convergence theorem for the deterministic integral lim r→r Λ s 1 s< τϵ(y) = exp

( -κs + η ∫ s 0 dB u 2 √ ν u - η 2 2 ∫ s 0 du 4ν u )
1 s<τϵ(y) = νs 1 s<τϵ(y) .

In the limit above, the proof of the convergence of the stochastic term comes from the equality ∫ s∧ τϵ(y)

0 dB u √ ν u + √ ν u = r ∫ s∧ τϵ(y) 0 dW 1 u √ ν u + √ ν u + √ 1 -r 2 ∫ s∧ τϵ(y) 0 dW 2 u √ ν u + √ ν u
and from the Doob's L 1maximal inequality for the convergence of each term of this sum.

As for (26), let us first prove that the second term vanishes. Indeed, using (27), we have

lim r→r ⟨ B, B -B r -r ⟩ u = lim r→r ( r + √ 1 -r 2 √ 1 -r 2 - √ 1 -r 2 r -r ) u = 0. ( 28
) Thus lim r→r C s∧ τϵ(y) = lim r→r ∫ s∧ τϵ(y) 0 η √ ν u Λ u d ( B u -B u r -r ) = lim r→r [ ∫ s∧ τϵ(y) 0 η √ ν u Λ u dW ⊥ u + ( r √ 1 -r 2 + √ 1 -r 2 - √ 1 -r 2 r -r ) ∫ s∧ τϵ(y) 0 η √ ν u Λ u dW 2 u ] , with W ⊥ = ( W 1 - r √ 1 -r 2 W 2
)

. The limit of

( r √ 1-r 2 + √ 1-r 2 - √ 1-r 2 r-r ) is null, consequently lim r→r C s∧ τϵ(y) = lim r→r ∫ s∧ τϵ(y) 0 η √ ν u Λ u d ( W 1 u - r √ 1 -r 2 W 2 u )
and thanks to the independence of ν and

W ⊥ = W 1 - r √ 1 -r 2 W 2
(fact that can be seen directly form ( 14)), we have

lim r→r ∫ s∧ τϵ(y) 0 η √ ν u Λ u dW ⊥ u = 1 νs∧τϵ(y) lim r→r ∫ s∧ τϵ(y) 0 η νs∧ τϵ(y) √ ν u Λ u dW ⊥ u .
Once more, we employ Doob's inequality on M r T = sup 0≤s≤T ∫ s∧ τϵ(y)

0 νs∧ τϵ(y) √ ν u [ νu -Λ u ] νu Λ u dW ⊥ u and for λ > 0 P ( M r T ≥ λ ) ≤ 1 λ 2 sup 0≤s≤T E ( ∫ s∧ τϵ(y) 0 ν2 s∧ τϵ(y) ν u [ νu -Λ u ] 2 ν2 u Λ 2 u du ) ≤ 1 λ 2 E ( ∫ T 0 ν2 τϵ(y) ν u∧ τϵ(y) [ νu∧ τϵ(y) -Λ u∧ τϵ(y) ] 2 ν2 u∧ τϵ(y) Λ 2 u∧ τϵ(y) du ) ≤ 1 λ 2 ∫ T 0 E ( ν u∧ τϵ(y) [ νu∧ τϵ(y) Λ u∧ τϵ(y) -1 ] 2 )
du.

Thus, one can choose a sequence r k that tends to r such that:

∑ k≥1 P(M r k T ≥ λ) < ∞ and
Borel-Cantelli Lemma allows us to conclude for the a.s. convergence

lim r→r C s∧ τϵ(y) = η ∫ s∧τϵ(y) 0 √ ν u νu ( dW 1 u - r √ 1 -r 2 dW 2 u
) .

Finally, we have for s < τ 0 (y)

∂ r ν s = lim ϵ→0 lim r→r ∆ s = η νs ∫ s 0 √ ν u νu ( dW 1 u - r √ 1 -r 2 dW 2 u
) .

Sensitivity using the infinitesimal generator

The presentation of this part is subdivided into two subparts: In section 3.1, we reuse the same operations performed in the introduction (section 1) but with stochastic volatility models. We also present the result of the formal computations to show the key tools that allow to extend the proven results obtained in section 3.2 for the bidimensional Heston model. Thus the last part of section 3.1 can be skipped for a first reading.

A general framework for stochastic volatility models

In this part, we suppose that the real price of the asset vector given by the market has the following multidimensional stochastic volatility dynamic

dS i t = S i t √ ν i t dZ i t , S i 0 = x i 0 , dν i t = b i (t, ν i t )dt + σ i (t, ν i t )d Z i t , ν i 0 = y i 0 , (29)
where (Z 1 t , ..., Z d t , Z 1 t , ..., Z d t ) is a vector of correlated Brownian motions. Let f be a payoff of a multidimensional European contract on the considered asset vector, the price F (t, x, y) of this contract is given by

F (t, x, y) = E ( f (S T ) S t = x, ν t = y ) = E t,x,y (f (S T )) . Thus F (t, x, y) satisfies the Black & Scholes PDE ∂F ∂t (t, x, y) + d ∑ i ∫ T 0 ∂F ∂y i (t, x, y)b i (t, y) + 1 2 2d ∑ i,j=1 ∂ 2 F ∂z i ∂z j (t, x, y)Γ ij (t, x, y) = 0, F (T, x, y) = f (x), with z i = x i if i ≤ d and z i = y i-d if i > d and Γ(t, x, y) has the following expression Γ(t, x, y) = ( ( ρ ′ ij x i x j √ y i y j ) 1≤i,j≤d ( ρ ij x i √ y i σ j (t, y j ) ) 1≤i,j≤d ( ρ ij x i √ y i σ j (t, y j ) ) 1≤i,j≤d ( ρ ′′ ij σ i (t, y i )σ j (t, y j ) ) 1≤i,j≤d ) , with R is the 2d × 2d correlation matrix of the vector (Z 1 t , ..., Z d t , Z 1 t , ..., Z d t ) of standard Brownian motions R = ( ρ ′ ρ ρ ρ ′′ ) . ( 30
)
We suppose now that the misspecified price of the asset vector has the dynamic (29) but with R ̸ = R and different volatility of the volatility parameters σ i ̸ = σ i , that is to say

dS i t = S i t √ ν i t dZ i t , S i 0 = x i 0 , dν i t = b i (t, ν i t )dt + σ i (t, ν i t )d Z i t , ν i 0 = y i 0 . ( 31
)
Using formally Ito calculus

F (T, S T , ν T ) = F (0, S 0 , ν 0 ) + d ∑ i=1 ∫ T 0 ∂F ∂x i (t, S t , ν t )dS i t + ∫ T 0 ∂F ∂t (t, S t , ν t )dt + d ∑ i=1 ∫ T 0 ∂F ∂y i (t, S t , ν t )dν i t + 1 2 2d ∑ i,j=1 ∫ T 0 ∂ 2 F ∂z i ∂z j (t, S t , ν t )Γ ij (t, S t , ν t )dt,
where

z i = x i if i ≤ d and z i = y i-d if i > d and the matrix Γ(t, x, y) has the following expression Γ(t, x, y) = ( ( ρ ′ ij x i x j √ y i y j ) 1≤i,j≤d ( ρ ij x i √ y i σ j (t, y j ) ) 1≤i,j≤d ( ρ ij x i √ y i σ j (t, y j ) ) 1≤i,j≤d ( ρ ′′ ij σ i (t, y i )σ j (t, y j ) ) 1≤i,j≤d ) , ( 32
)
and R is the 2d × 2d correlation matrix of the vector

(Z 1 t , ..., Z d t , Z 1 t , ..., Z d t ) of standard Brow- nian motions R = ( ρ ′ ρ ρ ρ ′′ ) . ( 33
)
Taking the expectation of the previous equality and using the localization for the local martingale term

E(F (T, S T , ν T )) = F (0, S 0 , ν 0 ) + 1 2 E { 2d ∑ i,j=1 ∫ T 0 ∂ 2 F ∂z i ∂z j (t, S t , ν t )Γ ij (t, S t , ν t )dt } + E { ∫ T 0 ∂F ∂t (t, S t , ν t )dt + d ∑ i ∫ T 0 ∂F ∂y i (t, S t , ν t )b i (t, ν t )dt } , where z i = x i if i ≤ d and z i = y i-d if i > d.
Combining the previous equality with the Black & Scholes PDE we get

E(F (T, S T , ν T )) = F (0, S 0 , ν 0 ) + 1 2 E { ∫ T 0 2d ∑ i,j=1 [ (Γ ij -Γ ij ) ∂ 2 F ∂z i ∂z j ] (t, S t , ν t )dt } .
When σ i = σ i and the misspecified SDE (31) is different from (29) only through a different correlation matrix R, then the difference (Γ -Γ)(t, x, y) is given by the expression

(Γ -Γ)(t, x, y) = Q(t, x, y) ( R -R ) Q(t, x, y), with Q(t, x, y) = ( ( δ i-j x i √ y i ) 1≤i,j≤n 0 0 (δ i-j σ i (t, y i )) 1≤i,j≤n ) (34)
and using the trace operator tr

E(F (T, S T , ν T )) -F (0, S 0 , ν 0 ) = 1 2 E {∫ T 0 tr [ Q ( R -R ) Q∂ 2 F ] (t, S t , ν t )dt } . ( 35
)
We give now the result of the formal computation of the matrix ∂ 2 F ∂z i ∂z j (t, x, y), with

z i = x i if i ≤ d and z i = y i-d if i > d.
An example of the mathematical justifications of the derivatives used and the permutation between the differentiation operator and the expectation depend on the model chosen and can be found in the section 3.2 for the bidimensional Heston model. The different terms of the Hessian matrix of the price ∂ 2 F ∂z i ∂z j (t, x, y), are given by

∂ 2 x i ,x j F (t, x, y) = E t,x,y (∂ 2 s i ,s j f (S T )∂ x i S i T ∂ x j S j T ), (36)
∂ 2 y i ,y j F (t, x, y) = E t,x,y (∂ 2 s i ,s j f (S T )∂ y i S i T ∂ y j S j T ) + E t,x,y (∂ s i f (S T )∂ 2 y i S i T δ i-j ), ( 37
)
∂ 2 x i ,y j F (t, x, y) = E t,x,y (∂ 2 s i ,s j f (S T )∂ x i S i T ∂ y j S j T ) + E t,x,y (∂ s i f (S T )∂ 2 x i ,y i S i T δ i-j ) (38)
with the notation

f (s) = f (s 1 , s 2 , ..., s d ) and ∂ s i+d f (s) = ∂f (s) ∂s i+d = ∂f (s) ∂s i = ∂ s i f (s 1 , ..., s i , ..., s d ). (39) If the function f is convex, M ij = { E t,x,y (∂ 2 s i ,s j f (S T )∂ z j S j T ∂ z i S i T ) } i,j
, with

z i = x i if i ≤ d and z i = y i-d if i > d
, is clearly a positive matrix. Consequently, if f is convex, we can rewrite the Hessian matrix of the price as a sum of a positive matrix M and a matrix N such that

∂ 2 F ∂z i ∂z j (t, x, y) = M (t, x, y) + N (t, x, y), z i = x i i f i ≤ d and z i = y i-d i f i > d with { M ij (t, x, y) = E t,x,y (∂ 2 s i ,s j f (S T )∂ z i S i T ∂ z j S j T ), N ij (t, x, y) = (δ i-j-d + δ i-j + δ i-j+d )E t,x,y (∂ s i f (S T )∂ 2 z i ,z j S i T ) (40)
where δ represents the Kronecker delta.

Let us now focus on models based on the Heston model like the multidimensional Heston model (dimension> 2) and the multidimensional double Heston model. The choice of these models is largely due to the fact that the results established in section 3.2 for the bidimensional Heston model can be easily extended to these models. However, the extension to a larger class of models is conceivable but will request other techniques to overcome some theoretical problems. For example, the assumption (A1) (in section 2) is an important point in the proofs given in sections 3.2 and 4.

As already mentioned, the correlation structure chosen for the bidimensional Heston model does not include all the configurations. The extension models considered here will have the same kind of correlation structure used for the bidimensional Heston model in ( 5), ( 6), [START_REF] Keller | Moment explosions and long-term behavior of affine stochastic volatility models[END_REF] and [START_REF] Del Bano Rollin | A new look at the Heston characteristic function[END_REF], that is to say, we correlate each pair of stocks (S i T , S j T ) independently by a coefficient ρ ij and we propagate this correlation on the volatilities (ν i T , ν j T ) thanks to ρ i and ρ j which are known from the one-dimensional calibration.

The idea here is first to check that trace

[(Q(t, S t , ν t )∆RQ(t, S t , ν t )) N ] = 0. A sufficient condition is to have a matrix ∆R orthogonal to the matrix O ij = (δ i-j-d +δ i-j +δ i-j+d ) 1≤i,j≤2d
in the sense of the bilinear symmetric form Φ(A, B) = trace(AB). This condition is fulfilled by all symmetric matrices that have zeros on the diagonal of the four blocks

∆R =                  0 × ... × × 0 × ... × × × 0 × ... × × 0 × ... × . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . × ... × 0 × × ... × 0 × × × ... × 0 × × ... × 0 0 × ... × × 0 × ... × × × 0 × ... × × 0 × ... × . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . × ... × 0 × × ... × 0 × × × ... × 0 × × ... × 0                  . ( 41
)
Regarding the multi-asset Heston model, as it will be done for the two-dimensional case, if we fix the correlation between each asset and its volatility we easily obtain a matrix ∆R similar to (41). Consequently, if the misspecified asset vector S differs only from the market asset S by ρ ij , the difference quotient (35) becomes

E(f (S T )) -E(f (S T )) ∆ρ ij = E {∫ T 0 S i t S j t √ ν i t ν j t E t,St,νt [ ∂ 2 s i ,s j f (S T ) Ṡi T Ṡj T α i t,T α j t,T ] dt } , ( 42
)
where α i t,T and α j t,T have similar values as α 1 t,T and α 2 t,T given later in (48). The same idea can be used for multi-asset models based on the double Heston model, indeed each stock i has the following dynamic

dS i t = S i t ( √ ν i1 t dZ i1 t + √ ν i2 t dZ i2 t ) , S i 0 = x i 0 , dν i1 t = κ i 1 (θ i 1 -ν i1 t )dt + η i 1 √ ν i1 t d Z i1 t , ν i1 0 = y i1 0 , dν i2 t = κ i 2 (θ i 2 -ν i2 t )dt + η i 2 √ ν i2 t d Z i2 t , ν i2 0 = y i2 0 , d ⟨ Z i1 , Z i2 ⟩ t = d ⟨ Z i1 , Z i2 ⟩ t = 0, d ⟨ Z i1 , Z i1 ⟩ t = ρ i1 dt, d ⟨ Z i2 , Z i2 ⟩ t = ρ i2 dt.
Because of the decorrelation of (Z 1i , Z 1i ) and (Z 2i , Z 2i ), if (ρ i1 , ρ i2 ) are already known using the one-dimensional calibration for each stock i, we obtain a matrix ∆R similar to (41) which allows to have a difference quotient analogous to (42).

Differentiability of the price and studying some specific cases

We suppose that the misspecified price of the asset vector is also given by ( 5), ( 6), ( 7) and ( 8) but with different inter-asset correlation ρ, that is to say, the only misspecified parameter is the inter-asset correlation. Thus, the difference (Γ -Γ)(t, x, y) is given as in (34) with

R -R = (ρ -ρ)     0 1 0 ρ 2 1 0 ρ 1 0 0 ρ 1 0 ρ 1 ρ 2 ρ 2 0 ρ 1 ρ 2 0     and
Q(t, x, y) =     x 1 √ y 1 0 0 0 0 x 2 √ y 2 0 0 0 0 η 1 √ y 1 0 0 0 0 η 2 √ y 2     , x = (x 1 , x 2 ) y = (y 1 , y 2 ) .
The matrix N given in ( 40) is orthogonal to R -R by the trace operator and thus it is also orthogonal to (Γ -Γ). In fact

N (R -R) ρ -ρ =     0 ρ 1 D 1 0 ρ 2 ρ 1 D 1 ρ 2 D 2 0 ρ 1 ρ 2 D 2 0 0 D ′ 1 0 ρ 2 D ′ 1 D ′ 2 0 ρ 1 D ′ 2 0     with D 1 = E t,x,y ( ∂ s 1 f (S T )∂ 2 x 1 ,y 1 S 1 T ) , D ′ 1 = E t,x,y ( ∂ s 1 f (S T ) [ ∂ 2 x 1 ,y 1 S 1 T + ρ 1 ∂ 2 y 1 ,y 1 S 1 T ]) , D 2 = E t,x,y ( ∂ s 2 f (S T )∂ 2 x 2 ,y 2 S 2 T ) and D ′ 2 = E t,x,y ( ∂ s 2 f (S T ) [ (∂ 2 x 2 ,y 2 S 2 T + ρ 2 ∂ 2 y 2 ,y 2 S 2 T ) ]) . Consequently, tr [ N (Γ -Γ)(t, x, y) ] = 0
and, with this model, (35) is reduced to

E(F (T, S T , ν T )) -F (0, S 0 , ν 0 ) = 1 2 E {∫ T 0 tr [ (Γ -Γ)M ] (t, S t , ν t )dt } , ( 43
)
where M given in (40).

Although we do get rid of the matrix N , we cannot obtain the uniqueness of ρ from (43). Indeed, even though we are happy that only the positive matrix M (positive when the payoff f is convex) remains in (43), the trace of the difference (Γ -Γ) is equal to zero which makes difficult the conclusion on the positivity of E(F (T, S T , ν T )) -F (0, S 0 , ν 0 ). This is why, in Proposition 3.2, we study only specific cases. The following proposition provides the difference quotient of the price according to ρ, here ∆ρ = ρ -ρ.

Proposition 3.1

We consider the model specified by ( 5), ( 6), ( 7) and ( 8), we make also the assumption (A1). Then, the flow derivatives

Ṡi s = S i s /x i , ( 44
)
∂ y 1 S 1 s = S 1 s γ 1 t,s , γ 1 t,s = ∫ s t ν1 u 2 √ ν 1 u dW 1 u - 1 2 ∫ s t ν1 u du, ( 45
)
∂ y 2 S 2 s = S 2 s γ 2 t,s , γ 2 t,s = ∫ s t ν2 u 2 √ ν 2 u ( ρdW 1 u + √ 1 -ρ 2 dW 2 u ) - 1 2
∫ s t ν2 u du. ( 46
)
where the CIR flow derivative νi s is either given in [START_REF] Glasserman | Monte Carlo Methods in Financial Engineering[END_REF] or replaced by its modification that vanishes once the volatility reaches zero.

Using these expressions, the difference quotient (35) becomes

E(f (S T )) -E(f (S T )) ∆ρ = E {∫ T 0 S 1 t S 2 t √ ν 1 t ν 2 t E t,St,νt [ ∂ 2 s 1 ,s 2 f (S T ) Ṡ1 T Ṡ2 T α 1 t,T α 2 t,T ] dt } ( 47
)
with S = (S 1 , S 2 ), ν = (ν 1 , ν 2 ), α 1 t,T and α 2 t,T provided by the equalities

α 1 t,T = 1 + η 1 ρ 1 γ 1 t,T = 1 + η 1 ρ 1 ( ∫ T t ν1 s 2 √ ν 1 s dW 1 s - 1 2
∫ T t ν1 s ds ) , α 2 t,T = 1 + η 2 ρ 2 γ 2 t,T = 1 + η 2 ρ 2 ( ∫ T t ν2 s 2 √ ν 2 s ( ρdW 1 s + √ 1 -ρ 2 dW 2 s ) - 1 2
∫ T t ν2 s ds ) . ( 48
)
Proof of Proposition 3.1:

The difference (Γ -Γ)(t, x, y) is equal to ∆ρ √ y 1 y 2     0 x 1 x 2 0 η 2 ρ 2 x 1 x 1 x 2 0 η 1 ρ 1 x 2 0 0 η 1 ρ 1 x 2 0 η 1 η 2 ρ 1 ρ 2 η 2 ρ 2 x 1 0 η 1 η 2 ρ 1 ρ 2 0     with x = (x 1 , x 2
) and y = (y 1 , y 2 ). Using this expression of (Γ -Γ)(t, x, y), the expression of M given in (40) and the value of the derivatives (44), ( 45) and ( 46) we get

tr [ (Γ -Γ)M ] (t, x, y) 2∆ρ √ y 1 y 2 = E t,x,y [ ∂ 2 s 1 ,s 2 f (S T )S 1 T S 2 T ( 1 + η 1 ρ 1 γ 1 t,T ) ( 1 + η 2 ρ 2 γ 2 t,T )] (49)
where the value of γ 1 and γ 2 are given in (45) and (46).

Based on the assumptions (A1) (section 2) and

(A2) |ρ| < 1, |ρ 1 | < 1, |ρ 2 | < 1,
the following theorem gives a sense to the differentiation ∂ 2 s 1 ,s 2 f (S T ) in (47) and it is based on the fact that the system of SDEs ( 7), ( 8), ( 5) and (6) driven by

A 1 t (W 1 , W 1 , W 2 , W 2 ) with A 1 =     ρ 1 √ 1 -ρ 2 1 0 0 ρρ 2 0 ρ 2 √ 1 -ρ 2 √ 1 -ρ 2 2 1 0 0 0 ρ 0 √ 1 -ρ 2 0     ,
can be rewritten thanks to the Brownian motion vector (β 1 , β 2 , β 3 , β 4 ) by setting the equality

A 1 t (W 1 , W 1 , W 2 , W 2 ) = C 1 t (β 1 , β 2 , β 3 , β 4 ) with C 1 =        1 0 0 0 ρρ 1 ρ 2 √ 1 -ρ 2 ρ 2 1 ρ 2 2 0 0 ρ 1 ρρ 2 (1-ρ 2 1 ) √ 1-ρ 2 ρ 2 1 ρ 2 2 √ 1-ρ 2 1 √ 1-ρ 2 ρ 2 2 √ 1-ρ 2 ρ 2 1 ρ 2 2 0 ρρ 1 ρ 2 (1-ρ 2 ρ 2 1 ) √ 1-ρ 2 ρ 2 1 ρ 2 2 ρ(1-ρ 2 2 ) √ 1-ρ 2 1 √ 1-ρ 2 ρ 2 1 ρ 2 2 √ 1-ρ 2 ρ 2 2 √ 1-ρ 2 √ 1-ρ 2 2 √ 1-ρ 2 ρ 2 2       
, this also implies that

t (W 1 , W 1 , W 2 , W 2 ) = A -1 1 C 1 t (β 1 , β 2 , β 3 , β 4 ). ( 50
)
Theorem 3.1 We suppose that the couples asset/volatility (S i T , ν i T ) i=1,2 have the dynamic given by ( 5), ( 6), [START_REF] Keller | Moment explosions and long-term behavior of affine stochastic volatility models[END_REF] and [START_REF] Del Bano Rollin | A new look at the Heston characteristic function[END_REF]. We also assume (A1), (A2) and f (s) = max(a 1 s 1 +a 2 s 2 -K, 0) with a 1 , a 2 ∈ (R * ) 2 . For a square integrable random variable X, the conditional expectation

E t,x,y,β 1 ,β 2 ( ∂ 2 s 1 ,s 2 f (S T )S 1 T S 2 T X ) = E t,x,y ( ∂ 2 s 1 ,s 2 f (S T )S 1 T S 2 T X (β 1 , β 2 ) t≤w≤T
) is equal to the two following values

E t,x,y,β 1 ,β 2 [ a 2 x 1 S 2 T ( K -a 2 S 2 T |a 1 | ) g 1 ( K -a 2 S 2 T a 1 x 1 S 2 T x 2 ) h ( K -a 2 S 2 T a 1 , S 2 T )] (51)
and

E t,x,y,β 1 ,β 2 [ a 1 x 2 S 1 T ( K -a 1 S 1 T |a 2 | ) g 2 ( K -a 1 S 1 T a 2 x 2 S 1 T x 1 ) h ( S 1 T , K -a 1 S 1 T a 2 )] , ( 52
) with h (s 1 , s 2 ) = E t,x,y,β 1 ,β 2 ( X S 1 T = s 1 , S 2 T = s 2 )
and

g 1 (v 1 |v 2 ) = exp ( -1 2(1-ρ 2 ) [ u 1 (v 1 ) σ 1 -ρ u 2 (v 2 ) σ 2 ] 2 ) √ 2πv 1 σ 1 √ 1 -ρ 2 1 v 1 >0 , ( 53
)
g 2 (v 2 |v 1 ) = exp ( -1 2(1-ρ 2 ) [ u 2 (v 2 ) σ 2 -ρ u 1 (v 1 ) σ 1 ] 2 ) √ 2πv 2 σ 2 √ 1 -ρ 2 1 v 2 >0 . ( 54
)
where

σ 1 = √ (1-ρ 2 1 )(1-ρ 2 ρ 2 2 ) 1-ρ 2 ρ 2 1 ρ 2 2 √ ∫ T t ν 1 s ds, σ 2 = √ (1-ρ 2 2 )(1-ρ 2 ρ 2 1 ) 1-ρ 2 ρ 2 1 ρ 2 2 √ ∫ T t ν 2 s ds, u 1 (v) = ln(v) + 1 2 ∫ T t ν 1 s ds -ρ 1 ∫ T t √ ν 1 s dβ 1 s - ρρ 2 (1-ρ 2 1 ) √ 1-ρ 2 ρ 2 1 ρ 2 2 ∫ T t √ ν 1 s dβ 2 s , u 2 (v) = ln(v) + 1 2 ∫ T t ν 2 s ds -ρρ 1 ∫ T t √ ν 2 s dβ 1 s - ρ 2 (1-ρ 2 ρ 2 1 ) √ 1-ρ 2 ρ 2 1 ρ 2 2 ∫ T t √ ν 2 s dβ 2 s , ρ = ρ √ (1-ρ 2 1 ) √ (1-ρ 2 2 ) √ (1-ρ 2 ρ 2 1 ) √ (1-ρ 2 ρ 2 2 ) ∫ T t √ ν 1 s ν 2 s ds σ 1 σ 2 .
The proof of this theorem is provided in the appendix.

Remark 3.1

1) According to section 2 (equality ( 11)), for i = 1, 2 and p ≥ 1, νi T or their modifications (once we reach τ 0 (y), we replace them by their modifications that vanish) are L p random variables and this is also the case for ( ∫ T t ν i s ds

) -1/2 thanks to the results developed in [START_REF] Dufresne | The integrated square-root process[END_REF]. Thus α 1 t,T α 2 t,T ∈ L 2 (Ω) and Theorem 3.1 tells us that the equality (49), expressed formally thanks to some elements of the matrix M (see (40)), is equal to E t,x,y (Λ t,x,y,β 1 ,β 2 ) where Λ t,x,y,β 1 ,β 2 is almost surely equal to both (51) and (52) with

h (s 1 , s 2 ) = E t,x,y,β 1 ,β 2 ( α 1 t,T α 2 t,T S 1 T = s 1 , S 2 T = s 2
) , which provides the sense of our previous use of the Dirac distribution without justification for the model specified by ( 5), ( 6), [START_REF] Keller | Moment explosions and long-term behavior of affine stochastic volatility models[END_REF] and [START_REF] Del Bano Rollin | A new look at the Heston characteristic function[END_REF].

2) The permutations of the differentiation and the expectation, that were done in the previous sections, are justified by the fact that for:

X ∈ { 1, √ y 1 γ 1 t,T , √ y 2 γ 2 t,T , √ y 1 y 2 γ 1 t,T γ 2 t,T
} either the expression:

a 2 x 1 S 2 T ( K-a 2 S 2 T |a 1 | ) h ( K-a 2 S 2 T a 1 , S 2 T ) g 1 ( K-a 2 S 2 T a 1 x 1 S 2 T x 2 ) or a 1 x 2 S 1 T ( K-a 1 S 1 T |a 2 | ) h ( S 1 T , K-a 1 S 1 T a 2 ) g 2 ( K-a 1 S 1 T a 2 x 2 S 1 T x 1
) can be dominated according to x = (x 1 , x 2 ) and y = (y 1 , y 2 ) by an L 1 -bounded random variable. Indeed, taking for example the first expression, we have first to get rid of S 2 T by a change of probability (S 2 T is a positive martingale and not only a local martingale, we refer the reader to [START_REF] Keller | Moment explosions and long-term behavior of affine stochastic volatility models[END_REF] and [START_REF] Del Bano Rollin | A new look at the Heston characteristic function[END_REF]), afterwards,

( K-a 2 S 2 T x 1 |a 1 |
) can be simplified with denominator of g 1 , finally, h can be easily dominated using the previous remark.

3 ) The assumption (A2) is necessary to have the two expressions ( 51 Now that we give a sense to all the formal expressions established previously, we provide the monotony result for some values of the products {η i ρ i } i=1,2 and {η i √ 1 -ρ 2 i } i=1,2 . Proposition 3.2 We suppose that the couples asset/volatility (S i T , ν i T ) i=1,2 have the dynamic given by ( 5), ( 6), [START_REF] Keller | Moment explosions and long-term behavior of affine stochastic volatility models[END_REF] and [START_REF] Del Bano Rollin | A new look at the Heston characteristic function[END_REF]. Assuming (A1), (A2) and a European option that has f (s) = max(a 1 s 1 + a 2 s 2 ± K, 0) as payoff, then the price is differentiable according to ρ and if:

c1) {η i ρ i } i=1,2 = 0 or c2) η 1 ρ 1 = 0, η 2 √ 1 -ρ 2 2 = 0 and 2κ 2 -η 2 ρ 2 > 0 or c3) η 2 ρ 2 = 0, η 1 √ 1 -ρ 2 1
= 0 and 2κ 1 -η 1 ρ 1 > 0, then the price is monotonous with respect to ρ. For these three cases, the price increases with respect to ρ if a 1 a 2 > 0 and decreases if a 1 a 2 < 0. Moreover, the prices of the one-dimensional calls and puts (a 1 a 2 = 0) do not depend on ρ.

Remark 3.2 • This result does not include the case {η

i √ 1 -ρ 2 i } i=1,2 = 0 because,

as we pointed out previously in Remark 3.1 3), one should have, at least, |ρ

1 | ̸ = 1 or |ρ 2 | ̸ =

to be able to use (51) or (52).

• Even though these choices are restrictive, in some cases, practitioners can found themselves using this kind of assumptions on the parameters. We refer the reader for example to [START_REF]The Volatility Surface[END_REF].

• Because the price is continuous according to η i and ρ i (because ν i is continuous according to these parameters and the payoff is continuous with respect to ν i ), we can replace the zeros in this proposition by "small values". However, we preferred not to announce this more general result because its proof is heavier and it does not help to clarify all the situations for which we have the monotony.

• From a numerical point of view, remark also that the condition (A1): 4κ i θ i > η 2 i is generally sufficient to have η i 2. Indeed, θ i represents the long term volatility and it is generally smaller than 0.4, also the mean reversion coefficient κ i used in applications can be considered smaller than 3. Subsequently,

η i |ρ i | 1 or η i √ 1 -ρ 2 i 1 is true.

Proof of Proposition 3.2:

According to (47), the domination remark 3.1.2) of the term under the double integral and the continuity of (-1, 1) ∋ r → ν s announced in Theorem 2.1 (here r = ρ), we have

∂ ρ E(f (S T )) = lim ρ→ρ E {∫ T 0 S 1 t S 2 t √ ν 1 t ν 2 t E t,St,νt [ ∂ 2 s 1 ,s 2 f (S T ) Ṡ1 T Ṡ2 T α 1 t,T α 2 t,T ] dt } = E {∫ T 0 lim ρ→ρ S 1 t S 2 t √ ν 1 t ν 2 t E t,St,νt [ ∂ 2 s 1 ,s 2 f (S T ) Ṡ1 T Ṡ2 T α 1 t,T α 2 t,T ] dt } = E {∫ T 0 S 1 t S 2 t √ ν 1 t ν 2 t E t,St,νt [ ∂ 2 s 1 ,s 2 f (S T ) Ṡ1 T Ṡ2 T α 1 t,T α 2 t,T ] dt } . ( 55
)
This then prove the differentiability of the price according to ρ when only (A1) is fulfilled.

Using formally the derivative 47), it is sufficient to prove the positivity of α i t,T . If η i ρ i = 0 then α i t,T = 1 which is sufficient to prove c1). Also if η 1 ρ 1 = 0 and η 2 √ 1 -ρ 2 2 = 0 then α 1 t,T = 1 and thanks to [START_REF] Glasserman | Monte Carlo Methods in Financial Engineering[END_REF],

∂ 2 s 1 ,s 2 f (s) = a 1 a 2 ε(a 1 s 1 +a 2 s 2 ±K) (ε is the Dirac distribution) in (
α 2 t,T = ν2 T + ( κ 1 - η 1 ρ 1 2
) ∫ T t νi s ds and provided that 2κ 2 -η 2 ρ 2 > 0, α 2 t,T > 0 which proves c2) and the proof of c3) is analogous.

Asymptotic approximation for short maturities

In this section, we remain working with the model specified by ( 5), ( 6), ( 7) and ( 8), we will establish, for short maturities, an asymptotic approximation of the derivative of the price with respect to ρ. For the sake of simplicity, we consider the option that has the following payoff

f (s 1 , s 2 ) = (s 1 -s 2 ) + .
However, the general result for the payoff of the exchange option f (s 1 , s 2 ) = (a 1 s 1 -a 2 s 2 ) + , with (a 1 , a 2 ) ∈ (R * + ) 2 , is given in Theorem 4.1 and a numerically good approximation for the spread options is given in (68). Provided that we can commute the derivative with respect to ρ and the expectation, and that the expression under the expectation is differentiable with respect to ρ (see the proof of Theorem 4.1, Step2), the derivative of the price with respect to ρ is given by

∂ ∂ρ E ( (S 1 T -S 2 T ) + ) = -E ( ∂ ρ S 2 T 1 S 1 T ≥S 2 T ) , ( 56
)
where 1 represents the indicator function. Provided that we can differentiate S 2 and ν 2 with respect to the correlation ρ (when the assumption (A0) of the section 2 is fulfilled)

∂ ρ S 2 T = S 2 T ( ∫ T 0 √ ν 2 s ( dW 1 s - ρ √ 1 -ρ 2 dW 2 s )) + S 2 T ( ∫ T 0 ∂ ρ ν 2 s 2 √ ν 2 s ( ρdW 1 s + √ 1 -ρ 2 dW 2 s ) - 1 2
∫ T 0 ∂ ρ ν 2 s ds ) . ( 57
)
Replacing the value of ∂ ρ S 2 T in (56), we get

∂ ∂ρ E ((S 1 T -S 2 T ) + ) = E ( 1 S 1 T ≥S 2 T S 2 T ( 1 2 ∫ T 0 ∂ ρ ν 2 s ds )) -E ( 1 S 1 T ≥S 2 T S 2 T ∫ T 0 ∂ ρ ν 2 s 2 √ ν 2 s ( ρdW 1 s + √ 1 -ρ 2 dW 2 s ) ) -E ( 1 S 1 T ≥S 2 T S 2 T ( ∫ T 0 √ ν 2 s ( dW 1 s - ρ √ 1 -ρ 2 dW 2 s )))
.

(58)

According to various works like the one presented in [START_REF] Keller | Moment explosions and long-term behavior of affine stochastic volatility models[END_REF] and [START_REF] Del Bano Rollin | A new look at the Heston characteristic function[END_REF], we know that S 2 T is a real positive martingale and not only a local martingale. This allows us to define a new probability measure P 2 whose density is given by

dP 2 dP = S 2 T S 2 0
. Under this new probability, Z 1 and Z 2 are two independent Brownian motions related to W 1 and W 2 by

dZ 1 t = dW 1 t -ρ √ ν 2 t dt, dZ 2 t = dW 2 t - √ 1 -ρ 2 √ ν 2 t dt.
Also, under the probability P 2 , the value of S 1 and S 2 are given by

S 1 T = x 1 exp ( ∫ T 0 √ ν 1 s dZ 1 s + ρ ∫ T 0 √ ν 1 s ν 2 s ds -1 2 ∫ T 0 ν 1 s ds ) , S 2 T = x 2 exp ( ∫ T 0 √ ν 2 s ( ρdZ 1 s + √ 1 -ρ 2 dZ 2 s ) + 1 2 ∫ T 0 ν 2 s ds
) .

By this change of probability and using (58), we obtain

E ( ∂ ρ S 2 T 1 S 1 T ≥S 2 T ) = S 2 0 E 2 ( 1 S 1 T ≥S 2 T ∫ T 0 √ ν 2 s ( dZ 1 s - ρ √ 1 -ρ 2 dZ 2 s )) + S 2 0 E 2 ( 1 S 1 T ≥S 2 T ∫ T 0 ∂ ρ ν 2 s 2 √ ν 2 s ( ρdZ 1 s + √ 1 -ρ 2 dZ 2 s ) ) . ( 59
)
For short maturities and under the assumption

(A3) ∃C ∈ R * such that ln ( a 2 x 2 a 1 x 1 ) = C √ T + o( √ T ), (a 1 , a 2 ) ∈ (R * + ) 2 ,
we will see in the proof of Theorem 4.1 that the second term of (59) can be neglected because it tends to zero with respect to T faster than the first one. Also, in Theorem 4.1, the asymptotic derivative of the price with respect to ρ is established thanks to the following lemma obtained by Ito isometry. Then

lim t→0 1 t E ( ∫ t 0 H s • dW s -H 0 • W t 2 ) = 0.
Theorem 4. [START_REF] Karoui | Robustness of the Black & Scholes formula[END_REF] We suppose that the couples asset/volatility (S i T , ν i T ) i=1,2 have the dynamic given by ( 5), ( 6), [START_REF] Keller | Moment explosions and long-term behavior of affine stochastic volatility models[END_REF] and [START_REF] Del Bano Rollin | A new look at the Heston characteristic function[END_REF]. We also make the assumptions (A0), (A2) and (A3). For short maturities, the derivative with respect to ρ of a European option that has f (s 1 , s 2 ) = (a 1 s 1 -a 2 s 2 ) + , with (a 1 , a 2 ) ∈ (R * + ) 2 as payoff can be asymptotically approximated by

∂ ∂ρ E ( (a 1 S 1 T -a 2 S 2 T ) + ) = T ∼ 0 -a 2 x 2 √ T ν 1 0 ν 2 0 2πλ exp ( - 1 2 [ C √ λ ] 2 ) + o( √ T ) (60) with λ = ν 1 0 + ν 2 0 -2ρ √ ν 1 0 ν 2 0 and the constant C comes from (A3).
From Theorem 4.1 and because a 2 > 0, it is clear that the price of an exchange option is decreasing according to ρ for short maturities.

Proof of Theorem 4.1:

We divide the proof of this theorem into two steps: In the first step, we detail the computations of (60). In the second step, we show that the commutation of the derivative with respect to ρ and the expectation in (56) is correct.

Step1: The use of the constants a 1 and a 2 is not restrictive because they can be included in the spot prices S 1 0 = x 1 and S 2 0 = x 2 . Defining the triplet of random variables (L 1 T , L 2 T , L 3 T ) by

L 1 T = 1 √ T ∫ T 0 √ ν 2 t ( dZ 1 t - ρ √ 1 -ρ 2 dZ 2 t ) , L 2 T = 1 √ T ∫ T 0 √ ν 1 t dZ 1 t , L 3 T = 1 √ T ∫ T 0 √ ν 2 t ( ρdZ 1 t - √ 1 -ρ 2 dZ 2 t ) + ln (S 2 0 /S 1 0 ) √ T - 1 √ T ∫ T 0 ( ρ √ ν 1 t ν 2 t - ν 1 t + ν 2 t 2 ) dt, (61)
as T → 0, Lemma 4.1 and the assumption (A3) allows us to have the convergence in probability of (L 1 T , L 2 T , L 3 T ) to (L 1 0 , L 2 0 , L 3 0 ) with

L 1 0 = √ ν 2 0 ( G 1 - ρ √ 1 -ρ 2 G 2 ) , L 2 0 = √ ν 1 0 G 1 L 3 0 = √ ν 2 0 ( ρG 1 - √ 1 -ρ 2 G 2 ) + C.
where G 1 and G 2 are two independent standard Normal random variables and C is the constant of the assumption (A3).

Moreover the first term of (59) is equal to

S 2 0 √ T E 2 ( L 1 T 1 L 2 T ≥L 3 T )
and thanks to both facts

P 2 (L 2 0 = L 3 0 ) = 0 and E 2 ( (L 1 T ) 2 ) = E 2 ( 1 T ∫ T 0 ν 2 t dt ) < ∞,
we obtain the convergence

E 2 ( L 1 T 1 L 2 T ≥L 3 T ) -→ T → 0 E 2 ( L 1 0 1 L 2 0 ≥L 3 0 ) . Let us compute E 2 ( L 1 0 1 L 2 0 ≥L 3 0 ) , E 2 ( L 1 0 1 L 2 0 ≥L 3 0 ) = √ ν 2 0 E 2 ( 1 A ( G 1 - ρ √ 1 -ρ 2 G 2
))

with

A = { G 1 - √ ν 2 0 (ρG 1 + √ 1 -ρ 2 G 2 ) √ ν 1 0 ≥ C } .
By the decomposition of G 1 into two independent standard Normal random variables G and G:

G 1 = ρ G + √ 1 -ρ 2 G, with G = ρG 1 + √ 1 -ρ 2 G 2 , G = √ 1 -ρ 2 ( G 1 -ρ √ 1-ρ 2 G 2 ) (62)
and A becomes

A = { G ≥ g( G) } , with g(u) = ( √ ν 2 0 -ρ √ ν 1 0 )u √ 1 -ρ 2 √ ν 1 0 + C. ( 63
)
The computation of this expectation provides

E 2 ( 1 A ( G 1 -ρ √ 1-ρ 2 G 2 )) = E 2 ( E 2 [ 1 A G √ 1 -ρ 2 G ]) = 1 √ 1 -ρ 2 E 2 ( 1 √ 2π ∫ ∞ g( G) ue -u 2 /2 du ) = 1 √ 1 -ρ 2 √ 2π E 2 ( exp { - [g( G)] 2 2
}) . ( 64
)
By finishing the calculation of the expectation and multiplying it by S 2 0 √ T ν 2 0 , we obtain the expression given in (60).

To conclude that the derivative with respect to ρ is asymptotically given by (60), it is sufficient to prove that the second term of (59) divided by √ T vanishes as T tends to zero. By Cauchy-Schwarz inequality

1 √ T E 2 1 S 1 T ≥S 2 T ∫ T 0 ∂ ρ ν 2 s 2 √ ν 2 s ( ρdW 1 s + √ 1 -ρ 2 dW 2 s ) ≤ √ I T √ E 2 ( 1 S 1 T ≥S 2 T ) (65)
with

I T = 1 T E 2 ( ∫ T 0 (∂ρν 2 s ) 2 4ν 2 s ds )
. Thanks to a conditioning with respect to B 2 and using Ito isometry, we get

I T = 1 T E 2 (∫ T 0 E 2 [(∂ ρ ν 2 s ) 2 |B 2 ] 4ν 2 s ds ) = η 2 2 ρ 2 2 4T (1 -ρ 2 ) ∫ T 0 ∫ s 0 E 2    e -κ 2 (s -u) - ( κ 2 θ 2 - η 2 2 4 ) ∫ s u dr ν 2 r    duds ≤ η 2 2 ρ 2 2 T 8(1 -ρ 2 )
.

(66)

Step2: The commutation of the derivative with respect to ρ and the expectation in (56) remains to be proven. When taking |ρ| < 1 -ϵ with 0 < ϵ ≪ 1, we can dominate the square of the random variables in the expectations E 2 of (59) by integrable random variables. The latter fact can be easily seen for the first term and regarding the second term, one should use the inequalities (65) and (66) to obtain it.

Remark 4.1

1) First, we point out that assumption (A0) is necessary to have the differentiability of ν 2 according to ρ in the strong sense which was needed in Theorem 4.1. However, it is sufficient to have 4κ 2 θ 2 > η 2 2 to use the boundedness of (66) in Step1 and Step2 of the previous proof. Also because of the differentiability of the price according to ρ (see proposition 3.2), we conjecture the validity of the asymptotic approximation when the assumption (A0) is replaced by (A1).

2) Rewriting the asymptotic approximation (60) without the constant C of the assumption (A3), we get

-a 2 x 2 √ T ν 1 0 ν 2 0 2πλ exp    - 1 2   ln ( a 2 x 2 a 1 x 1 ) √ T λ   2    , with λ = ν 1 0 + ν 2 0 -2ρ √ ν 1 0 ν 2 0 .
Although this approximation works well for T ≤ 0.2, we find out numerically that the expression

-a 2 x 2 √ T ν 1 0 ν 2 0 2πλ exp    - 1 2   ln ( a 2 x 2 a 1 x 1 ) √ T λ + √ T λ 2   2    (67)
allows us to have good results even when T = 0.3. The term

√ T λ 2 comes from the finite variation process 1 √ T ∫ T 0 ( ρ √ ν 1 t ν 2 t - ν 1 t +ν 2 t 2
) dt in the expression of L 3 T in (61).

3) As we will see in section 5.1, the expression (67) provides a good estimation of the derivative with respect to ρ for exchange options with maturities T ≤ 0.3. For short maturity, using the following approximations

E ((a 1 S 1 T -a 2 S 2 T + K) + ) ≃ 0.5 * E ( (a 1 S 1 T -a 2 S 2 T ) + ) + 0.5 * E ( (a 1 S 1 T -a 2 S 2 T ) + ) with S 1 T = ( x 1 /x 1 )S 1 T , S 2 T = ( x 2 /x 2 )S 2 T and x 1 = x 1 + K a 1 , x 2 = x 2 -K a 2 .
and applying (67) on these approximations, we obtain another good estimation of the derivative of the spread options with respect to ρ, given by

-a 2 x 2 2 √ T ν 1 0 ν 2 0 2πλ exp   - 1 2 [ ln (a 2 x 2 /a 1 x 1 ) √ T λ + √ T λ 2 ] 2   + -a 2 x 2 2 √ T ν 1 0 ν 2 0 2πλ exp   - 1 2 [ ln (a 2 x 2 /a 1 x 1 ) √ T λ + √ T λ 2 ] 2   , ( 68
)
with λ = ν 1 0 + ν 2 0 -2ρ √ ν 1 0 ν 2 0 .
4) Finally, for short maturities we point out that using models based on Heston like the twodimensional double Heston model and driving the same computations as the one done in this section, one can also obtain an approximation of the derivative of the price of an exchange option with respect to ρ.

Numerical results

From a practitioner's point of view, it is interesting to figure out the interval of maturities for which the approximation (68) (or (67)) is acceptable and to estimate, thanks to a Monte Carlo simulation, the value of the errors produced by this approximation. In addition to that, because the monotony result is established for some values of η i , ρ i and √ 1 -ρ 2 i , it is important to show, at least numerically, that the practical values of these parameters ensure the monotony.

When using Monte Carlo, in order to check the monotony of the price according to ρ, one has to decrease significantly the variance of the simulations by using as many trajectories as possible. The latter fact is even more true for the approximation of the derivative with respect to ρ using Monte Carlo. In all the implemented simulations we make sure that the obtained results are, at least, ten times bigger than the error induced by the 95% confidence interval1 . To reach this high accuracy Monte Carlo simulation in an acceptable execution time, we simulated M = 2 22 trajectories on an Nvidia 480 GTX GPU (Graphics Processing Unit).

The reader may have noticed that the correlation structure, used in ( 5), ( 6), ( 7) and ( 8), does not allow the model to be affine. Consequently, we cannot use, for instance, the Alfonsi discretization scheme [START_REF] Alfonsi | High order discretization schemes for the CIR process: Application to affine term structure and heston models[END_REF] for the Monte Carlo simulations. Nevertheless, for the volatilities, we implement the Milstein scheme because it is known to provide good results. Indeed, as already mentioned in [START_REF]The Volatility Surface[END_REF], when the assumption 4κ 2 θ 2 ≥ η 2 2 is fulfilled, by setting

ν t k+1 = ( √ ν t k + η 2 √ ∆tG ) 2 + κ(θ -ν t k )∆t - η 2 4 ∆t, ∆t = t k+1 -t k , G ∼ N (0, 1)
then ν t k+1 > 0 when ν t k = 0 which reduces considerably the cases when ν t k+1 < 0. If the simulation provides ν t k+1 < 0, then it is sufficient to set ν t k+1 = 0 (for more details on the choice of discretization schemes, we refer the reader to [START_REF] Glasserman | Monte Carlo Methods in Financial Engineering[END_REF]). Besides, the assets are simulated by an Euler scheme and the discretization time δt = 0.01. Consequently, in both sections 5.1 & 5.2, the parameters of the performed simulations fulfill the assumption (A1).

Results for short maturities

This section is exclusively dedicated to testing the asymptotical derivative (68) thanks to a Monte Carlo simulation. We will consider spread options with maturities T = 0.1, 0.2, 0.3. We take the correlations ρ i ∈ {-0.85, -0.8, ..., 0.8, 0.85} such that ∆ρ = ρ i+1 -ρ i = 0.05 and we approach the derivative of the price with respect to ρ by the expression

∂ ρ F (ρ i ) = F (ρ i+1 ) -F (ρ i ) ∆ρ ( 69
)
where F (ρ i+1 ) and F (ρ i ) are the prices obtained by Monte Carlo. The resulted error between (68) and (69) will be quantified in percentage:

Error P ercentage = 100 * Expression(68) -Expression(69) Expression(69) . ( 70
)
We point out that the assumption |ρ| < 1, in Theorem 4.1, plays an important role in the precision of the approximation (69). In addition, because simulating M = 2 22 trajectories with a discretization time step δt = 0.01 is already time consuming, we have chosen to restrict ourselves to the values ρ i ∈ {-0.85, -0.8, ..., 0.8, 0.85}. Besides, after a large number of simulations, we have decided to present only the most important numerical results related to the precision of the expression (68). For example, after a large set of simulations, we concluded that ρ 1 and ρ 2 do not intervene a lot in the accuracy of the approximation (68) and we took for all figures ρ 1 = ρ 2 = -0.5 that is also a reasonable choice in practice.

We first study the impact of the model parameters on the error. This allows us to derive the "worst" cases for which the error is big. We then examine the error behavior of the approximation (68) as a function of the maturity.

The parameters that deteriorate the most the asymptotic approximation

According to Figure 1, η 1 and η 2 change barely the error produced by (68). In fact, for short maturities, using small values of η i creates bigger errors when the value of ρ is close to -1, but the average value of errors remains the same.

According to figures 2, 3 and 4, we notice that the precision of ( 68) is altered much more when κ i is big and when θ i is very different form ν i 0 . The latter fact can be explained heuristically by the mean reversion characteristic of the Heston model and because (68) does not include the action of θ i which plays quickly an important role when κ i is big.

κ 1 = κ 2 = 3.0, θ 1 = θ 2 = 0.1, η 1 = η 2 = 0.1, ν 1 0 = ν 2 0 = 0.5, x 1 = x 2 =
100 and the strike K = 0.

κ 1 = κ 2 = 3.0, θ 1 = θ 2 = 0.1, η 1 = η 2 = 0.1, ν 1 0 = ν 2 0 = 0.5, x 1 = x 2 =
100 and the strike K = 0.

κ 1 = κ 2 = 3.0, θ 1 = θ 2 = 0.1, η 1 = η 2 = 0.1, ν 1 0 = ν 2 0 = 0.5 and a 1 x 1 = a 2 x 2 = 100.
κ 1 = κ 2 = 2.0, θ 1 = θ 2 = 0.2, η 1 = η 2 = 1.2, ν 1 0 = ν 2 0 = 0.4 and x 1 = x 2 = 100.

The maturities for which the asymptotic approximation can be accepted

Now that we know the model parameters that reduce the most the accuracy of the approximation (68), we want to study the action of the payoff parameters a 1 , a 2 , S 1 0 = x 1 , S2 0 = x 2 and the strike K on the precision of the approximation (68). In figures 5, 6, 7 and 8, we have tested an extreme choice of model parameters in order to be pretty sure that the error obtained, more or less, dominates the errors gotten with standard market parameters.

From these figures, when the option is In The Money (ITM) or Out of The Money (OTM), we remark that the error increases quickly when ρ is close to 1. Although a small part 2 of the error is due to the approximation (69), the other part tells us that, when T = 0.3, ρ > 0.8 and the payoff is 20% ITM or OTM, one has to have small values of κ i (κ i ≤ 1) or small difference between θ i and ν i 0 (θ i /ν i 0 ≥ 1/2), otherwise the approximation (68) is strongly wrong. When T = 0.3, we have found out that the error percentage is always lower than 18% when either κ i ≤ 1.5 and θ i /ν i 0 ≥ 1/3 or κ i ≤ 2 and θ i /ν i 0 ≥ 1/2. The maximum error percentage associated to all these cases is lower than 18% and the average error is lower than 10%.

To sum up, with |ρ| ≤ 0.9 and ν i 0 ≤ 0.5, when

• T ≤ 0.1 and the payoff is less than 20% ITM or OTM, the approximation (68) can be accepted when θ i /ν i 0 ≥ 1/4 and κ i ≤ 3.

• T ≤ 0.2 and the payoff is less than 20% ITM or OTM, the approximation (68) can be accepted when θ i /ν i 0 ≥ 1/4 and κ i ≤ 1.5.

• T ≤ 0.3 and the payoff is less than 10% ITM or OTM, the approximation (68) can be accepted when θ i /ν i 0 ≥ 1/5 and κ i ≤ 3.

• κ i ≤ 1.5 and θ i /ν i 0 ≥ 1/3 or κ i ≤ 2 and θ i /ν i 0 ≥ 1/2, the approximation (68) can be accepted for maturities T ≤ 0.3 and payoffs less than 20% ITM or OTM.

In Figure 9, we give an example of a standard choice of parameters when η 1 and η 2 do not fulfill the Feller assumption, but we remark that we still obtain good numerical results.

Results for medium and large maturities

We have already seen, in section 3.2, that the monotony of the price according to ρ is fulfilled when η i , ρ i or √ 1 -ρ 2 i are sufficiently small. As far as ρ i and √ 1 -ρ 2 i are concerned in our successive simulations, changing the value of ρ 1 and ρ 2 did not change much numerically the rate of the monotony of the price according to ρ. Consequently, we took for all figures ρ 1 = ρ 2 = -0.5. Nevertheless, we noticed that the monotony is much stronger for small values of η i than when η i is close to 2 √ κ i θ i . What we call "Relative Increment %" in these figures is the quantity defined by

100 * F (ρ i ) -F (ρ i+1 ) F (ρ i ) , ( 71
)
where ρ i ∈ {-0.9, -0.8, ..., 0.8, 0.9} and F (ρ i ) is the price obtained by Monte Carlo.

We have preferred to simulate the value of (71), instead of the price or its derivative, for two reasons. The first one is due to the heaviness of the simulation of the derivative of the price. In fact, for T ≥ 5, to have a good Monte Carlo approximation of the derivative of the price according to ρ, one should simulate M = 2 24 trajectories and preferably use ∆ρ = 0.05 instead of 0.1. In addition to a maturity T ≥ 5 and a discretization δt = 0.01, the simulations would take an enormous time even on a GPU. The second reason comes from the fact that the monotony of the price when ρ > 0.5 is much bigger than for the other values of ρ. This behavior makes the curves almost flat when ρ ≤ -0.5 which deteriorates the monotony information.

Even though figures 10 & 11 are only illustrative, we remarked that for a maturity T ∈ [0, 10] all the prices are monotonous. In addition, the speed of this monotony decreases according to the maturity. Indeed, for maturities T ≥ 10 and ρ < -0.5, the monotony can be barely seen from prices when simulating less than M = 2 20 trajectories.

We conclude that, even though the conditions of Proposition 3.2 can be considered as restrictive, the simulation results strengthen our faith in the global monotony result of the multidimensional Heston model.

Conclusion

In this work, we tried to present, as consistent as possible, the study of the price according to the correlation. We provided a good approximation of the derivative of the price with respect to ρ for short maturities. We also saw theoretically that the monotony is fulfilled for special choices of the parameters of the model. When compared to the simulation results, the theoretical ones are a bit frustrating because we remarked numerically the clear monotony of the price according to ρ. However, only from the proofs, one can identify the important difficulties that one can face when dealing with this kind of problem. In contrast to the simulation heaviness for which the parallel GPU implementation provides serious advantages that allowed us to have solid numerical study of the monotony of the price.

x 2

)

conditionally to (β 1 w , β 2 w ) t≤w≤T , then

g(v 1 , v 2 ) = exp ( - 1 2(1 -ρ 2 ) [ u 2 1 (v 1 ) σ 2 1 + u 2 2 (v 2 ) σ 2 2 -2 ρ u 1 (v 1 )u 2 (v 2 ) σ 1 σ 2 ]) 2πv 1 v 2 σ 1 σ 2 √ 1 -ρ 2 1 v 1 >0 1 v 2 >0 . g 1 (
) = 1 x 1 x 2 ∫ R 2 + ∂ 2 v 1 ,v 2 f (x 1 v 1 , x 2 v 2 ) Φ β 1 ,β 2 (x 1 v 1 , x 2 v 2 ) g(v 1 , v 2 )dv 1 dv 2 ( * ) = - a 2 x 1 ∫ R 2 + 1 a 1 x 1 v 1 +a 2 x 2 v 2 ≥K ∂ v 1 [Φ β 1 ,β 2 (x 1 v 1 , x 2 v 2 ) g(v 1 , v 2 )] dv 1 dv 2 ( * * ) = - a 1 x 2 ∫ R 2 + 1 a 1 x 1 v 1 +a 2 x 2 v 2 ≥K ∂ v 2 [Φ β 1 ,β 2 (x 1 v 1 , x 2 v 2 ) g(v 1 , v 2 )] dv 2 dv 1 .
If a 1 > 0, equality ( * ) provides )) .

E t,x,y,β 1 ,β 2 ( ∂ 2 s 1 ,s 2 f (S T )Φ β 1 ,β 2 (S 1 T , S 2 T ) ) = a 2 x 1 ∫ R + Φ β 1 ,β 2 ( K -a 2 x 2 v 2 a 1 , x 2 v 2 ) g ( K -a 2 x 2 v 2 a 1 x 1 , v 2 ) dv 2 . Denoting g(v 2 ) = ∫ R + g(v 1 , v 2 )
If a 1 < 0, equality ( * ) provides

E t,x,y,β 1 ,β 2 ( ∂ 2 s 1 ,s 2 f (S T )Φ β 1 ,β 2 (S 1 T , S 2 T ) ) = - a 2 x 1 ∫ R + Φ β 1 ,β 2 ( K -a 2 x 2 v 2 a 1 , x 2 v 2 ) g ( K -a 2 x 2 v 2 a 1 x 1 , v 2 ) dv 2 = - a 2 x 1 ∫ R + Φ β 1 ,β 2 ( K -a 2 x 2 v 2 a 1 , x 2 v 2 ) g ( K -a 2 x 2 v 2 a 1 x 1 , v 2 ) g(v 2 ) g(v 2 ) dv 2 = -E t,x,y,β 1 ,β 2 ( E t,x,y,β 1 ,β 2 [ a 2 x 1 Φ β 1 ,β 2 ( K -a 2 S 2 T a 1 , S 2 T ) g 1 ( K -a 2 S 2 T a 1 x 1 S 2 T x 2 ) S 2 T ])
= -E t,x,y,β 1 ,β 2 )) .

The expression (51) comes from a combination of this result with the one obtained when a 1 > 0.

In the same fashion, using equality ( * * ) and (54), we obtain (52).

Figure 1 :

 1 Figure 1: The error according to η 1 and η 2 , the other parameters used are:κ 1 = κ 2 = 2.25, θ 1 = θ 2 = 0.1, ν 1 0 = ν 2 0 = 0.5, a 1 x 1 = a 2 x 2 = 100, T = 0.2 and K = 0.

Figure 2 :

 2 Figure 2: The error according to θ 1 and θ 2 , the other parameters used are:κ 1 = κ 2 = 2.25, η 1 = η 2 = 0.4, ν 1 0 = ν 2 0 = 0.5, a 1 x 1 = a 2 x 2 = 100, T = 0.2 and K = 0.

  Error Percentage κ 1 = 2.25, κ 2 = 2.25 κ 1 = 2.25, κ 2 = 1.00 κ 1 = 1.00, κ 2 = 1.00

Figure 3 :

 3 Figure 3: The error according to κ 1 and κ 2 , the other parameters used are: θ 1 = θ 2 = 0.1, η 1 = η 2 = 0.4, ν 1 0 = ν 2 0 = 0.5, a 1 x 1 = a 2 x 2 = 100, the maturity T = 0.2 and K = 0.

1 =

 1 κ 2 = 1.00, θ 1 = θ 2 = 0.30κ 1 = κ 2 = 2.25, θ 1 = θ 2 = 0.10

Figure 4 :

 4 Figure 4: The error according to {θ i } i=1,2 and {κ i } i=1,2 , the other parameters used are: η 1 = η 2 = 0.4, ν 1 0 = ν 2 0 = 0.5, a 1 x 1 = a 2 x 2 = 100, the maturity T = 0.2 and the strike K = 0.

Figure 5 :

 5 Figure5: The error percentage for a maturity T = 0.1 when changing a 1 /a 2 , the parameters used are:κ 1 = κ 2 = 3.0, θ 1 = θ 2 = 0.1, η 1 = η 2 = 0.1, ν 1 0 = ν 2 0 = 0.5, x 1 = x 2 =100 and the strike K = 0.

Figure 6 :

 6 Figure6: The error percentage for a maturity T = 0.3 when changing a 1 /a 2 , the parameters used are:κ 1 = κ 2 = 3.0, θ 1 = θ 2 = 0.1, η 1 = η 2 = 0.1, ν 1 0 = ν 2 0 = 0.5, x 1 = x 2 =100 and the strike K = 0.

Figure 7 :

 7 Figure 7: The error percentage for a maturity T = 0.1 when changing the strike K, the parameters used are:κ 1 = κ 2 = 3.0, θ 1 = θ 2 = 0.1, η 1 = η 2 = 0.1, ν 1 0 = ν 2 0 = 0.5 and a 1 x 1 = a 2 x 2 = 100.

Figure 8 :

 8 Figure 8: The error percentage for a maturity T = 0.3 when changing the strike K, the parameters used are:κ 1 = κ 2 = 3.0, θ 1 = θ 2 = 0.1, η 1 = η 2 = 0.1, ν 1 0 = ν 2 0 = 0.5 and a 1 x 1 = a 2 x 2 = 100.

1 Figure 9 :

 19 Figure9: The error percentage for 20% ITM or OTM contracts, the parameters used are:κ 1 = κ 2 = 2.0, θ 1 = θ 2 = 0.2, η 1 = η 2 = 1.2, ν 1 0 = ν 2 0 = 0.4 and x 1 = x 2 = 100.

Figure 10 :

 10 Figure10: The relative increment % for a maturity T = 5 when changing a 1 /a 2 , the parameters used are:η 1 = η 2 = 1.5, κ 1 = κ 2 = 3.0, θ 1 = θ 2 = 0.2, ν 1 0 = ν 2 0 = 0.4, x 1 = x 2 =100 and the strike K = 0.

Figure 11 :

 11 Figure11: The relative increment % for a maturity T = 5 when changing the strike K, the parameters used are:η 1 = η 2 = 1.5, κ 1 = κ 2 = 3.0, θ 1 = θ 2 = 0.2, ν 1 0 = ν 2 0 = 0.4 and a 1 x 1 = a 2 x 2 = 100.

1 ( K -a 2 S 2 T a 1 x 1 S 2 T

 1212

  ) and (52). Indeed, for instance if |ρ 1 | = 1, |ρ 2 | < 1 and |ρ| < 1 then the expression (52) still can be used but (51) cannot.

4) Although Theorem 3.1 considers that

f (s) = max(a 1 s 1 + a 2 s 2 -K, 0) with a 1 , a 2 ∈ (R * ) 2 , the result for f (s) = max(a 1 s 1 + a 2 s 2 + K, 0) with a 1 , a 2 ∈ (R * )

2 can be easily derived in the same way. When dealing with f (s) = max(a 1 s 1 + a 2 s 2 -K, 0), a 1 and a 2 can be both positive and, subsequently, the result of Theorem 3.1 can be applied on contracts beyond the spread options.

  Lemma 4.1 On R d , we define a Brownian motion W t and [0, ∞) ∋ t → H t ∈ L 2 (R d ) an adapted random process such that lim

	t↓0	E	(	||H t -H 0 || 2 )	= 0, where || • || is the Euclidean norm.

  v 1 |v 2 ) and g 2 (v 2 |v 1 ), given in (53) and (54), are the conditional densities respectively to ,s 2 f(S T )Φ β 1 ,β 2 (S 1 T , S 2 T )|{β 1 w , β 2 w } t≤w≤T

	S 2 T x 2 = v 2 and to	S 1
				T S 2 T X	)	is equal to
	E t,x,y	(	∂ 2 s 1

T x 1 = v 1 . Besides, if we denote Φ β 1 ,β 2 (s 1 , s 2 ) = s 1 s 2 E ( X S 1 T = s 1 , S 2 T = s 2 , {β 1 w , β 2 w } t≤w≤T ) then, setting v i = s i /x i , E t,x,y,β 1 ,β 2 ( ∂ 2 s 1 ,s 2 f (S T )S 1

  dv 1 , we obtainE t,x,y,β 1 ,β 2 ( ∂ 2 s 1 ,s 2 f (S T )Φ β 1 ,β 2 (S 1 T , S 2 T )

																	)
	=	a 2 x 1	∫ R +	Φ β 1 ,β 2	(	K -a 2 x 2 v 2 a 1	, x 2 v 2	)	g	(	K -a 2 x 2 v 2 a 1 x 1	, v 2	)	g(v 2 ) g(v 2 )	dv 2
	= E t,x,y,β 1 ,β 2	(	E t,x,y,β 1 ,β 2	[	a 2 x 1	Φ β 1 ,β 2	(	K -a 2 S 2 T a 1	, S 2 T	)	g 1 (	K -a 2 S 2 T a 1 x 1	T x 2 S 2	)	T S 2	])
	= E t,x,y,β 1 ,β 2	(	a 2 x 1	S 2 T	(	K -a 2 S 2 T a 1	)	g 1 (	K -a 2 S 2 T a 1 x 1	S 2 T x 2	)	h	(	K -a 2 S 2 T a 1	, S 2 T

The difference F (ρ i+1 ) -F (ρ i ) defined in (69) is at least, ten times bigger than the error induced by the 95% confidence interval.

When simulating M = 2 26 trajectories and using ∆ρ = we found out that the maximum error attained in Figure5is 24% instead of 28%.

Appendix Proof of Theorem 3.1:

All the computations will be done thanks to the fact that the couple (