Eric Benhamou
email: eric.benhamou@pricingpartners.com

Emmanuel Gobet
email: emmanuel.gobet@imag.fr.

Mohammed Miri
email: mohammed.miri@pricingpartners.com.

Analytical formulas

Keywords: asymptotic expansion, local volatility model, HJM framework, Hull and White model, Malliavin calculus, small diffusion process, CEV model JEL Classification: G13 Mathematics Subject Classification (2000): 60J75, 60HXX

published or not. The documents may come

Introduction

Long term callable path dependent equity options have generated new modeling challenges as the path dependency requires consistency in the asset diffusion while the early exercise on long period suggests interest rates risk. To appropriately account both for the asset diffusion consistency and the interest rates risk, we consider in this paper a local volatility model with stochastic interest rates. Recent works have mostly focused on extending stochastic volatility models to stochastic interest rates, as described in Piterbarg (2005a), [START_REF] Balland | Stoch-vol model with interest rate volatility[END_REF], [START_REF] Andreasen | Closed form pricing of FX options under stochastic interest rates and volatility[END_REF] or [START_REF] Haastrecht | Pricing long-maturity equity and FX derivatives with stochastic interest rates and stochastic volatility[END_REF]. However, very few works have been done on extending local volatility models to stochastic interest rates, except some work on the explicit bias between the Dupire local volatilities in a stochastic and deterministic interest rate models (see [START_REF] Benhamou | Stochastic interest rates for local volatility hybrid models[END_REF]).

Local volatility enables to infer a diffusion process, which is consistent with the whole volatility surface as explained in [START_REF] Dupire | Pricing with a smile[END_REF]. But the introduction of stochastic interest rates makes the calibration process much harder: indeed, the forward PDE approach is now much more computationally expensive as the forward PDE to solve includes an additional stochastic factor due to interest rates. In order to achieve real-time pricing computations, we revisit our perturbation approach (see Benhamou et al. (2009c)) to derive approximation formulae in the case of stochastic interest rates.

As a preliminary to our computations, we briefly discuss the choice of the model for the volatility of the spot process. First, owing to the absence of arbitrage, we know that if the spot process (S t ) t and the interest rate instruments follow Itô-type dynamics, then necessarily dS t S t = r t dt + σ t dW 1 t ,

r t = f (0,t) - t 0 γ(s,t).Γ (s,t)ds + t 0 γ(s,t)dB s ,
where (r t ) t is the short term interest rate, ( f (0,t)) t is the forward rate curve at time 0 (deduced from the initial yield curve), (σ t ) t is the instantaneous volatility process, (Γ (t, T )) t is the volatility of the zero coupon bond (B(t, T )) t paying 1e at time T , γ(t, T ) = -∂ T Γ (t, T ) is the volatility of the forward rates (this is the HJM framework). The above general decomposition is written under the risk-neutral measure Q, under which W 1 and B = (B 1 , • • • , B n ) are respectively a linear and a ndimensional standard Brownian motions. So far, we have not defined the model for (σ t ) t , (Γ (t, T ) t ) t and the correlation between W 1 and B: this is the topic of the following discussion, with a focus on the pricing and calibration features. If our pricing problem were only for European options, a standard market practice is to model the forward process and to perform a change of measure choosing the forward measure as numéraire. Namely, consider the forward process F T t for the maturity date T given by F T t = S t B(t,T ) . As shown in [START_REF] Geman | Changes of numeraire, changes of probability measure and option pricing[END_REF], the pricing of a European option with final payoff ϕ(S T ) can be reformulated in the forward measures as follows

E[e -T 0 r s ds ϕ(S T )] = B(0, T )E T [ϕ(F T T )] (1.1)
where E T is the expectation under the forward measure Q T . Interestingly, the forward process (F T t ) 0≤t≤T is a martingale under the forward measure Q T , meaning that only its volatility needs to be specified (it is known that it is equal to the difference of the volatilities of (S t ) t and (B(t, T )) t up to correlation factors between W 1 and B). Thus, the equation (1.1) is illuminating as it shows that the stochastic interest rates risk seems to be eliminated from the pricing/calibration problem. However, this approach that models each forward process F T under Q T , leads to as many volatility models as the number of maturities (in addition, each volatility model is written under a specific forward measure). Furthermore, in the case of path dependent options, it is not enough to model each F T under Q T , since we can not extend the representation (1.1) using only (F T t ). Thus, we are forced to model the full dynamics of S under Q. This supports the choice of a local volatility model for S with stochastic interest rates. Being inspired by the previous arguments on modeling the forward that is martingale under the suitable forward probability, we choose to define a model on the discounted price process:

S d t = e -t 0 r s ds S t
which is also a martingale (under Q). We assume that

dS d t S d t = σ d (t, S d t )dW 1 t .
Equivalently, we study the log discounted process X t = log(S d t ) = log(S t ) -t 0 r s ds, whose dynamics is

dX t = σ (t, X t )dW 1 t - σ 2 2 (t, X t )dt, X 0 = x 0 , (1.2)
where σ (t, X t ) is the volatility term that can be related to the local volatility of the discounted process σ (t, X t ) = σ d (t, S d t ). Taking the log discounted process as a local volatility model is not very conventional as the local volatility function is now a function of the log discounted process and not of the log process itself. However, this new approach has the great advantage to remove the influence of stochastic interest rates in the local volatility function and presumably to lead to intuitive approximations. In addition, this approach leads to similar types of local volatility as the one for the forward process. Both approaches model the local volatility functions of martingale processes. For a trader accustomed to quote local volatility for the forward process, it becomes easy to shift to our approach.

To complete our preliminary discussion related to the choice of the model, it remains to specify the assumptions of the volatility of interest rates and the correlation. We consider Gaussian model for interest rates, by assuming that Γ , γ : R + × R + R n are deterministic functions (n is the number of Gaussian factors). The Brownian motions

W 1 and B = (B 1 , • • • , B n ) are correlated using deterministic functions (ρ S,r i,t ) i,t : d W 1 , B i t = ρ S,r i,t dt 1 ≤ i ≤ n.
Now, our aim is to give an analytical accurate approximation of any European option price, written as the expected value under the risk neutral probability measure of a payoff function h evaluated at the maturity time T :

A = E[e -T 0 r s ds h( T 0 r s ds + X T )], (1.3)
where h(x) = ϕ(e x ). The important cases are related to call/put for which h(x) = (e x -K) + and h(x) = (Ke x ) + . Using the zero coupon B((t, T )) t as a numéraire, one has

A = B(0, T )E T [h( T 0 r s ds + X T )], (1.4)
where E T is the expectation under the forward neutral probability Q T . The process (X t ) t has the following dynamics under the probability Q T :

dX t = σ (t, X t )dW 1,T t + (ρ S,r t .Γ (t, T )σ (t, X t ) - σ 2 2 (t, X t ))dt, X 0 = x 0 , (1.5)
where (W 1,T ) t is a Brownian motion under Q T . Black formula. An important case in our study is associated to time dependent volatility σ (t, x) = σ t for which the price is given explicitly by the Black formula for call/put payoffs. This feature relies on the Gaussian property of t 0 r s ds + X T under Q T . One has

T 0 r s ds + X T = 1 2 T 0 |Γ (t, T )| 2 dt - T 0 Γ (t, T )dB t + T 0 f (0,t)dt + log(S 0 ) + T 0 σ t dW 1 t - 1 2 T 0 σ 2 t dt = log( S 0 B(0, T ) ) + T 0 σ t dW 1,T t - T 0 Γ (t, T )dB T t - 1 2 T 0 |Γ (t, T )| 2 dt - 1 2 T 0 σ 2 t dt + T 0 σ t ρ S,r t .Γ (t, T )dt, (1.6)
where (W 1,T , B T ) is a Q T -Brownian motion (with the same correlation than under Q). Denote by σ Black the equivalent Black volatility defined by

(σ Black ) 2 T = T 0 [σ 2 t + |Γ (t, T )| 2 -2σ t ρ S,r t .Γ (t, T )]dt.
Thus, T 0 r s ds + X T is Q T -distributed as a Gaussian r.v. with mean log( S 0 B(0,T ) ) -1 2 (σ Black ) 2 T and variance (σ Black ) 2 T . In particular for call options (h(x) = (e x -K) + ), it follows that

A = B(0, T )[ S 0 B(0, T ) N ( 1 σ Black √ T log( S 0 B(0, T )K ) + 1 2 σ Black √ T ) -KN ( 1 σ Black √ T log( S 0 B(0, T )K ) - 1 2 σ Black √ T )].
General local volatility models. Then, to obtain the analytical approximation for general local volatility models, we follow the ideas from the paper Benhamou et al. (2009c) and we introduce a parameterized process given by:

dX ε t = ε(σ (t, X ε t )dW 1,T t + (ρ S,r t .Γ (t, T )σ (t, X ε t ) - σ 2 2 (t, X ε t ))dt), X ε 0 = x 0 , (1.7)
where the parameter ε lies in the range [0, 1]. Obviously, this parameterized process is equal to the initial one for ε = 1. Remarkably, it is much easier to calculate the price (1.4) as an expansion formula with respect to ε (this is related to Black formula).

Once we have derived all the terms of the expansion, we see that the price of the European option is obtained by taking ε = 1 in the expansion.

Before giving the details of this expansion, in Section 2 we expose the general methodology underlying our approximations. This is based on the idea of using a proxy model (in the current paper, this is the Black model). Actually, the proxy model is strongly connected to the expansion of the model w.r.t. some parameter ε (see Equation (1.7)). From this point of view, our approach seems to be very close to Watanabe's expansion [START_REF] Watanabe | Analysis of Wiener functionals (Malliavin calculus) and its applications to heat kernels[END_REF] of Wiener functionals. However, this is different and in Section 2, we emphasize these differences. In Section 3, we state the assumptions used in the paper and define the notation. Then we give the main results (Theorems 3.3 and 3.6). An extension to commodities (local volatility with stochastic convenience yield) is performed in Section 4. In Section 5, we present some numerical results to illustrate the accuracy of the formulas. Proofs of main results are given in Section 6.

2 Our methodology and comparison with Watanabe's approach

Methodology

There are many expansion methods. See [START_REF] Labordère | A general asymptotic implied volatility for stochastic volatility models[END_REF] for geodesic expansion, see [START_REF] Fouque | Derivatives in financial Markets with stochastic volatility[END_REF] for ergodic approach, see [START_REF] Lee | The moment formula for implied volatility at extreme strikes[END_REF] for extremes strikes, see Piterbarg (2005b) for parameter averaging. We refer to Benhamou et al. (2009a) for a discussion and more detailed references. Our approach is different and presumably can be applied in a wider framework.

To simplify the exposure of our methodology, we assume in this section that the quantity of interest is E[h(X T )] (i.e. compared to (1.4) we do not write the terms T 0 r s ds and omit to specify the probability measure). To ensure that a perturbation method can be efficiently used in practice, the quantity of interest (i.e. the price) should be decomposed as a summation of explicit terms, namely a principal part plus some correction terms. The main idea is that the principal part is given by the price in a proxy model. The distance to the proxy model is represented by a small parameter ∆ . Below, we cite some examples of proxy and ∆ that we used in our previous works.

-In Heston model Benhamou et al. (2009a), a possible proxy could be a Black Scholes model by taking the vol of vol ξ equal to 0. In that case, we could define ∆ = ξ . -In local volatility models Benhamou et al. (2009b) Benhamou et al. (2009c), the Black Scholes model could be a proxy by freezing the local volatility function σ at its initial value. In that case, ∆ could be the sup norm of the derivatives of σ .

Actually, this measure of the distance to the proxy is too rough and it needs to be specified for each model. Furthermore, the stochastic representation of the distance to the proxy is obtained by using a suitable parameterization of the model. Denoting by (X t ) t the true model for the underlying asset price (or log-asset price) and by (X P t ) t the proxy model, assume the existence of a smooth parameterization

(X ε t ) t≥0 (0 ≤ ε ≤ 1) such that X ε t | ε=1 = X t .
Then usually the proxy coincides with the first order Taylor expansion

X P t = X ε t | ε=0 + ∂ ε X ε t | ε=0 .
This step is very model specific, it is not possible to write a general theory and we refer for the details to the introduction of our works Benhamou et al. (2009b), Benhamou et al. (2009c) Benhamou et al. (2009a), the PhD manuscript of the third author (Miri (2009)) or the next section for the current model. Now, at least formally we can write a Taylor expansion for the price of a vanilla option h on the asset X maturing at time T :

E[h(X T )] = E[h(X P T )] + E[h (1) (X P T )(X T -X P T )] + • • • + E[h ( j) (X P T ) (X T -X P T ) j j! ] + Resid j .
(2.1)

The leading term E[h(X P T )] is explicit, because usually the prices in the proxy model are explicit (actually, it is somehow a constraint in the choice of our proxy model).

The second term

E[h (1) (X P T )(X T -X P T )
] plays the role of the first correction term, but it is not explicit for general processes and has to be approximated. The trick is to decompose this term as a summation of Greeks (in the proxy model) plus a residual term:

E[h (1) (X P T )(X T -X P T )] = n 1 ∑ i=1 a (1) i ∂ i x E[h(X P T + x)]| x=0 + O(∆ k ).
(2.2)

The above error magnitude O(∆ k ) (for a given k > 0) is the expected global accuracy of our final formula. The Greeks in the proxy model arising in (2.2) have to be explicit as well. The derivation of this decomposition is strongly model-dependent and we derive it using Malliavin calculus combined with the parameterization (X ε t ) t≥0 . For the convenience of the reader, we shall explain why Greeks appear naturally in (2.2). The identification using Greeks can be seen as an inverse procedure used in the literature about integration by parts formula and Malliavin calculus [START_REF] Fournié | Applications of Malliavin calculus to Monte Carlo methods in finance[END_REF]). Indeed, we know that

∂ i x E[h(X P T + x)]| x=0 = E[h (i) (X P T )] = E[h(X P T )H i T ]
for some "Malliavin weights" H i T and to get (2.2), we should identify (a

(1) i ) i such that ∑ n i=1 a (1) i H i T ≈ (X T -X P T )
. Then, we repeat the decomposition (2.2) for each term of

the type E[h ( j) (X P T ) (X T -X P T ) j j! ], by writing E[h ( j) (X P T ) (X T -X P T ) j j! ] = ∑ n j i=1 a ( j) i ∂ i x E[h(X P T + x)]| x=0 + O(∆ k ).
Finally, we bring together all these contributions in (2.1) up to the minimal order j k verifying Resid j k = O(∆ k ). This gives

E[h(X T )] = E[h(X P T )] + j k ∑ j=1 [ n j ∑ i=1 a ( j) i ∂ i x E[h(X P T + x)]| x=0 ] + O(∆ k ) = E[h(X P T )] + max j≤ j k n j ∑ i=1 ( j k ∑ j=1 a ( j) i ✶ i≤n j )∂ i x E[h(X P T + x)]| x=0 + O(∆ k ).
The general form is the following:

E[h(X T )] = E[h(X P T )] + weighted sum of Greeks ∂ i x E[h(X P T + x)]| x=0 + error.
(2.3) Of course, these formal arguments need to be mathematically clarified with many respects:

the payoff h is not smooth (the second derivative of a call payoff does not exist in the classical sense, which avoids to write directly (2.1)).

one has to carefully use Malliavin calculus to derive explicit expansions for the coefficients (a

( j) i ) i .
the estimate of error terms is a very difficult task and actually, it depends on the model and the payoff.

Comparison with Watanabe's approach

It is worth emphasizing that the Watanabe approach (Watanabe (1987)) using an asymptotic expansion of Wiener functionals is a very important reference and there exist similarities and significant differences between Watanabe's approach and ours. In his work, Watanabe considers a family (F ε ) ε≥0 of random variables defined on the Wiener space and smooth in the Malliavin sense. Below, we follow his notation. Suppose that we can write an asymptotic expansion of F ε in powers of ε at any order k:

F ε -(F 0 + εF 1 + • • • + ε k F k ) = O(ε k+1 ) as ε 0 + (2.4)
The random variables (F i ) i are smooth in the Malliavin sense and the above equality with Landau symbol has to be understood w.r.t. the Sobolev norms • D l,p . Assume additionally a uniform non degeneracy condition:

lim sup ε 0 + 1/ det(γ F ε ) p < ∞ for all p ≥ 1 (2.5)
where γ F ε is the Malliavin covariance matrix of F ε . Then, Watanabe shows that, for any function h with polynomial growth, we have:

E[h(F ε )] = E[h(F 0 )] + εE[h(F 0 )π 1 ] + • • • + ε k E[h(F 0 )π k ] + O(ε k+1 ) (2.6)
where (π i ) i is a sequence of random variables. This equality can be extended to distributions h (Theorem 2.3 in [START_REF] Watanabe | Analysis of Wiener functionals (Malliavin calculus) and its applications to heat kernels[END_REF]). At first sight, (2.3) and (2.6) are of the same type by taking X T = F ε and F 0 = X P T . Within this analogy, we may try to relate the Greeks in (2.3) with the terms (E[h(F 0 )π i ]) i in (2.6), but actually this identification is not straightforward at all.

An important feature of Watanabe's approach is that the expansion accuracy is written in terms of the small parameter ε. The parameterization w.r.t. ε is crucial (for instance, see small noise and small time expansion in [START_REF] Watanabe | Analysis of Wiener functionals (Malliavin calculus) and its applications to heat kernels[END_REF]). Last but not least, the impact of other model parameters do not enter in the estimates. And this point is a significant drawback of such computations. To illustrate this, consider the toy model

F ε = σW 1 + √ εB 1
where (W, B) is a two-dimensional Brownian motion, and σ is positive. Let us expand

E[h(F ε )] in powers of ε for h(x) = x 2 and h(x) = x + . We use that F ε is distributed as a N (0, σ 2 + ε). Clearly F 0 = σW 1 . 1. Case h(x) = x 2 .
We have

E[h(F ε )] = E[(F ε ) 2 ] = σ 2 + ε = E[h(F 0 )] + ε. 2. Case h(x) = x +
. By a scaling argument, we have:

E[h(F ε )] = σ 2 + εE[(W 1 ) + ] = σ 2 + ε σ 2 E[(σW 1 ) + ] = E[h(F 0 )] + 1 2 ε σ E[(W 1 ) + ] + O( ε 2 σ 3 ).
These two computations are coherent with the Watanabe expansion results:

E[h(F ε )] = E[h(F 0 )] + εc 1 + O(ε 2 ).
But the coefficient c 1 depends strongly on the model and on the function h (in the second case,

εc 1 = 1 2 ε σ E[(W 1 ) + ]).
In case where σ is small as well, the expansion accuracy is heavily related to ratio of ε σ and thus not only to ε. This is another strong argument against a direct application of the Watanabe asymptotic expansion. In our works, we provide non-asymptotic estimates, which crucially enables us to deduce the domain of validity of our formulas regarding all the model parameters. In other words, a significant part of our work emphasises the impact of model parameters on the approximation. This is confirmed by numerous numerical experiments (see Benhamou et al. (2009b), Benhamou et al. (2009c), Benhamou et al. (2009a)). We also show that the magnitude of the error is impacted by the payoff smoothness. The previous toy example is a convincing illustration of this phenomenon, whereas within Watanabe's approach, the payoff regularity does not play any role in the estimates.

We now discuss in more details the differences in the proofs and methodologies.

1. In our approach, we quantify the error according to the payoff smoothness: We perform a Taylor expansion for smooth options as we show before. Then, we use a regularisation method and integration by parts (Malliavin calculus) in order to upper bound the errors and express the truncated terms as a combination of Greeks of the first term. 2. As a difference, Watanabe gives asymptotic expansions of the density of the process by expressing the density as an expectation of a Dirac distribution. In Theorem 2.3 in [START_REF] Watanabe | Analysis of Wiener functionals (Malliavin calculus) and its applications to heat kernels[END_REF] for the pull-back of distributions by Wiener functionals, he uses the integration by parts formula (Malliavin calculus) applied to the distributions in order to retrieve smooth test functions. Then, he uses Taylor expansions in order to express the corrections as a function of the distribution derivatives. We guess that the expansion of Malliavin integration by parts formula is less tractable compared to our direct approach. In other words, in view of having closed formulas, it is easier to "expand and integrate by parts" than "integrate by parts and expand". 3. As a common fact within the two approaches, it is both assumed that the proxy model (X P T or F 0 ) is related to a Gaussian process or a Gaussian random variable. Actually, in our case, we also use log-normal proxys. 4. Another technical difference lies in the assumptions used for the expansions. In Watanabe results, C ∞ smoothness of the model coefficients is usually required.

In our framework, we assume optimal regularity assumptions to achieve a given accuracy (for instance, see Assumption (R N ) in Section 3). In Benhamou et al. (2009a), we also handle non smooth cases like Heston model (because of the square root in the diffusion coefficient). Finally, to complete our comparative discussions, we refer to the works [START_REF] Siopacha | Weak and strong Taylor methods for numerical solutions of stochastic differential equations[END_REF], Yoshida (1992a), Yoshida (1992b), [START_REF] Kunitomo | The asymptotic expansion approach to the valuation of interest rate contingent claims[END_REF] where Watanabe's approach is applied to computational finance.

3 Notation and main results

Notations, definitions, assumptions

The following notation will be used extensively throughout the paper.

Notation 3.1 Differentiation.

If these derivatives have a meaning, we write:

-

ψ (i) t (x) = ∂ i ψ ∂ x i (t, x) for any function ψ of two variables. -σ t = σ (t, x 0 ), σ (i) t = σ (i) (t, x 0 ). -X ε i,t = ∂ i X ε t ∂ ε i is the i th derivative of the parameterized process with respect to ε. -X i,t = ∂ i X ε t ∂ ε i | ε=0 .
These processes play a crucial role in this work.

The following notation of Greeks will be useful for interpreting the expansion terms. Notation 3.2 Greeks. Let Z be a random variable. Given a payoff function h, we define the i th Greek for the variable Z by the quantity (if it has a meaning) :

Greek h i (Z) = ∂ i E T [h( T 0 r s ds + Z + x)] ∂ x i | x=0 .
Assumptions. In order to derive tight upper bounds for our expansions, we assume that the coefficient σ is smooth enough. In what follows, N is an integer greater than 4.

-Assumption (R N ). The function σ is bounded and of class C N w.r.t x. Its derivatives up to order N are bounded.

This assumption may be restrictive because σ has to be bounded as well its derivatives. Actually, this statement is made only to simplify a bit our analysis, but we can prove that our approximation remains valid if some boundedness requirements are partially relaxed.

Notation 3.3 Function amplitudes. Under (R N ), we set M 0 =max(|σ | ∞ , • • • , |σ (N) | ∞ ), (3.1)
M 1 =max(|σ (1) | ∞ , • • • , |σ (N) | ∞ ). (3.2)
Although M 0 and M 1 may depend on N, we remove this dependence in our notation, for the sake of simplicity. In our expansion, we expect these quantities to be small.

Remark 3.1 The constant M 0 measures the amplitude of the objective function σ and its derivatives w.r.t. the second variable, whereas M 1 measures only the amplitude of its derivatives. Notice that M 1 ≤ M 0 and in case of deterministic function σ , one has M 1 = 0.

In real market, the correlation between the asset and the short rate is close to zero (see table 5.1). Therefore, the following assumption is consistent with real market data.

-Assumption (Rho). The asset is not perfectly correlated (positively or negatively)

to the interest rate:

|ρ S,r | ∞ = sup t∈[0,T ] |ρ S,r t | < 1.
To perform the infinitesimal analysis, we rely on smoothness properties not related to the payoff function itself but rather to the law of the underlying stochastic models.

-Assumption (E). The function σ does not vanish and its oscillation is bounded,

meaning 1 ≤ |σ | ∞ σ in f ≤ C E where σ in f = inf (t,x)∈R + ×R σ (t, x).
The assumption (E) is commonly called an ellipticity assumption.

Definition 3.2 As usual, we define C ∞ 0 (R) as the space of real infinitely differentiable functions h with compact support. We also define H as the space of functions with exponential growth.

As in Benhamou et al. (2009c), our analysis depend on the payoff smoothness. We split our analysis into three cases.

-Assumption (H 1 ). h belongs to C ∞ 0 (R). This case corresponds to smooth payoffs. -Assumption (H 2 ). h and h (1) belong to H . This case corresponds to vanilla options (call/put). -Assumption (H 3 ). h belongs to H . This is the case of binary options (digital).

Main results

Our perturbation approach relies on the Taylor expansion of the parameterized process (X ε t ) defined in (1.7). We have paved the way in our previous work Benhamou et al. (2009c). In the quoted reference, the parameterized process has the form

dX ε t = ε(µ(t, X ε t )dt + σ (t, X ε t )dW t )
and the aim was to approximate E[h(X 1 T )]. Hence, compared to the current study, we take a specific form for µ, namely µ(t, x) = ρ S,t t .Γ (t, T )σ (t, x) -σ 2 (t,x) 2 : with this respect, the expansion on the process (X ε t ) is very similar to that of Benhamou et al. (2009c). On the other hand, in our case, the quantity of interest is E T [h( T 0 r s ds+X 1

T )] and the extra term T 0 r s ds induces significant differences when the correction terms are computed. For the convenience of the reader, we briefly expose the computations when similar to Benhamou et al. (2009c), and we detail the arguments when new compared to Benhamou et al. (2009c).

From the definitions, X i,t ≡

∂ i X ε t ∂ ε i | ε=0 ,
we can expand the perturbed process X ε T as follows:

X ε T = X ε T | ε=0 + εX 1,T + ε 2 2! X 2,T + . . . (3.3)
Indeed, under the assumption (R 5 ), almost surely for any t, X ε t is C 4 w.r.t ε (see Theorem 2.3 in [START_REF] Kunita | Stochastic flows of diffeomorphisms. École d'[END_REF]). The diffusion dynamics of (X ε i,t ≡

∂ i X ε t ∂ ε i ) t≥0 is obtained by
a straight differentiation of the parameters of the diffusion equation of X ε . The first order term X ε 1,t is easily obtained as follows:

X ε 1,0 =0, dX ε 1,t =σ t (X ε t )dW 1,T t + (ρ S,r t .Γ (t, T )σ t (X ε t )) - σ 2 t (X ε t ) 2 )dt +εX ε 1,t (σ (1) t (X ε t )dW 1,T t + (ρ S,r t .Γ (t, T )σ (1) t (X ε t ) -σ t (X ε t )σ (1) t (X ε t ))dt). (3.4)
From the definitions, we have σ t ≡ σ (t, x 0 ) and σ (i) t ≡ σ (i) (t, x 0 ). Then, we obtain

dX 1,t =σ t dW 1,T t + (ρ S,r t .Γ (t, T )σ t - σ 2 t 2 )dt, X 1,0 = 0, (3.5) dX 2,t =2X 1,t (σ (1) t dW 1,T t + (ρ S,r t .Γ (t, T )σ (1) t -σ t σ (1) t )dt), X 2,0 = 0. (3.6)
Applying the expansion (3.3) at ε = 1, we conclude that x 0 + X 1,T is a proxy for X T . It follows the notation:

X B T = x 0 + X 1,T = x 0 + T 0 (ρ S,r s .Γ (s, T )σ s - σ 2 s 2 )ds + T 0 σ s dW 1,T s , (3.7)
where the exponent B stands for Black, which is the proxy. To obtain an approximation formula as in Benhamou et al. (2009c), we assume that h is smooth and then, we obtain approximations which are valid even if h is not smooth, which allows us to handle finally the case of arbitrary payoffs. Use the Taylor formula twice: first, for X 1 T at the second order w.r.t ε around x 0 , secondly for smooth function h at the first order w.r.t x around X B T . This leads to:

A = B(0, T )E T [h( T 0 r s ds + X T )] = B(0, T )E T [h( T 0 r s ds + X B T + X 2,T 2 + ...)] = B(0, T )(E T [h( T 0 r s ds + X B T )] + E T [h (1) ( T 0 r s ds + X B T ) X 2,T 2 ] + ...).
Note that the first term is explicit for call/put options since it is given by the Black formula previously mentioned. To achieve a fully explicit formula, it remains to transform the correction term involving X 2,T into a summation of Greeks computed in the Black proxy. This is performed using the Malliavin calculus.

Theorem 3.3 (Second order approximation price formula).

Assume that the model fulfills (R 5 ), (E) and (Rho), and that the payoff function fulfills one of the assumptions (H 1 ), (H 2 ) or (H 3 ). Then

E[e -T 0 r s ds h( T 0 r s ds + X T )] =B(0, T ) E T [h( T 0 r s ds + X B T )] + 3 ∑ i=1 α i,T Greek h i ( T 0 r s ds + X B T ) + Resid 2 , (3.8)
where

α 1,T = - T 0 (ρ S,r t .Γ (t, T )σ t - σ 2 t 2 )( T t a s σ (1) s ds)dt, α 2,T = -α 1,T -α 3,T , α 3,T = T 0 a t σ t ( T t a s σ (1)
s ds)dt, a t =σ t -ρ S,r t .Γ (t, T ).

Additionally, estimates of the error term Resid 2 is analysed according to the payoff smoothness.

-For smooth payoff (assumption (H 1 )), one has:

|Resid 2 | ≤ C sup |h ( j) | ∞ 1≤ j≤⌊ m 2 ⌋ M 1 M 2 0 ( √ T ) 3 .
-For vanilla payoff (assumption (H 2 )), one has:

|Resid 2 | ≤C( h (1) ( T 0 r s ds + X B T ) 2 + sup v∈[0,1] h (1) ( T 0 r s ds + vX T + (1 -v)X B T ) 2 ) M 0 σ in f 1 -|ρ S,r | 2 ∞ M 1 M 2 0 ( √ T ) 3 .
-For binary payoff (assumption (H 3 )), one has:

|Resid 2 | ≤C( h( T 0 r s ds + X B T ) 2 + sup v∈[0,1] h( T 0 r s ds + vX T + (1 -v)X B T ) 2 ) ( M 0 σ in f 1 -|ρ S,r | 2 ∞ ) 2 M 1 M 0 ( √ T ) 2 .
In the above estimates, the constant C depends (in an increasing way) on the bounds of the model parameters and the maturity, and the norm . 2 is the L 2 norm under the probability measure Q T .

The proof of the above Theorem is postponed to Subsection 6.2.

Remark 3.4

The above approximation is a summation of the leading term and a combination of some Greeks of the leading term:

1. B(0, T )E T [h( T 0 r s ds + X B T )] is the leading order, corresponding to the price when the parameters σ is deterministic. In the case of call/put option, it is given by the Black formula previously mentioned. For other payoffs, we can use numerical integration because the law of the random variable T 0 r s ds + X B T is Gaussian with known parameters.

B(0, T )Greek

h i ( T 0 r s ds + X B T )
is the i th derivative of the leading term w.r.t. the initial value x 0 = log(S 0 ).

3.

The coefficients α i,T are explicit and depend on the function σ , its derivative at the point x 0 , the zero coupon volatility Γ , its derivative γ and the correlation ρ S,r . In the next subsection, these constants will be made specified for some important case to enlighten more their simple expressions w.r.t. the model parameters.

As in Benhamou et al. (2009c), the accuracy of the approximation has different justifications which can be related from the error estimates on Resid 2 . Either σ is only time dependent (M 1 = 0) and the formula is exact (Black formula). Or the shorter the maturity T or the smaller the volatility (measured by M 0 ), the more accurate the approximation.

Remark 3.5 The case of the one factor Hull and white model plus time homogeneous diffusion Here, we consider the case of γ 1 (t, T ) = ξ e -κ(T -t) , a constant correlation ρ and a homogeneous volatility σ (t, x) = σ (x). Then here σ (t, x 0 ) = σ (x 0 ) ≡ σ and σ (0,1) (t, x 0 ) = σ (1) (x 0 ) ≡ σ (1) . Using Mathematica, we can compute exactly the correction coefficients. Their expressions are

α 1,T = e -2κT σ σ (1) 4κ 4 (2ρ 2 ξ 2 + 2e κT ρ(κσ (2κT + 1) + 2ρ(κT -1)ξ )ξ + e 2κT σ 2 T 2 κ 4 + ρσ (κT (3κT -2) -2)ξ κ + 2ρ 2 (κT -1) 2 ξ 2 , α 2,T = -α 1,T -α 3,T , α 3,T = e -2κT σ σ (1) ρξ + e κT σ T κ 2 + ρT ξ κ -ρξ 2 2κ 4 .

Third order approximation formula

Notice also that in real market the amplitude of volatility of the Hull and White model is ξ ≈ 1% while for the asset σ ≈ 20%. Therefore, one has presumably |Γ (.)| ∞ = O(M 2 0 ). Thus, we expect that the third order approximation formula w.r.t. M 0 does not yield additional interest rate corrections. The proof follows the arguments of Benhamou et al. (2009c) and we skip the details. Therefore, in the following higher order, we neglect terms related to additional stochastic rate corrections.

Theorem 3.6 (Third order approximation price formula). Assume that the model fulfills (R 7 ), (E) and (Rho) with |Γ (.)| ∞ = O(M 2 0 ), and that the payoff function is a vanilla payoff (assumption (H 2 )). Then

E[e -T 0 r s ds h( T 0 r s ds + X T )] =B(0, T ) E T [h( T 0 r s ds + X B T )] + 6 ∑ i=1 β i,T Greek h i ( T 0 r s ds + X B T ) + Resid ′ 3 , (3.9)
where

β 1,T =α 1,T - C 2,T 2 - C 3,T 2 - C 4,T 4 - C 5,T 4 - C 6,T 2 , β 2,T =α 2,T + C 2,T 2 + C 3,T 2 + 5C 4,T 4 + 5C 5,T 4 + 7C 6,T 2 + C 7,T 2 + C 8,T 4 , β 3,T =α 3,T -2C 4,T -2C 5,T -6C 6,T -3C 7,T - 3C 8,T 2 , β 4,T =C 4,T +C 5,T + 3C 6,T + 13C 7,T 2 + 13C 8,T 4 , β 5,T = -6C 7,T -3C 8,T , β 6,T =2C 7,T +C 8,T , and
C 2,T = ω(σ 2 , (σ (1) ) 2 ) T 0 , C 3,T = ω(σ 2 , σ σ (2) ) T 0 , C 4,T = ω(σ 2 , σ 2 , (σ (1) ) 2 ) T 0 , C 5,T = ω(σ 2 , σ 2 , σ σ (2) ) T 0 , C 6,T = ω(σ 2 , σ σ (1) , σ σ (1) ) T 0 , C 7,T = ω(σ 2 , σ 2 , σ σ (1) , σ σ (1) ) T 0 , C 8,T = ω(σ 2 , σ σ (1) , σ 2 , σ σ (1) ) T 0 .
In the above definition of the constants, the notation ω() is defined by

ω( f 1 , • • • , f k ) = T 0 f 1 (r 1 ) T r 1 f 2 (r 2 ) • • • T r k-1 f k (r k )dr 1 • • • dr k , k ≥ 0.
In addition, the error term

Resid ′ 3 is estimated as follows |Resid ′ 3 | ≤C h (1) ( T 0 r s ds + X B T ) 2 + sup v∈[0,1] h (1) ( T 0 r s ds + vX T + (1 -v)X B T ) 2 × ( M 0 σ in f 1 -|ρ S,r | 2 ∞ ) 2 M 1 M 3 0 ( √ T ) 4 +C h (1) ( T 0 r s ds + X B T ) 2 ( M 0 σ in f 1 -|ρ S,r | 2 ∞ ) 5 M 1 M 3 0 ( √ T ) 3 ,
where the constant C depends (in an increasing way) on the bounds of the model parameters and the maturity.

Proof Using an adaptation of Theorem 4.3 in Benhamou et al. (2009c), one has

E[e -T 0 r s ds h( T 0 r s ds + X T )] =B(0, T )(E[h( T 0 r s ds + X B T )] + E[h (1) ( T 0 r s ds + X B T ) X 2,T 2 ] + E[h (1) ( T 0 r s ds + X B T ) X 3,T 3! ] + E[h (2) ( T 0 r s ds + X B T ) ( X 2,T 2 ) 2 2 ] + Resid 3 ).
The error Resid 3 is estimated using an adaptation of Theorem 4.3 by :

|Resid 3 | ≤C( h (1) ( T 0 r s ds + X B T ) 2 + sup v∈[0,1] h (1) ( T 0 r s ds + vX T + (1 -v)X B T ) 2 ) ( M 0 σ in f 1 -|ρ S,r | 2 ∞ ) 2 M 1 M 3 0 ( √ T ) 4 . The first correction term E[h (1) ( T 0 r s ds + X B T ) X 2,T 2 ] is made explicit in Theorem 3.3. The other correction terms E[h (1) ( T 0 r s ds + X B T ) X 3,T 3! ] + E[h (2) (X B T ) ( X 2,T 2 ) 2 2
] are computed using integration by parts (Malliavin calculus) and a truncation argument of the weights. The truncation argument consists in neglecting the additional weights which are related to Γ since |Γ (.)| ∞ = O(M 2 0 ). Hence, the computation of these corrections terms is reduced to the computation of the additional weights of the third order correction by taking Γ ≡ 0 (this is done in Theorem 2.4 of Benhamou et al. (2009c)). By easy but tedious computations that are not detailed, we can prove that the related truncation error is estimated by:

|TruncationError| ≤ C h (1) ( T 0 r s ds + X B T ) 2 ( M 0 σ in f 1 -|ρ S,r | 2 ∞ ) 5 M 1 M 3 0 ( √ T ) 3 .
Note that the above error does not decrease to 0 as quickly as before w.r.t. the maturity T (the power is equal to 3 instead of 4).

Remark 3.7 In the case of homogeneous volatility σ (t, x) = σ (x). We write σ (t, x 0 ) ≡ σ , σ (0,i) (t, x 0 ) ≡ σ (i) . Then

C 2,T = σ 2 (σ (1) ) 2 T 2 2 , C 3,T = σ 3 σ (2) T 2 2 , C 4,T = σ 4 (σ (1) ) 2 T 3 6 = C 6,T , C 5,T = σ 5 σ (2) T 3 6 , C 7,T = σ 6 (σ (1) ) 2 T 4 24 , C 8,T = C 7,T .

Extension to stochastic dividend and convenience yield

The current framework can be easily adapted to deal with stochastic dividends or stochastic convenience yield in a local volatility model applied to commodity markets. This can be seen as an extension to Gibson Schwartz model to handle local volatility functions for example. We recall the SDE of the underlying spot in the Gibson Schwartz model (see [START_REF] Gibson | Stochastic convenience yield and the pricing of oil contingent claims[END_REF]):

dS t S t = (r t -y t )dt + σ dW 1 t , dy t = κ(α t -y t )dt + ξ t dW 2 t , d W 1 ,W 2 t = ρ t dt.
Here, the interest rate (r t ) t is deterministic. (α t ) t and (ξ t ) t are time dependent functions. Therefore, using similar modeling like for stochastic rates, we have the following framework:

dX t = σ (t, X t )dW 1 t - σ 2 (t, X t ) 2 dt, dy t = κ(α t -y t )dt + ξ t dW 2 t , d W 1 ,W 2 t = ρ t dt,
where S t = e X t e t 0 (r s -y s )ds and σ (t, x) is the local volatility function. Hence, our aim is to estimate:

e -T 0 r s ds E[h( T 0 (r s -y s )ds + X T )]. (4.1)
Analogously, the proxy X B t is:

dX B t = σ t dW 1 t - σ 2 t 2 dt, X B 0 = x 0 .
Hence, we obtain analogous corrections results:

Theorem 4.1 (Second order approximation price formula).

Assume that the model fulfills (R 5 ), (E) and (Rho), and that the payoff function fulfills one of the assumptions (H 1 ), (H 2 ) or (H 3 ). Then

e -T 0 r s ds E[h( T 0 (r s -y s )ds + X T )] = e -T 0 r s ds E[h( T 0 (r s -y s )ds + X B T )] + 3 ∑ i=1 λ i,T Greek h i ( T 0 (r s -y s )ds + X B T ) + Resid 2 , (4.2)
where

λ 1,T = T 0 σ 2 t 2 ( T t σ s σ (1) s ds)dt, λ 2,T = - T 0 b t σ t ( T t σ s σ (1) s ds)dt - T 0 σ 2 t 2 ( T t b s σ (1) s ds)dt, λ 3,T = T 0 b t σ t ( T t b s σ (1) s ds)dt, b t =σ t -ρ t ξ t T t e -κ(s-t) ds.
The error term Resid 2 is estimated as in Theorem 3.3.

Numerical Experiments

Here we give numerical examples for the accuracy of our approximation formula. As a benchmark, we use Monte Carlo methods with a variance reduction technique. We consider the one factor Hull and White model for interest rates, the CEV diffusion for the spot and constant correlation ρ. Then, γ(t, T ) = ξ e -κ(T -t) , σ (t, x) = νe (β -1)x .

In this case the correction coefficients are computed in Remark 3.5. We consider the call h(x) = (e x -K) + , ensuring that the price and the Greeks in the Black proxy are explicit.

For the correlation, we take ρ = 15%, which is a realistic value (see Table 5.1).

Monte Carlo with control variate

Using the HJM framework for the Hull and White short rate (r t ), the integrated T 0 r s ds is a Gaussian variable with mean m and variance v (see [START_REF] Brigo | Interest Rate Models: Theory and Practice[END_REF])

m = T 0 f (0,t)dt + ξ 2 2κ 2 (T + 2 κ (e -κT -1) - 1 2κ (e -2κT -1)), v = ξ 2 κ 2 (T + 2 κ e -κT - 1 2κ e -2κT - 3 2κ
).

The simulated random variable is e -T 0 r t dt (e T 0 r t dt+X T -K) + . In order to reduce the statistical error, we use a control variate method. Namely, the control variate is

e -T 0 f (0,t)dt (e T 0 f (0,t)dt+X T -K) + -E[e -T 0 f (0,t)dt (e T 0 f (0,t)dt+X T -K) + ]
. The latter expectation is approximated analytically using the third order accurate formula using lognormal proxy approximation derived in Benhamou et al. (2009c).

We take forward rates ( f (0,t) = 2%). Indeed, this choice is arbitrary, and it does not influence the accuracy or the correction terms calculus.

5.2 Accuracy of the second and third order approximation formulas (3.8), (3.9)

In Tables 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9 (corresponding to the maturities 6M, 1Y , 5Y , 10Y for a small and a large skew β = 0.8 and β = 0.2), we give detailed numerical results about the accuracy of the second order formula (3.8) and the third order formula (3.9). MC-and MC+ are the bounds of the 95%-confidence interval of the Monte Carlo estimator. Remark also that we increase the range of strike according to maturity in order to test our approximation formula for real quoted strikes. Therefore, we see that our formula (3.8) is very accurate (errors in implied volatilities are smaller1 than few bps) for β close to 1. For various values of β , we remark that our third order formula (3.9) is extremely accurate.

Table 5.2 Implied Black-Scholes volatilities for the second order formula (3.8), the third order formula (3.9) and the Monte Carlo simulations (3 millions simulations using Euler scheme with 50 time steps per year) expressed as a function of strikes at the expiry T = 6M. Parameters: β = 0.8, ν = 0.2, ξ = 0.7%, κ = 1%, ρ = 15% and x 0 = 0. 6 Proofs

Strikes

Here, we bring together the results (and their proofs) which allow us to derive the explicit terms in the formula (3.8). In the following, (u t ) (resp. (v t )) is a square integrable and predictable (resp. deterministic) process and l is a smooth function with compact support. We recall that a t = σ t -ρ S,r t .Γ (t, T ).

E T [( T 0 u t dW 1,T t )l( T 0 a t dW 1,T t )] = E T [( T 0 a t u t dt)l (1) ( T 0 a t dW 1,T t )].
If u is deterministic then

E T [( T 0 u t dW 1,T t )l( T 0 a t dW 1,T t )] = ( T 0 a t u t dt)∂ x E T [l( T 0 a t dW 1,T t + x)] x=0 .
Table 5.9 Implied Black-Scholes volatilities for the second order formula (3.8), the third order formula (3.9) and the Monte Carlo simulations (3 millions simulations using Euler scheme with 50 time steps per year) expressed as a function of strikes at the expiry T = 10Y . Parameters: β = 0.2, ν = 0.2, ξ = 0.7%, κ = 1%, ρ = 15% and x 0 = 0. where we have used Lemma 6.2 for the last equality.

Strikes

6.2 Proof of Theorem 3.3

Using Equation (1.6), the r.v. T 0 r s ds + X B T can be projected on the Q T -Brownian motion W 1,T as follows: where we have used Equation (3.5) for second Equality and Lemma 6.2 for last one. An application of Lemma 6.3 gives immediately the Equality (3.8).

T 0 r s ds + X B T = T 0 (σ t -ρ S,
Error analysis. For the smooth case, we only need estimates on ε X ε t and its derivatives, in terms of M 0 and M 1 . This is very similar to our previous work (see Benhamou et al. (2009c)) and we skip the details (use our discussion in the beginning of Section 3). For the call/put case or digital case, once again we follow the lines of the proof of Benhamou et al. (2009b) and Benhamou et al. (2009c). The computations and estimates remain the same, except for Lemma 5.5 in Benhamou et al. (2009c) and Lemma 5.3 in Benhamou et al. (2009b) which now writes as follows. Here, appears the correlation assumption (Rho). Proof We follow the same approach from the quoted references. First, we prove that a suitable Malliavin covariance matrix is invertible and we estimate the L p -norm of its inverse. Second, we apply the integration by parts from Malliavin calculus to get the existence of Z ν k , and finally, we provide estimates for its L p -norm. Only the first step is a bit different and is worth being detailed. For the other steps, we refer to the proof of Lemma 5.3 in Benhamou et al. (2009b). Let us denote F ν = T 0 r s ds + νX T + (1 -ν)X B T . All the calculus of stochastic variations will be performed relatively to the

Lemma 6. 4

 4 Assume (E),(Rho) and (R k+1 ) for a given k ≥ 1 . Let Z belong to ∩ p≥1 D k,p . For any v ∈ [0, 1], there exists a random variable Z v k in any L p (p ≥ 1) such that for any function l ∈ C ∞ 0 (R), we haveE T [l (k) ( T 0 r s ds + vX T + (1v)X B T )Z] = E T [l(v S,r | 2 ∞ σ in f √ T ) k , uniformly in v, where C depends (in an increasing way) on the bounds of the model coefficients and the maturity.

Table 5 .

 5 1 Historical correlation between assets and short term interest rate EUR. Period: 23-Sep-2007 to 22-Sep-09.

	Asset	Historical correlation
	ADIDAS	18.32%
	BELGACOM	4.09%
	CARREFOUR	7.08%
	DAIMLER	-0.94%
	DANONE	7.23%
	LVMH	4.53%
	NOKIA	4.37%
	PHILIPS	5.23%

Table 5 .

 5 3 Implied Black-Scholes volatilities for the second order formula (3.8), the third order formula (3.9) and the Monte Carlo simulations (3 millions simulations using Euler scheme with 50 time steps per year) expressed as a function of strikes at the expiry T = 6M. Parameters: β = 0.2, ν = 0.2, ξ = 0.7%, κ = 1%, ρ = 15% and x 0 = 0. Implied Black-Scholes volatilities for the second order formula (3.8), the third order formula (3.9) and the Monte Carlo simulations (3 millions simulations using Euler scheme with 50 time steps per year) expressed as a function of strikes at the expiry T = 1Y . Parameters: β = 0.8, ν = 0.2, ξ = 0.7%, κ = 1%, ρ = 15% and x 0 = 0.

	Strikes	80%	90%	100%	110%	120%
	Second Order formula	21.70%	20.94%	20.11%	19.34%	18.57%
	Third Order formula	21.95%	20.97%	20.12%	19.36%	18.70%
	MC with control variate	21.95%	20.98%	20.12%	19.37%	18.70%
	MC-	21.76%	20.91%	20.09%	19.34%	18.68%
	MC+	22.14%	21.04%	20.16%	19.39%	18.72%
	Table 5.4 Strikes	60%	80%	100%	120%	140%
	Second Order formula	20.96%	20.53%	20.10%	19.73%	19.40%
	Third Order formula	21.12%	20.55%	20.10%	19.74%	19.43%
	MC with control variate	21.14%	20.56%	20.11%	19.74%	19.44%
	MC-	18.51%	20.43%	20.07%	19.72%	19.41%
	MC+	22.56%	20.68%	20.14%	19.77%	19.47%

Table 5

 5

	Strikes	60%		80%		100%	120%	140%
	Second Order formula	22.78%		21.88%	20.22%	18.73%	16.91%
	Third Order formula	24.11%		22.08%	20.24%	18.81%	17.71%
	MC with control variate	24.13%		22.08%	20.25%	18.82%	17.72%
	MC-	23.23%		21.98%	20.21%	18.80%	17.69%
	MC+	24.90%		22.19%	20.28%	18.84%	17.74%
	6.1 Technical results related to explicit correction terms
	The two first lemmas are proved in Appendix of Benhamou et al. (2009b).
	Lemma 6.1 For any continuous (or piecewise continuous) function f , any continu-
	ous semimartingale Z vanishing at t=0, one has:
	T		T		T
	0	f t Z t dt =	0	(	t	f s ds)dZ t .

.5 Implied Black-Scholes volatilities for the second order formula (3.8), the third order formula (3.9) and the Monte Carlo simulations (3 millions simulations using Euler scheme with 50 time steps per year) expressed as a function of strikes at the expiry T = 1Y . Parameters: β = 0.2, ν = 0.2, ξ = 0.7%, κ = 1%, ρ = 15% and x 0 = 0.

Table 5 .

 5 6 Implied Black-Scholes volatilities for the second order formula (3.8), the third order formula (3.9) and the Monte Carlo simulations (3 millions simulations using Euler scheme with 50 time steps per year) expressed as a function of strikes at the expiry T = 5Y . Parameters: β = 0.8, ν = 0.2, ξ = 0.7%, κ = 1%, ρ = 15% and x 0 = 0. Implied Black-Scholes volatilities for the second order formula (3.8), the third order formula (3.9) and the Monte Carlo simulations (3 millions simulations using Euler scheme with 50 time steps per year) expressed as a function of strikes at the expiry T = 5Y . Parameters: β = 0.2, ν = 0.2, ξ = 0.7%, κ = 1%, ρ = 15% and x 0 = 0.

	Strikes	40%	70%	100%	140%	180%
	Second Order formula	21.99%	21.23%	20.55%	19.90%	19.38%
	Third Order formula	22.36%	21.27%	20.56%	19.91%	19.43%
	MC with control variate	22.44%	21.31%	20.60%	19.95%	19.47%
	MC-	21.61%	21.20%	20.55%	19.91%	19.44%
	MC+	23.17%	21.42%	20.65%	19.98%	19.50%
	Table 5.7 Strikes	40%	70%	100%	140%	180%
	Second Order formula	25.17%	23.69%	21.15%	18.50%	15.97%
	Third Order formula	28.26%	24.27%	21.27%	18.71%	17.01%
	MC with control variate	28.35%	24.29%	21.30%	18.73%	17.02%
	MC-	28.00%	24.20%	21.25%	18.70%	17.00%
	MC+	28.69%	24.38%	21.34%	18.75%	17.04%

Table 5

 5

	Strikes	30%	60%	100%	160%	220%
	Second Order formula	22.99%	22.16%	21.25%	20.38%	19.77%
	Third Order formula	23.54%	22.25%	21.27%	20.40%	19.84%
	MC with control variate	23.66%	22.32%	21.34%	20.47%	19.91%
	MC-	22.87%	22.18%	21.28%	20.43%	19.87%
	MC+	24.37%	22.47%	21.40%	20.51%	19.94%
	Lemma 6.2 One has:

.8 Implied Black-Scholes volatilities for the second order formula (3.8), the third order formula (3.9) and the Monte Carlo simulations (3 millions simulations using Euler scheme with 50 time steps per year) expressed as a function of strikes at the expiry T = 10Y . Parameters: β = 0.8, ν = 0.2, ξ = 0.7%, κ = 1%, ρ = 15% and x 0 = 0.

  Proof Applying first Lemma 6.1 to f (t) = v t and Z t = X 1,t , we obtain:

															30%	60%	100%	160%	220%
	Second Order formula			27.03%	25.69%	22.40%	18.91%	16.09%
	Third Order formula				31.69%	27.04%	22.69%	19.26%	17.35%
	MC with control variate	31.88%	27.07%	22.71%	19.27%	17.30%
	MC-														31.56%	26.96%	22.66%	19.24%	17.27%
	MC+														32.20%	27.17%	22.76%	19.30%	17.32%
	Lemma 6.3
	E T [l(		0	T	a t dW 1,T t	)(	0	T	v t X 1,t dt)]
	= (	0	T	a t σ t (		t	T	v s ds)dt)E T [l (1) (	0	T	a t dW 1,T t	)]
		+ (		0	T	(ρ S,r t .Γ (t, T )σ t -	σ 2 t 2	)(	t	T	v s ds)dt)E T [l(	0	T	a t dW 1,T t	)].
	E T [(		0	T	v t X 1,t dt)l(	0	T	a t dW 1,T t	)]
	=E T [(		0	T	(	t	T	v s ds)dX 1,t )l(	0	T	a t dW 1,T t	)]
	=E T [(		0	T	(	t	T	v s ds)(σ t dW 1,T t	+ (ρ t Γ 1 (t, T ) -	σ 2 t 2	)dt)l(	0	T	a t dW 1,T t	)]
	=(	0	T	a t σ t (	t	T	v s ds)dt)E T [l (1) (	0	T	a t dW 1,T t	)]
	+ (	0	T	(ρ t Γ 1 (t, T ) -	σ 2 t 2	)(	t	T	v s ds)dt)E T [l(	0	T	a t dW 1,T t	)],

  r t .Γ (t, T ))dW 1,T t + D T ,where D T is a Gaussian random variable independent on (W 1,T t ) t . Then using notation a t = σ t -ρ S,r t .Γ (t, T ), one getsE T [h(

	0	T	r s ds + X B T )	X 2,T 2	]
	= E T [h(	0	T	a t dW 1,T t	+ D T )	X 2,T 2	]
	= E T [h(	0	T	a t dW 1,T t	+ D T )	0	T	X 1,t (σ t dW 1,T (1) t	+ (ρ S,r t .Γ (t, T )σ t -σ t σ (1) t )dt)] (1)
	= E T [h (1) (	0	T	a t dW 1,T t	+ D T )	0	T	a t σ t X 1,t dt)] (1)
	+ E T [h(	0	T	a t dW 1,T t		+ D T )	0	T	(ρ S,r t .Γ (t, T )σ t -σ t σ (1) t )X 1,t dt], (1)

1 bp on implied volatilities is equal to 0.01%.

(n + 1)-dimensional Brownian motion ( W 1,T , B 1,T , • • • , B n,T ). We define W 1,T by the relation

from which we deduce that

The key feature in this choice is that the first component of the Malliavin derivative of T 0 r s ds is zero: D W 1,T ( T 0 r s ds) ≡ 0. Hence, we have:

and thus

It is easy to deduce that the Malliavin covariance matrix is bounded from below by:

from which it readily follows (for p ≥ 1)