Romain Bompis
email: romain.bompis@credit-agricole-sa.fr

Weak approximations for arithmetic means of geometric Brownian motions and applications to Basket options *

Keywords: Weak approximation, Geometric Brownian motion, Arithmetic mean, Malliavin calculus, Basket options AMS subject classifications. 34E10, 60Hxx

In this work we derive new analytical weak approximations for arithmetic means of geometric Brownian motions using a scalar log-normal Proxy with an averaged volatility. The key features of the approach are to keep the martingale property for the approximations and to provide new integration by parts formulas for geometric Brownian motions. Besides, we also provide tight error bounds using Malliavin calculus, estimates depending on a suitable dispersion measure for the volatilities and on the maturity. As applications we give new price and implied volatility approximation formulas for basket call options. The numerical tests reveal the excellent accuracy of our results and comparison with the other known formulas of the literature show a valuable improvement.

Introduction.

£ Motivation. In Mathematical Finance, the quick and efficient approximation of sums of Geometric Brownian motions (GBMs in short) is very useful in many problems. Among them we cite the pricing of basket options in equity models where the basket is an average of assets, the pricing of swaptions in the Libor market interest rate model where the swap rate is a stochastic convex combination of Libor rates modelled with GBMs or the forward models for commodities. For these previous examples, the multivariate log-normal model may seem restrictive because practitioners have adopted more sophisticated models as local or/and stochastic volatility models with many factors. However the classical Black-Scholes model remains still extremely popular and one can indeed be really interested in the use of equivalent multivariate formulas. Although the mathematical theory does not present any particular difficulty, the lack of tractability for the multidimensional log-normal distribution does not allow the computation of prices and hedges in closed form any more. Hence practitioners have to resort to numerical methods : numerical integration/PDE methods for low dimension or Monte Carlo simulations for high dimension. However, these methods may be to slow for real time-issues while many areas of computational finance have a crucial need of real-time, robust and accurate pricing/calibration algorithms. Then an alternative is to use analytical approximations and this is the aim of this work.

£ Methodology and main results. In this paper, adapting the Proxy principle developed in [START_REF] Benhamou | Smart expansion and fast calibration for jump diffusion[END_REF]- [START_REF] Benhamou | Expansions formulas of European options in a local volatility model[END_REF]- [START_REF] Bompis | Asymptotic and non asymptotic approximations for option valuation[END_REF]- [START_REF] Gobet | Weak approximation of averaged diffusion processes[END_REF], we provide approximations in law of arithmetic means of GBMs, using as proxy, a suitable one-dimensional GBM with an averaged volatility. We apply the results for the pricing of basket call options and obtain analytical formulas which are very accurate while remaining very fast and easy to implement. In addition, we also derive new implied volatility approximation formulas. The approximation formula for basket call options (See Theorem 13) takes the form of an explicit one-dimensional log-normal representation:

Call Basket (T, K) = Call BS (T, K)

+ i C i,T Greek BS i + error,
where the main term is the Black-Scholes formula, C i,T are weights depending on the time-dependent volatility structure of the GBMs, Greek BS i are Greeks in the Black-Scholes model well defined as soon as the log-normal proxy is non degenerate and the error is rigorously estimated w.r.t. the maturity and a dispersion measure of the volatilities. Roughly speaking, error is expected to be small when the GBMs have a close volatility structure or/and when the maturity T is short and these features are encoded in the non asymptotic error estimates of Theorems 5, 13 and 14 which emphasize the role played by the coefficients in the approximation accuracy. Besides the numerical results show an excellent accuracy, much higher than main formulas of the literature. Despite the inspiration of [START_REF] Benhamou | Smart expansion and fast calibration for jump diffusion[END_REF]- [START_REF] Gobet | Weak approximation of averaged diffusion processes[END_REF] where Gaussian proxies are used, we develop new tools to perform an expansion directly around a log-normal proxy without working on the logarithms of prices. In particular, we use a suitable interpolation between the Basket dynamic and the log-normal proxy preserving martingale properties (see subsection 2.2.1) and provide new Malliavin integration by parts formulas to compute the corrective terms (see Lemma 7).

£ Literature review. The literature is very profuse and we only summarize the main ideas related to approximations/expansions in the context of options pricing.

First regarding the asymptotic expansions, we begin with papers related to price approximation formulas. For perturbation of the PDE pricing, see the valuable reference [START_REF]Interest rate model calibration using semidefinite programming[END_REF] where approximation formulas for basket options with log-normal assets are provided though the error analysis is not handled. Regarding the stochastic analysis point of view, we refer to [START_REF] Shiraya | Pricing multi-asset cross-currency options[END_REF] in which the authors establish formally a third order asymptotic formula for the pricing of multi-asset cross-currency options using Malliavin calculus. The resulting approximation reads as an expansion around the Bachelier (Gaussian) model. We refer to [START_REF] Ju | Pricing asian and basket options via taylor expansion[END_REF] for the pricing of Asian and Basket options using Taylor expansions (without error analysis) of characteristic functions in small volatility using a one dimensional lognormal proxy matching the two first moments. Passing now to works providing implied volatility approximations, see [START_REF] Durrleman | From implied to spot volatilities[END_REF] for the shape of the implied volatility surface at the money in the small maturity asymptotic using links between the spot and implied volatility dynamics. or the strike asymptotic, see [START_REF] Avellaneda | Application of large deviation methods to the pricing of index options in finance[END_REF] where approximations of the around the money implied volatility of index options in a multidimensional time-homogeneous local volatility model are obtained using large-deviation techniques. Extensions of this work are provided in a series of papers of Bayer and Laurence [START_REF] Bayer | Small-time asymptotics for the at-the-money implied volatility in a multi-dimensional local volatility model[END_REF] (ATM options), [START_REF] Bayer | On the probability density function of baskets[END_REF] (short maturity and small volatility expansions) and [START_REF] Bayer | Asymptotics beats Monte Carlo: the case of correlated local vol baskets[END_REF] (small/large strike regime).

Last we focus on the non asymptotic expansions. In [START_REF] Carmona | Generalizing the black-scholes formula to multivariate contingent claims[END_REF], Carmona et al. take advantage of the convexity of x → x + to compute explicitly quite tight lower and upper bounds for the basket price. Another valuable reference is Piterbarg [START_REF] Piterbarg | Markovian Projection Method for Volatility Calibration[END_REF] who uses Markovian projections to approximate the dynamic of the Basket with a one dimensional local volatility model. Next an important reference is the paper of Gobet and Miri [START_REF] Gobet | Weak approximation of averaged diffusion processes[END_REF] for the approximation in law of general averaged diffusion processes using the geometric mean as proxy for the mean of exponentials. Here we provide an enhancement of the results by considering a more accurate martingale proxy and by providing a third order formula. We finally cite some works related to spread options: see [START_REF] Alexander | Analytic approximations for multi-asset option pricing[END_REF] for the pricing of rainbow options using a recursive procedure, [START_REF] Deng | Multi-asset spread option pricing and hedging[END_REF] for approximations of the exercise boundary, [START_REF] Alos | On the short-maturity behaviour of the implied volatility skew for random strike options and applications to option pricing approximation[END_REF] for the short-maturity skew behaviour using Malliavin calculus techniques or again the Thesis of Landon [START_REF] Landon | Stratégies de couverture presque optimale : théorie et applications[END_REF]Chapter 8] where the proxy principle is used in the general local volatility case but without error analysis.

£Contribution. As a comparison with these existing works, our contribution is threefold. First, for the arithmetic mean of GBMs, we provide new high order explicit approximations: (1) using a suitable one-dimensional log-normal proxy and new integration by parts formulas, [START_REF] Alos | On the short-maturity behaviour of the implied volatility skew for random strike options and applications to option pricing approximation[END_REF] preserving martingale properties, (3) with a rigorous non-asymptotic error analysis using relevant dispersion measures under a local non-degeneracy condition. Second, we apply the results to basket call options by providing explicit formulas of price and implied volatility which the shape at the money for short maturity is consistent with the results of [START_REF] Piterbarg | Markovian Projection Method for Volatility Calibration[END_REF], [START_REF] Durrleman | From implied to spot volatilities[END_REF] and [START_REF] Alos | On the short-maturity behaviour of the implied volatility skew for random strike options and applications to option pricing approximation[END_REF]. Finally, numerical investigations indicate that our formulas are very accurate with a valuable improvement in comparison to the approximations in [START_REF]Interest rate model calibration using semidefinite programming[END_REF]- [START_REF] Carmona | Generalizing the black-scholes formula to multivariate contingent claims[END_REF]- [START_REF] Landon | Stratégies de couverture presque optimale : théorie et applications[END_REF]- [START_REF] Gobet | Weak approximation of averaged diffusion processes[END_REF]- [START_REF] Ju | Pricing asian and basket options via taylor expansion[END_REF]- [START_REF] Shiraya | Pricing multi-asset cross-currency options[END_REF]. In addition, the review of the benchmark formulas and their comparison postponed to Appendix B may interest the reader.

£ Organization of the paper. The paper is organized as follows. We state the weak expansion results in section 2: first we define the setting, notations and assumptions used throughout our work, second we expose the methodology and we finally provide a third order approximation formula Theorem 5). The explicit derivation of the expansion coefficients and the error analysis are given in section 3 with complementary proofs in Appendix A. Section 4 is devoted to the application of the results to basket call options with price formulas (Theorem 13) and approximations of the implied volatility (Theorem 14). Then numerical experiments are gathered in section 5 to illustrate the excellent accuracy of our approximation formulas, taking as a benchmark the Monte Carlo method and making the comparison with the previously quoted formulas of the literature which are made explicit in Appendix B.

Weak expansion.

2.1. Setting.

Framework and proxy process.

Given a fixed time horizon T > 0 and d ≥ 1, let consider d driftless GBMs S 1 ,. . . ,S d starting from the initial value 1 to model financial assets, which are solutions of the SDEs, for any i ∈ {1, . . . , d} :

dS i,t S i,t = σ i (t) | dW t , S i,0 = 1, (1)
where (W t ) t∈[0,T ] is a standard Brownian motion (GM in short) in R q on a filtered probability space (Ω, F , (F t ) 0≤t≤T , P) with the usual assumptions on the filtration (F t ) 0≤t≤T , (σ i (t)) t∈[0,T ] is a bounded measurable function taking values in R q representing the volatility vector of the i-th GBM and . | . denotes the inner product on R q . Introducing the weights (α i ) i∈{1,...,d} ∈ [0; 1] d , such that d i=1 α i = 1, we consider the associated discrete probability measure µ defined for any x ∈ R d by µ[x . ] = d i=1 α i x i . Then we define the arithmetic mean process of the S i :

S := µ[S . ] = d i=1 α i S i , S 0 = 1, (2)
and our aim is to derive analytical approximations of E[h(S T )] for a given payoff function h at least a.e. once time differentiable with polynomial growth and a given maturity T > 0.

Remark 1. Actually our framework is equivalent to the general case with positive weights a i ≥ 0 and GBMs with drift and any positive initial value:

d S i,t / S i,t = σ i (t) | dW t + β i (t)dt, S i,0 > 0. Indeed one has g d i=1 a i S i,T = h d i=1 α i S i,T with h(S ) = g(Q T S ), Q T = d i=1 a i S i,0 e
T 0 β i (t)dt and α i := (a i S i,0 e T 0 β i (t)dt )/Q T . Hence we keep our initial setting more convenient for the subsequent analysis.

To approximate the law of S , we are looking for a proxy process. Computing the dynamic of the mean process (2), one has :

dS t = d i=1 α i S i,t σ i (t) | dW t ⇒ dS t = S t σ(t) | dW t , with: σ(t) = d i=1 S i,t S t α i σ i (t)
Then one obtains a GBM dynamic by freezing at their initial values all the stochastic processes involved in the above stochastic Basket volatility σ(t) :

dS t S t ≈ d i=1 S i,0 S 0 α i σ i (t) | dW t = d i=1 α i σ i (t) | dW t
This is a method adopted by the practitioners in the financial industry and similarly to [START_REF]Interest rate model calibration using semidefinite programming[END_REF] and [START_REF] Landon | Stratégies de couverture presque optimale : théorie et applications[END_REF]Chapter 7], one considers a martingale proxy given by the following log-normal process:

dS P t = S P t σ(t) | dW t , S P 0 = S 0 = 1 (3)
where σ(t) is the arithmetic mean of the volatility vectors :

σ(t) := µ[σ . (t)] = d i=1 α i σ i (t). (4)
But we go further because we will also provide corrective terms and we give thereafter the hypotheses stressing the limits of applicability of the weak expansion presented in subsection 2.2.1.

Remark 2. In [START_REF] Gobet | Weak approximation of averaged diffusion processes[END_REF], the geometric mean

S GM T = d i=1 S α i i,T = e T 0 σ(t)|dW t -1 2 T 0 d i=1 α i |σ i (t)| 2 dt
is used as proxy. In addition to systemically underestimating the arithmetic mean of exponentials and

S P T = e T 0 σ(t)|dW t -1 2 T 0 | d i=1 α i σ i (t)| 2 dt , S GM
T is not a martingale in the general case. That means for instance that the approximation of E S T = 1 could be inexact with the proxy S GM T , what is an undesirable feature.

Notations, definitions and assumptions.

£ Vectors, scalar products and norms.

• For a column vector v ∈ R q , we denote by v * its transpose which is a row vector.

• We recall that . | . stands for the inner product on R q and | • | for the Euclidean norm on R q .

• For a r.v. Y ∈ R m (m ≥ 1) and for p ≥ 1, ||Y|| p = (E|Y| p ) 1 p denotes its L p -norm in the L p space. £ Functions. We introduce H 1 P (R) the space of real-valued functions a.e. differentiable, having, with their first derivative, a polynomially growth. Hence, φ ∈ H

1 P (R) if ∃ C φ , C φ (1) ≥ 0 and p φ ≥ 1 such that : |φ(x)| ≤ C φ (1 + |x| p φ ), for any x ∈ R. |φ (1) (x)| ≤ C φ (1) (1 + |x| p φ -1 ), for any x ∈ R. Assuming h ∈ H 1 P (R) ensures that h(S T ) + h (1) (S T ) ∈ L 1 .
In addition, for n ∈ N, C n P (R) denotes the space of real-valued functions n-times continuously differentiable with polynomially bounded derivatives.

£ Assumption (H σ ) on the volatilities. (H σ )-i) For any i ∈ {1 . . . , d}, σ i is a bounded measurable function from [0, T ] to R q . We set:

|σ| ∞ := max i∈{1,...,d} |σ i | ∞ = max i∈{1,...,d} sup t∈[0,T ] |σ i (t)|.
One has obviously | σ| ∞ ≤ |σ| ∞ . Then one introduces two dispersion measures M σ,i and M σ:

M σ,i = |σ i -σ| ∞ and M σ = max i∈{1,...,d} M σ,i .
Clearly the following inequalities hold: M σ,i ≤ M σ and M σ ≤ 2|σ| ∞ . Errors will be quantified in terms of the control µ[(M σ,.

) k ] = d i=1 α i (M σ,i ) k , for k ≥ 1,
whose meaning is the following: if µ[M σ,. ] = 0, then, for any i ∈ {1, . . . , d}, either α i = 0 (and the contributions of S i on S and of σ i on σ are null) or σ i = σ. In all these cases, S T = S P T and the approximation is exact.

(H σ )-ii) We define V P T = T 0 | σ(t)| 2 dt
the variance of the log-normal r.v. S P T and we assume that:

V P T = T 0 | σ(t)| 2 dt > 0 (5)
In other words, the proxy S P T follows a non degenerate log-normal law. From (H σ )-ii), we easily deduce the next proposition: Proposition 3. Assume (H σ ) and that h ∈ H 1 P (R). One defines 1 the n-th Greek for h:

G h n := ∂ n S n 0 E[h(S P T )] = ∂ n ε n E h (1 + ε)S P T ε=0 . (6)
One has the following control, ∀n ≥ 1 :

|G h n | ≤ C n C h (1) [V P T ] -n-1 2 , (7)
for a constant C n > 0 depending on a non-decreasing way on |σ| ∞ , T and on the ellipticity ratio

|σ| 2 ∞ T V P T .
Proof. Consider the Gaussian density D(y) := expy 2 /(2V P T ) / 2πV P T for any y ∈ R, to get readily, transferring the higher derivatives into the Gaussian density:

G h n = ∂ n ε n R h (1 + ε)e y-1 2 V P T D(y)dy ε=0 = R h (1) (e y-1 2 V P T )e y-1 2 V P T ∂ n-1 ε n-1 D y -ln(1 + ε) ε=0 dy.
The proof is completed using standard upper bounds for the derivatives of D.

1 well defined as soon as h has polynomial growth and S P T is non-degenerate, i.e. our assumption (H σ ).

£ Generic constants and upper bounds. In the derivation of error estimates, we keep the same notation c for all non-negative constants depending on: universal constants, the index p considered for L p -norms, the growth parameters C h , C h (1) , p h of h, in a non decreasing way on the model parameters M σ, |σ| ∞ , the ellipticity ratio |σ| 2 ∞ T/V P T and T . The constants c remain bounded as these dependence parameters go to 0 and do not depend on µ and d. To state the error estimate, one uses the following notations: A ≤ c B for positive A, meaning A ≤ cB for a generic constant c and "A = O(B)" standing for |A| ≤ c B.

Approximation

dS η i,t = S η i,t σ η i (t) | dW t , S η i,0 = 1, σ η i (t) = ησ i (t) + (1 -η) σ(t), (8)
where η is an interpolation parameter lying in [0, 1] averaging the volatilities σ i and σ. One also defines the convex combination of the interpolated processes S η = d i=1 α i S η i which is a martingale for any η ∈ [0, 1] and equals S for η = 1 on the one hand, and coincide with S P for η = 0 on the other hand.

Remark 4. Our parametrization differs from the one used in [START_REF] Landon | Stratégies de couverture presque optimale : théorie et applications[END_REF]Chapter 7] denoted by S η t :

S η t = d i=1 α i e t 0 σ η i (s)|dW s -1 2 t 0 |ησ i (t)+(1-η) σ(t)| 2 ds S η t = d i=1 α i e t 0 σ η i (s)|dW s -1 2 t 0 η|σ i (s)| 2 +(1-η)| σ(s)| 2 ds ,
which is not a martingale in general for η ∈]0, 1[. Consequently the resulting approximation may suffer from numerical arbitrage. We think that preserving the martingale property (serving as a base for call/put parity relationship for instance) is crucial and this a benefit of our work.

Regarding the closed form of S η i,t in [START_REF] Benhamou | Expansions formulas of European options in a local volatility model[END_REF], one has that almost surely for any t,

η → S η i,t is C ∞ ([0, 1], R).
Besides it is more convenient to differentiate the SDEs satisfied by S η i (see [START_REF] Kunita | Stochastic flows and stochastic differential equations[END_REF]) and then to resolve the resulting SDEs. One obtains using the Leibniz formula for the derivatives of S η i,t w.r.t. η denoted by S η,(k) i,t (with the convention S η,(0) i,t

:= S η i,t ) that for any k ≥ 1: dS η,(k) i,t =S η,(k) i,t σ η i (t) | dW t + kS η,(k-1) i,t σ i (t) -σ(t) | dW t , S η,(k) i,0 = 0
One solves this linear system using [START_REF] Protter | Stochastic integration and differential equations[END_REF] and a readily induction leads to, for any k ≥ 1:

S η,(k) i,t =k!S η i,t t 0 t k-1 0 . . . t 1 0 dZ η i,t 0 . . . dZ η i,t k-2 dZ η i,t k-1 (9)
with Z η i,. the Gaussian process defined by:

dZ η i,t = σ i (t) -σ(t) | dW t -σ η i (t)dt . ( 10
)
When considering the above derivatives and processes at η = 0, we use the notations S (k) i,t := S 0,(k) i,t and Z i,t := Z 0 i,t . In addition, S η,(k) t

= d i=1 α i S η,(k) i,t
stands for the derivatives of S η t w.r.t. η. With (9), the processes useful for the next calculations are defined, for any t ∈ [0, T ], by :

           S (k) t = d i=1 α i S (k) i,t , S (k) i,t = k!S P t t 0 t k-1 0 . . . t 1 0 dZ i,t 0 . . . dZ i,t k-2 dZ i,t k-1 , i ∈ {1, . . . , d}, dZ i,t = σ i (t) -σ(t) | dW t -σ(t)dt . (11)
Owing to the identities d i=1 α i σ i = σ and d i=1 α i = 1, observe with (11) that, ∀t ∈ [0, T ] :

S (1) t = d i=1 α i S (1) i,t = d i=1 α i S P t t 0 σ i (s) -σ(s) | dW s -σ(s)ds = 0 (12)

Taylor expansions.

Assume that h ∈ C 3 P (R). To obtain a third order approximation formula, perform Taylor expansions twice : first for the function h at S = S T around S = S P T , second for the interpolated process S η T at η = 1 around η = 0 :

E h(S T ) =E h(S P T ) + E h (1) (S P T )(S T -S P T ) + 1 2 E h (2) (S P T )(S T -S P T ) 2 + . . . =E h(S P T ) + E h (1) (S P T ) S (2)
T 2 + E h (1) (S P T )

S (3) T 3! + 1 2 E h (2) (S P T ) S (2) T 2 2 + Error 3,h , (13)
using the identity [START_REF] Carmona | Generalizing the black-scholes formula to multivariate contingent claims[END_REF]. The explicit calculus of the corrective terms E h (1)

(S P T ) S (2)
T 2 , E h (1) (S P T )

S (3) T 3! and 1 2 h (2) (S P T ) S (2)
T 2 2 is provided in Proposition 8 of subsection 3.1 whereas estimate of Error 3,h is given in subsection 3.2. This leads to the Theorem (stated for only h ∈ H 1 P (R)) given in the next subsection.

Third order weak approximation.

Theorem 5 (Third order weak approximation using the log-normal proxy). Assume (H σ ) and suppose that h ∈ H 1 P (R). Then we have the following weak approximation:

E h(S T ) =E h(S P T ) + Cor 3,h + Error 3,h , (14)
where the corrective term Cor 3,h is given by:

Cor 3,h = 1 2 G h 2 + G h 3 d i=1 α i C 2 i + G h 2 + 3G h 3 + G h 4 1 6 d i=1 α i C 3 i + 1 2 i, j∈{1,...,d} α i α j C i C j C i, j + 1 4 G h 2 i, j∈{1,...,d} α i α j C 2 i, j + 1 8 G h 2 + 15G h 3 + 25G h 4 + 10G h 5 + G h 6 d i=1 α i C 2 i 2 , ( 15
)
with the Greeks G h defined in [START_REF] Bayer | Small-time asymptotics for the at-the-money implied volatility in a multi-dimensional local volatility model[END_REF] and with the coefficients:

C i = T 0 σ i (t) -σ(t) | σ(t) dt, C i, j = T 0 σ i (t) -σ(t) | σ j (t) -σ(t) dt. ( 16
)
The error term is estimated as follows:

|Error 3,h | ≤ c C h (1) µ[(M σ,. ) 4 ]T 2 . ( 17
)
Corollary 6. If one prefers to restrict to a second order approximation, it simply writes:

E h(S T ) =E h(S P T ) + 1 2 G h 2 + G h 3 d i=1 α i C 2 i + O C h (1) µ[(M σ,. ) 3 ]T 3 2 (18)
Proof. It is sufficient to show that additional corrective terms of the expansion ( 14) are of order O C h [START_REF] Alexander | Analytic approximations for multi-asset option pricing[END_REF] µ[(M σ,. ) 3 ]T 3 2 using Proposition 3 for the Greeks and regarding the magnitudes of the coefficients C i and C i, j defined in [START_REF] Durrleman | From implied to spot volatilities[END_REF]. We let this verification to the reader.

We make several additional remarks: a) As announced, the expansion involves only a scalar log-normal r.v. with coefficients C i an C i, j independent of the payoff function h.

b) The result states that for a large class of test functions h, the error is of order four w.r.t. the standard deviation M σ √ T and is null if M σ = 0 (meaning that all the r.v. S i,T are identical) or again if C h (1) = 0 (meaning that h is constant). On the other hand [START_REF] Gobet | Weak approximation of averaged diffusion processes[END_REF] provides an error estimate only in terms of |σ| ∞ √ T (magnitude of the volatilities), what does not take into account the dispersion of the volatilities.

c) The call/put parity relationship is preserved within these approximations. Indeed, remark that the above expansion formula is exact for the payoff function h(S ) = (S -K) as E h(S T ) = E h(S P T ) = 1 -K and as the Greeks, all equal to or greater than order 2, vanish.

d) Notice that if, all the r.v. S i,T are i.i.d. with the common volatility coefficient σ and if α i := 1/d for any i ∈ {1, . . . , d}, then C i = 0 2 for any i ∈ {1, . . . , d}. Then the second order expansion [START_REF] Ju | Pricing asian and basket options via taylor expansion[END_REF] reduces to the main term. To achieve a better accuracy, it is thus necessary to use the third order formula whose the corrective term reduces to Cor 3,h = 1 4 G h 2 i, j∈{1,...,d} α i α j C 2 i, j = 1 4 G h 2 (d -1)(V P T ) 2 . 3 e) As the generic constants do not depend on d, one could under technical assumptions pass to the limit when the discrete probability measure µ converges to some probability measure (see the arguments in [START_REF] Gobet | Weak approximation of averaged diffusion processes[END_REF]Section 2.5]). These investigations with applications to asian options are left for further research.

Proof of Theorem 5.

Computation of the expansion coefficients.

3.1.1. Notations. In all the following, ϕ is a smooth function with bounded derivatives. 2 In this case σ i = (0, . . . , 0, σ i , 0, . . . , 0) * and σ = σ/d(1, . . . , 1) * . Thus it comes

C i = T 0 σ i (t) | σ(t) -σ(t) | σ(t) dt = T 0 σ 2 (t)/d -σ 2 (t)/d dt = 0. 3 In this case, a similar calculus gives that C i, j = T 0 σ i (t) -σ(t) | σ j (t) -σ(t) dt equals C 1,1 := (d -1)/d T 0 σ 2 (t)dt = (d -1)V P T if i = j or C 1,2 := -1/d T 0 σ 2 (t)dt = -V P T if i j. Then it comes i, j∈{1,...,d} α i α j C 2 i, j = 1/d 2 dC 2 1,1 + d(d -1)C 2 1,2 = (d -1)(V P T ) 2 .
£ Integral operator. We adopt the following notation for any n ≥ 1 and any l 1 , . . . , l n measurable and bounded functions:

ω(l 1 , . . . , l n ) T 0 := T 0 l 1 (t 1 ) T t 1 l 2 (t 2 ) . . . T t n-1 l n (t n )dt n . . . dt 2 dt 1 . (19)
£ Differential operator. We introduce a differential operator L intensively used in the following:

Lϕ : S → S ϕ (1) (S ). ( 20
)
Its compound of order n is denoted by

L n := L • • • • • L n times
. By convention L 0 is the identity operator.

Integration by parts formulas.

First we provide integration by parts formulas for S P T in the following Lemma which proof is postponed to Appendix A.1.

Lemma 7.

Let N ≥ 1 be fixed, and consider for j = 1, . . . , N, measurable and bounded deterministic functions t → l j (t) taking values in R and t → L j (t) taking values in R q . Then, introduce for i = 1, . . . , N, the following 2N uni-dimensional processes defined by

dY 0 i (t) = l i (t)dt, dY 1 i (t) = L i (t) | dW t .
Using the notation N := N k=1 I k for any (I 1 , . . . , I N ) ∈ {0, 1} N , the following identity holds:

E ϕ S P T T 0 t N 0 . . . t 2 0 dY I 1 1 (t 1 ) . . . dY I N-1 N-1 (t N-1 )dY I N N (t N ) = ω(Λ 1 , . . . , Λ N ) T 0 E L N ϕ S P T , (21)
where

Λ k (t) := l k (t) if I k = 0, σ(t) | L k (t) if I k = 1.
and where the integral operator ω and the differential operator L are defined respectively in [START_REF] Kunita | Stochastic flows and stochastic differential equations[END_REF] and [START_REF] Landon | Stratégies de couverture presque optimale : théorie et applications[END_REF]. We make explicit E L k ϕ S P T up to k = 6:

E L k ϕ S P T =                          G ϕ 1 for k = 1, G ϕ 1 + G ϕ 2 for k = 2, G ϕ 1 + 3G ϕ 2 + G ϕ 3 for k = 3, G ϕ 1 + 7G ϕ 2 + 6G ϕ 3 + G ϕ 4 for k = 4, G ϕ 1 + 15G ϕ 2 + 25G ϕ 3 + 10G ϕ 4 + G ϕ 5 for k = 5, G ϕ 1 + 31G ϕ 2 + 90G ϕ 3 + 65G ϕ 4 + 15G ϕ 5 + G ϕ 6 for k = 6,
where

G ϕ i = E ϕ (i) S P T S P T i = ∂ i ε i E ϕ (1 + ε)S P T ε=0 = ∂ i S i 0 E ϕ(S P T )
owing to the regularity of ϕ.

The explicit calculus of the corrective terms is given in the next Proposition proven in Appendix A.2.

Proposition 8. One has for the coefficients C i and C i, j defined in Theorem 5:

1 2 E ϕ (1) (S P T ) S (2) T 2 = 1 2 G ϕ 2 + G ϕ 3 d i=1 α i C 2 i , (22) 1 6
E[ϕ (1) (S P T )S (3) T ] =

1 6 G ϕ 2 + 3G ϕ 3 + G ϕ 4 d i=1 α i C 3 i , (23) 1 2 E ϕ (2) (S P T ) S (2) T 2 2 = 1 2 G ϕ 2 + 3G ϕ 3 + G ϕ 4 i, j∈{1,...,d} α i α j C i C j C i, j + 1 4 G ϕ 2 i, j∈{1,...,d} α i α j C 2 i, j + 1 8 G ϕ 2 + 15G ϕ 3 + 25G ϕ 4 + 10G ϕ 5 + G ϕ 6 d i=1 α i C 2 i 2 . ( 24
)
Remind that the Greeks in the r.h.s. are well defined even if ϕ is not smooth.

Error analysis.

Representation of Error 3,h for smooth h and outline of the proof.

Definition 9. Assume (H σ ). We introduce for any k ∈ N, the S-residual processes defined by: R S ,k = S -

k j=0 S ( j) j! = 1 0 S η,(k+1) (1 -η) k k! dη. ( 25
)
where by convention, S (0) = S P . Owing to [START_REF] Carmona | Generalizing the black-scholes formula to multivariate contingent claims[END_REF], one has the useful identity R S ,0 = R S ,1 .

Assume that h ∈ C 3 P (R) to obtain using ( 13) the following representation for Error 3,h :

Error 3,h =E h (1) (S P T )R S ,3 T + 1 2 E h (2) (S P T )R S ,2 T (R S ,2 T + S (2) T ) +E R S ,1 T 3 1 0 h (3) S P T + ηR S ,1 T 1 -η 2 2 dη . ( 26
)
One proves the estimate [START_REF] Gobet | Weak approximation of averaged diffusion processes[END_REF] for the above error term in three steps:

1. L p norm estimates of the interpolated process S η and of its derivatives; 2. Small Gaussian noise perturbation to smooth the function h; 3. Careful use of Malliavin integration by parts formulas to achieve the proof.

Approximation of S and related error estimates.

Lemma 10. Assume (H σ ). We have the following estimates ∀p ≥ 1, and for any k ≥ 1:

sup t∈[0,T ],η∈[0,1] ||S η,(k) i,t || p ≤ c M σ,i √ T k , ∀i ∈ {1, . . . , d}, (27) sup t∈[0,T ],η∈[0,1] ||S η,(k) t || p ≤ c µ[(M σ,. ) k ]T k 2 , (28) sup t∈[0,T ] ||R S ,k t || p ≤ c µ[(M σ,. ) k+1 ]T k+1 2 , (29)
Proof. W.l.o.g. assume p ≥ 2. For (27), starting from [START_REF] Berestycki | Asymptotics and calibration of local volatility models[END_REF], one shows by induction that we have the following bound for any k ≥ 1:

sup t∈[0,T ],η∈[0,1] t 0 t k-1 0 . . . t 1 0 dZ η i,t 0 . . . dZ η i,t k-2 dZ η i,t k-1 p ≤ c M σ,i √ T k
with the Gaussian process dZ

η i,t = σ i (t) -σ(t) | dW t -σ η i (t)dt .
The result stands for k = 1 owing to standard computations involving the Burkholder-Davis-Gundy and Hölder inequalities:

t 0 dZ η i,s p p ≤ c t p 2 M p σ,i + t p M p σ,i |σ| p ∞ ≤ c t p 2 M p σ,i (30)
Then if the results stands for a certain rank k ≥ 1, the same calculations give the estimate:

t 0 t k 0 . . . t 1 0 dZ η i,t 0 . . . dZ η i,t k-1 dZ η i,t k p p ≤ c t p 2 -1 M p σ,i + t p-1 M p σ,i |σ| p ∞ t 0 t k 0 . . . t 1 0 dZ η i,t 0 . . . dZ η i,t k-1 p p dt k ≤ c M σ,i √ T p(k+1) .
We are done. The second estimate (28) is handled with the Minkowski inequality and ( 27). The proof of ( 29) is straightforward using the representation ( 25) and the previous estimate (28).

3.2.3.

Regularization of h with a small noise perturbation. For smooth function h, estimate of Error 3,h is straightforward using its representation (26) and Lemma 10. To account for non-smooth functions and to overcome some degeneracy in the Malliavin sense, we suitably regularize the function h using small noise perturbation. This standard scheme of proof has been successfully employed in [START_REF] Gobet | Weak approximation of averaged diffusion processes[END_REF]- [START_REF] Bompis | Stochastic Approximation Finite Element method for analytical approximation of multidimensional diffusion[END_REF] and we adapt the arguments for GBMs processes instead of Gaussian processes. Let W ⊥ be an extra independent scalar BM and consider the C ∞ P (R)-function:

h δ (S ) = E h S × exp(δW ⊥ T ) = E h δ/ √ 2 S × exp(δW ⊥ T/2 ) , (31)
with the small parameter δ assumed positive:

δ = µ[(M σ,. ) 4 ]T 3 2 > 0. (32)
Replacing h by h δ in our expansion analysis induces extra errors quantified below.

Lemma 11. Assume (H σ ) and suppose that h ∈ H 1 P (R). Then, for any n ∈ N we have the estimates:

sup η∈[0,1] E h δ (S η T ) -E h(S η T ) ≤ c C h (1) µ[(M σ,. ) 4 ]T 2 , G h δ n -G h n ≤ c C h (1) µ[(M σ,. ) 4 ]T 2 V P T -n 2 ,
with h δ defined in (31) and G h δ n the Greeks for h δ .

Proof. For the first estimate, write:

h δ (S η T ) -h(S η T ) = E S (e δW ⊥ T -1) 1 0 h (1) S 1 + λ(e δW ⊥ T -1) dλ S =S η T .
and conclude with the estimate

E e δW ⊥ T -1 = E δW ⊥ T 1 0 e λδW ⊥ T dλ = O(δ √ T ).
For the second estimate, denoting by D(y) := expy 2 /(2V P T ) / 2πV P T for any y ∈ R the Gaussian density, write for ε > -1:

G h δ n -G h n = R E h e y-1 2 V P T +δW ⊥ T -h e y-1 2 V P T ∂ i ε i D y -ln(1 + ε) ε=0 dy
Complete the proof by combining the first estimate with standard upper bounds for derivatives of D.

Regarding the magnitudes of the coefficients C i and C i, j defined in ( 16), Lemma 11 readily gives the estimate:

Error 3,h = E h(S T ) -E h(S P T ) -Cor 3,h ≤ Error 3,h δ + O C h (1) µ[(M σ,. ) 4 ]T 2 .
Hence, proving the error estimate ( 17) is reduced to prove the following Proposition:

Proposition 12. Assume (H σ ) and suppose that h ∈ H 1 P (R). We have the estimate:

Error 3,h δ = O C h (1) µ[(M σ,. ) 4 ]T 2 .
This is performed in Appendix A.3 using Malliavin integration by parts formulas.

Application to the Basket Call options.

In all the following, we consider basket call options with payoff h(S ) = (S -K) +4 , for a strike K > 0. Hence h ∈ H 1 P (R) and Theorem 5 is applicable. In order to obtain more tractable and hopefully more accurate formulas, we also provide approximations of the implied volatility.

Notations.

£ Call options. We denote by Call(T, K) the price at time 0 of a basket call option with maturity T and strike K, written on the asset

S = d i=1 α i S i that is Call(T, K) = E[(S T -K) + ].
As usual, ATM (At The Money) Call refers to K ≈ S 0 = 1, ITM (In The Money) to K 1, OTM (Out The Money) to K 1. £ Black-Scholes Call price function. For the sake of completeness, we give the Black-Scholes Call price function depending on the spot S 0 , the total variance V and the strike K :

Call BS (S 0 , V, K) = S 0 N d 1 (S 0 , V, K) -KN d 2 (S 0 , V, k) (33)
where:

N(x) = x -∞ N (y)dy, N (y) = e -y 2 /2 √ 2π , d 1 (S 0 , V, K) = ln S 0 /K √ V + 1 2 √ V, d 2 (S 0 , V, K) =d 1 (S 0 , V, K) - √ V. One has E[(S P T -K) + ] = Call BS (1, V P T , K) which equals Call(T, K) = E[(S T -K) + ] when M σ = 0. When unambiguous, we also write d 1 = d 1 (1, V P T , K) and d 2 = d 2 (1, V P T , K)
and use the notation k := ln(K) for the log-moneyness.

£ Implied Black-Scholes volatility. For (T, K) given, the implied Black-Scholes volatility of a basket call option price Call(T, K) is the unique non-negative volatility parameter σ I (T, K) such that:

(34) Call BS 1, σ I (T, K) 2 T, K = Call(T, K).
£Quadratic mean of the volatility on [0, T ].

Finally we introduce σ(0, T ) =

1 T T 0 | σ(t)| 2 dt = 1 T V P
T the quadratic mean of the proxy volatility parameter on [0, T ].

Price approximation formulas.

Price approximation formulas for basket call options are easily obtained using Theorem 5, Corollary 6 and the closed forms of Greeks (see Appendix B.1):

Theorem 13 (Second and third order basket call option price approximations). Assume (H σ ). Then we have the basket call price approximations:

Call(T, K) =Call BS (1, V P T , K) - 1 2 N (d 1 ) H 1 (d 1 ) V P T d i=1 α i C 2 i + O µ[(M σ,. ) 3 ]T 3 2 , (35) =Call BS (1, V P T , K) + N (d 1 ) V P T - 1 2
H 1 (d 1 )

V P T d i=1 α i C 2 i + H 2 (d 1 ) V P T 1 6 d i=1 α i C 3 i + 1 2 i, j∈{1,...,d} α i α j C i C j C i, j + 1 4 i, j∈{1,...,d} α i α j C 2 i, j + 1 8
H 4 (d 1 )

V P T 2 d i=1 α i C 2 i 2 + O µ[(M σ,. ) 4 ]T 2 , ( 36
)
with C i and C i, j defined in Theorem 5 and with the Hermite polynomials H i defined in Appendix B.1.

Expansion formulas for the implied volatility.

To obtain approximations of the implied volatility, write expansions (35)-(36) in terms of the sensitivities Vega BS and Vomma BS (defined in Appendix B.1) w.r.t. the volatility to obtain the approximations:

Call(T, K) =Call BS (1, V P T , K) + Vega BS 1, σ(0, T ) 2 T, K Σ I 2 (T, K) + O µ[(M σ,. ) 3 ]T 3 2 =Call BS (1, V P T , K) + Vega BS 1, σ(0, T ) 2 T, K [Σ I 2 (T, K) + Σ I 3 (T, K)] + 1 2 Vomma BS Σ I 2 (T, K) 2 + O µ[(M σ,. ) 4 ]T 2
where Σ I 2 (T, K) and Σ I 3 (T, K) are defined by:

Σ I 2 (T, K) = - H 1 (d 1 ) 2 σ2 (0, T )T 3 2 d i=1 α i C 2 i , (37)
Σ I 3 (T, K) = H 2 (d 1 ) σ3 (0, T )T 2 1 6 d i=1 α i C 3 i + 1 2 i, j∈{1,...,d} α i α j C i C j C i, j + 1 σ(0, T )T 1 4 i, j∈{1,...,d} α i α j C 2 i, j + 1 8 (d 3 1 σ(0, T ) √ T -6d 2 1 + 3) σ5 (0, T )T 3 d i=1 α i C 2 i 2 (38)
These formulas read as expansions of the implied volatility. The next Theorem states without proof (derivation of error estimates is quite standard, we refer to [10, Section 3.4] for full details) that the approximated implied volatility approximations write as polynomials of degrees 1 or 3 w.r.t. the log-moneyness k := ln(K).

Theorem 14 (Second and third order implied volatility approximations). Assume (H σ ). We have the next implied volatility approximations:

Call(T, K) =Call BS ( σ I 2 (T, K) 2 T, K) + O µ[(M σ,. ) 3 ]T 3 2 , (39) Call(T, K) =Call BS ( σ I 3 (T, K) 2 T, K) + O µ[(M σ,. ) 4 ]T 2 . ( 40
)
with σ I 2 (T, K) and σ I 3 (T, K) such that:

σ I 2 (T, K) = σ(0, T ) + Σ I 2 (T, K), σ I 3 (T, K) = σ I 2 (T, K) + Σ I 3 (T, K)
, and Σ I 2 (T, K) and Σ I 3 (T, K) defined in (37)-( 38). The next Corollary analyses the behaviour of the implied volatility ATM for short maturity by formally differentiating5 σ I 3 (T, K) w.r.t. k := ln(K) at (T, k) = (0, 0). Corollary 15 (Shape of the implied volatility ATM for short maturity). Assume that σ(t) is continuous at t = 0. Then the level, the slope and the curvature ATM for short maturity are approximated by the next formulas:

σ I (T, k)| (T =0,k=0) ≈ σ(0, 0) = | σ(0)|, ∂ k σ I (T, k)| (T =0,k=0) ≈ 1 2 σ(0, 0) 3 d i=1 α i C i (0) 2 , ∂ 2 k 2 σ I (T, k)| (T =0,k=0) ≈ 1 σ(0, 0) 5 1 3 d i=1 α i C i (0) 3 + i, j∈{1,...,d} α i α j C i (0)C j (0)C i, j (0) - 3 2 σ(0, 0) 7 d i=1 α i C i (0) 2 2 . with C i (0) = σ i (0) -σ(0) | σ(0) and C i, j (0) = σ i (0) -σ(0) | σ j (0) -σ(0) for any i, j ∈ {1, . . . , d}.
Our level and skew exactly coincide with those of [START_REF] Piterbarg | Markovian Projection Method for Volatility Calibration[END_REF], [START_REF] Durrleman | From implied to spot volatilities[END_REF] and [START_REF] Alos | On the short-maturity behaviour of the implied volatility skew for random strike options and applications to option pricing approximation[END_REF] (see Appendix B.3).

Numerical experiments.

In the context of basket call options, we numerically illustrate the accuracy of the second and third order approximation formulas on both prices and implied volatilities provided in Theorems 13 and 14. First, we consider two 10-dimensional examples : firstly an i.i.d. case and secondly a correlated case. Then we consider a 2-dimensional toy model and perform an impact analysis of the correlation. Additional tests for different volatilities and maturities are left for further works. In all the tests, we approximate the implied volatility for various strikes denoted by K, chosen to approximately equal e q| σ| √ T where q takes the value of various quantiles of the standard Gaussian law (1%-5%-10%-20%-30%-40%-50%-60%-70%-80%-90%-95%-99%) to consider far ITM and far OTM options. As a benchmark, we use Monte Carlo simulations (denoted by (MC)) with 10 9 simulations and the Proxy as control variate to speed the convergence. This allows a 95%-confidence interval width of one bp 6 or less for all the following Monte Carlo estimates. Then we compute the approximative implied volatilities using our second and third order approximation formulas on both prices and implied volatilities (denoted respectively by Vol(App2), App2Vol, Vol(App3) and App3Vol) and the following benchmark formulas of the literature (see Appendix B.2): D'Aspremont [START_REF]Interest rate model calibration using semidefinite programming[END_REF] (order 2), Carmona and Durrleman [START_REF] Carmona | Generalizing the black-scholes formula to multivariate contingent claims[END_REF] (a priori order 2), Landon [START_REF] Landon | Stratégies de couverture presque optimale : théorie et applications[END_REF] (order 2), Gobet-Miri [START_REF] Gobet | Weak approximation of averaged diffusion processes[END_REF] (order 2), Ju [START_REF] Ju | Pricing asian and basket options via taylor expansion[END_REF] (order 3) and Shiraya-Takahashi [START_REF] Shiraya | Pricing multi-asset cross-currency options[END_REF] (order 3) formulas denoted respectively by (DAF), (CDF) inf , (CDF) sup , (LF), (GMF), (JUF) and (STF). In addition we have proposed and tested a common modification of the D'Aspremont and Landon formulas denoted by (MDALF) (see equation ( 67)). All the following computations are performed using C++ on a Intel(R) Core(TM) i5 CPU@2.40GHz with 4 GB of ram.

I.i.d. case in dimension 10.

In the first test, we take for any i ∈ {1, . . . , d} the parameters::

d =10, S 0,i =1, T =1Y, α i =10%, Vols =30% × I 10 ,
where Vols is the matrix which contains (in row) all the volatility vectors of the assets and I 10 the identity matrix in dimension 10. The values of the parameters allowing to quantify the error are the following: | σ| = 30%/ √ 10 ≈ 9.49% and M σ = µ[M σ,. ] = 27%. Although the proxy volatility is quite small due to the diversification, the dispersion parameters M σ and µ[M σ,. ] are quite important. The execution time of (MC) is close to 12 hours while the execution time of the various analytical approximations is almost instantaneous: less than 0.1s. Results are given in Table 1. First we represent the results of order 2 formulas and then those of order 3 in order to make easier comparisons. As all the second order formulas, excepting (CDF) sup and (GMF), lead exactly to the same results 7 with a uniform approximated implied volatility equalling 9.49%, we only indicate results for Vol(App2), (CDF) sup and (GMF).

Regarding the results, first remark that the true implied volatility estimated by (MC) has small variations (5 bps of amplitude). Error for Vol(App2) is around 20 bps what is quite important and it is necessary to use a third order formula to improve the accuracy. (GMF) is quite inaccurate (more Estimates of the implied volatility (%) in the 10-dimensional i.i.d. case.

than 100 bps of error) and approximated prices are sometimes outside arbitrage bounds (indicated by ND in Table 1). The upper bound (CDF) sup equals the volatility of the assets and is too large. Then for the third order, our formulas are very closed to (MC) with an error less than 1 bp for all the strikes excepting the largest strike for which the error is 4 bps. (JUF) [START_REF] Ju | Pricing asian and basket options via taylor expansion[END_REF] is the most accurate method with errors less than 1 bp for all strikes, but the formula remains a little bit more complicated than ours.

Notice that the proxy volatility used in [START_REF] Ju | Pricing asian and basket options via taylor expansion[END_REF] (see (69)) equals 9.68%, i.e. the ATM implied volatility, whereas | σ| = 9.49%. (STF) is very accurate ATM, but less accurate ITM and OTM (errors up to 16 bps). Finally comparing Vol(App3) and App3Vol, notice that results are very close. For its simplicity, it seems better to use App3Vol.

As a conclusion the accuracy of our formulas is very satisfying and our third order formulas surpass in accuracy all the benchmark approximations excepting (JUF) for high strikes which is slightly more accurate.

Correlated case in dimension 10

. We now consider assets having initial values different from 1 and we use the normalisation described in Remark 1. We use a volatility matrix generated with a uniform random sampling 8 . The parameters of the model are the following (with α 2 the weights after normalisation) :

d =10, S 0
=                                               34
                                              × 1% d i=1 α 2,i |σ i | ≈ 32%, | σ| ≈ 19.45%, M σ ≈ 26.78%, µ[M σ,. ] ≈ 21.79%
The mean of the volatilities of the assets is close to 30%, the volatility of the proxy is quite important, around 20%, and the dispersion parameters are large (M σ close to 27% and µ[M σ,. ] close to 22%). We give the results with this new set of parameters in

Table 2

Estimates of the implied volatility (%) in the 10-dimensional correlated case. The true implied volatility estimated by (MC) has now large variations (from 17.32% to 21.31%). The best price second order formula is the lower bound (CDF) inf : errors are roughly of 15 bps for all strikes. Regarding App2Vol, errors are close to 20 bps ATM what is quite important but the fitting is curiously quite good ITM and OTM (7 bps and 2 bps for the extreme strikes). The upper bound (CDF) sup is close to the mean of the volatility of the assets and is too large. Notice that as expected, Vol(App2) and (MDALF) give very close results, mostly ATM : 20 bps of error ATM and even less accurate results for the extreme strikes (more than 30 bps). Observe however that Vol(App2) gives slightly better results ITM and OTM. Regarding (LF), results are excellent ATM but bad ITM and OTM. (GMF) leads to clearly less accurate results, ATM as well as OTM and ITM. Finally, for the third order, observe that results are excellent OTM and very close whatever is the considered approximation. Results for (STF) are bad for the largest strike and middling for (JUF) for the smallest strike (error of 27 bps) while remaining satisfying for Vol(App3) and App3Vol (5 and 9 bps OTM and 15 and 12 bps ITM). Accuracy is very good for Vol(App3) and App3Vol ITM and OTM, slightly better for Vol(App3) OTM and slightly better for App3Vol ITM. Notice that results are biased for (JUF) ATM (5 bps). The proxy volatility of (JUF) (see(69)) is around 19.83% whereas | σ| ≈ 19.45% which is closer to the ATM implied volatility 19.57%. Here again, the accuracy of our third order formulas is excellent, showing an improving accuracy in comparison to the literature formulas.

Impact of the correlation in dimension 2.

Here we consider a toy model in dimension 2 with :

d =2, S 0,1 =S 0,2 = 1, T =1Y, α 1 =α 2 = 50%, Vols = 30% 0 ρ × 30% 1 -ρ 2 × 30% ,
where the correlation ρ takes the values {-90%, -75%, -50%, -25%, -10%, 0%, 10%, 25%, 50%, 75%, 90%}. For the second order, we only indicate the results for Vol(App2) and (GMF) because results are the same for all the methods excepting (GMF). Whatever is the correlation, (CDF) sup is flat and equals 30% what is a strong overestimation except for correlation close to 100%. For the third order, Vol(App3) and App3Vol give same results for ρ ≥ 25%. Results are given in Tables 3 to 13. We also indicate the values of the parameters | σ|, M σ, µ[M σ,. ] and σ (the proxy volatility σ used in [START_REF] Ju | Pricing asian and basket options via taylor expansion[END_REF], see (69)).

First observe that the more the correlation is close to -100% (anti-correlated case):

• The more the volatility of the proxy | σ| (from 6.71% to 29.24%) has small magnitude.

• The more the true implied volatility has variations: one observes a positive skew whatever is the correlation. For instance the implied volatility varies from 7.66% to 10.41% for ρ = -90%, from 21.39% to 21.53% for ρ = 0% and is flat, equal to 29.24% for ρ = 90%. • The more the dispersion parameters M σ and µ[M σ,. ] (which are always equal) are important (from 28.50% to 6.54%). • The less the approximation formulas are accurate. Whatever is the correlation, Vol(App2) is flat, equals | σ| and always underestimates the (MC) volatility. Results are quite bad for negative correlations (from 100 bps to 370 bps for ρ = -90%, and from 20 bps to more than 30 bps for ρ = -10%). Errors are around 10 bps for ρ = 25%, 5 bps for ρ = 50% and 0 bp for ρ = 90% (in this case, all the estimations coincide and are flat). (GMF) is clearly less accurate with estimations outside the arbitrage bounds for very negative correlation values. Regarding the third order, App3Vol gives flat estimation of the volatility whatever is the correlation. Results are bad for very negative correlations (from 170 bps to 100 bps for ρ = -90% and from 80 bps to 90 bps for ρ = -75%), are average for negative correlations (from 25 bps to 45 bps for ρ = -50% and from 15 bps to 20 bps for ρ = -25%) and are acceptable when approaching ρ = 0% (around 10 bps for ρ = -10% and around 5 bps for ρ = 0%). For positive correlations, errors are always less than 5 bps and results are excellent ATM. Vol(App3) is not flat with an increasing behaviour ITM and a decreasing behaviour OTM. The formula coincides with App3Vol ATM, is globally less accurate OTM and more accurate ITM with errors from 60 bps to 210 bps for ρ = -90%. The approximation starts to be flat from ρ = -25%. For its robustness, we prefer use App3Vol. (STF) is also not flat with an increasing behaviour ITM and a decreasing behaviour OTM. Results are similar to App3Vol ATM but globally quite bad ITM and OTM whatever is the correlation (from 170 bps to 250 bps for ρ = -90%, from 30 bps to 70 bps for ρ = 0% and from 140 bps to 150 bps for ρ = 75%). This is due to the fact that even for the very correlated case (ρ close to 100%), the underlying model of the method (STF) is Gaussian whereas the true model is almost log-normal, what induces a high sensitivity on the wings. Finally regarding (JUF), one remarks that results are globally less bad than all the other formulas for very negative correlation (ND for ρ = -90%), excepting for the smallest strike where our third order approximations give better results. Results are satisfying from the correlation -25% and very good from the correlation -10% with smaller errors than our formulas. In particular the fitting ATM is excellent. This is due to the fact that the proxy volatility σ used in [START_REF] Ju | Pricing asian and basket options via taylor expansion[END_REF] defined in (69) is very close to the ATM implied volatility, whatever is the correlation (excepting for the correlation -90% with an error of around 50 bps). 20.12 20.12 20.12 20.12 20.12 20.12 20.12 20.12 20.12 20.12 20.12 20.12 20.12 (GMF)

Table 13

Estimates of the implied volatility (%) in the 2-dimensional case with ρ = 90%. On has | σ| ≈ 29.24%, M σ = µ[M σ,. ] ≈ 6.54% and σ ≈ 29.24%.

Lemma 16. For any c : [0, T ] → R q square integrable and predictable process, one has the following identity:

E ϕ S P T T 0 c(t) | dW t = E Lϕ S P T T 0 σ(t) | c(t) dt ,
where the differential operator L is defined in [START_REF] Landon | Stratégies de couverture presque optimale : théorie et applications[END_REF]. The case (N = 1, I 1 = 0) is obvious whereas the case (N = 1, I 1 = 1) exactly corresponds to Lemma 16. Now suppose that the formula (21) holds for N ≥ 2. Then apply Lemma 16 if I N+1 = 1 to get:

E ϕ S P T T 0 t N+1 0 . . . t 2 0 dY I 1 1 (t 1 ) . . . dY I N N (t N )dY I N+1 N+1 (t N+1 ) =E L I N+1 ϕ S P T T 0 Λ N+1 (t N+1 ) t N+1 0 . . . t 2 0 dY I 1 1 (t 1 ) . . . dY I N N (t N )dt N+1 =E L I N+1 ϕ S P T T 0 T t N Λ N+1 (s)ds dY I N N (t N ) t N 0 . . . t 2 0 dY I 1 1 (t 1 ) . . . dY I N-1 N-1 (t N-1 ) ,
where we applied at the last equality the identity

T 0 f (t)Z t dt = T 0
T t f (s)ds dZ t for any continuous semi-martingale Z starting from 0 and any measurable, bounded and deterministic function f (apply the Itô's Lemma to the product Z t T t f (s)ds). We conclude without difficulty.

A.2. Proof of Proposition 8. For [START_REF] Nualart | Malliavin calculus and related topics[END_REF], write E ϕ (1) (S P T ) (1) (S P T )S (2) i,T ] and use [START_REF] Bompis | Stochastic Approximation Finite Element method for analytical approximation of multidimensional diffusion[END_REF] to obtain the following calculus:

S (2) T 2 = d i=1 α i E[ϕ
1 2 E[ϕ (1) (S P T )S (2) i,T ] =E ϕ (1) (S P T )S P T T 0 t 0 σ i (s) -σ(s) | σ(s)ds σ i (t) -σ(t) | σ(t)dt . -E ϕ (1) (S P T )S P T T 0 t 0 σ i (s) -σ(s) | σ(s)ds σ i (t) -σ(t) | dW t -E ϕ (1) (S P T )S P T T 0 t 0 σ i (s) -σ(s) | dW(s) σ i (t) -σ(t) | σ(t)dt +E ϕ (1) (S P T )S P T T 0 t 0 σ i (s) -σ(s) | dW(s) σ i (t) -σ(t) | dW t
Then write ϕ (1) (S P T )S P T = Lϕ(S P T ) and apply Lemma 7 to get the identity:

1 2 E[ϕ (1) (S P T )S (2) i,T ] =E {L -2L 2 + L 3 }ϕ S P T ω σ i -σ | σ , σ i -σ | σ T 0 = G ϕ 2 + G ϕ 3 ω σ i -σ | σ , σ i -σ | σ T 0 .
Conclude with the identity

ω σ i -σ | σ , σ i -σ | σ T 0 1 2 T 0 σ i (t) -σ(t) | σ(t) dt 2 = 1 2 C 2 i .
Then for [START_REF] Piterbarg | Markovian Projection Method for Volatility Calibration[END_REF], write E[ϕ (1) (S P T )S (3) T ] = d i=1 α i E[ϕ (1) (S P T )S (3) i,T ] and use [START_REF] Bompis | Stochastic Approximation Finite Element method for analytical approximation of multidimensional diffusion[END_REF] to obtain for S (3) i,T :

1 6 S (3) i,T = S P T j 1 , j 2 , j 3 ∈{0,1} T 0 t 0 s 0 σ i (u) -σ(u) | dX j 1,u σ i (s) -σ(s) | dX j 2,s σ i (t) -σ(t) | dX j 3,t
with (X j,. ) j∈{0,1} defined by: dX j,t = dW t for j = 0, σ(t)dt for j = 1. (41) Then use again Lemma 7 to obtain the identity:

1 6 E ϕ (1) (S P T )S (3) i,T =E {-L + 3L 2 -3L 3 + L 4 }ϕ S P T ω σ i -σ | σ , σ i -σ | σ , σ i -σ | σ T 0 = 1 6 G ϕ 2 + 3G ϕ 3 + G ϕ 4 C 3 i ,
and achieve the proof. Finally for ( 24), one has using [START_REF] Bompis | Stochastic Approximation Finite Element method for analytical approximation of multidimensional diffusion[END_REF] the following result:

1 2 E ϕ (2) (S P T ) S (2) T 2 2 = 1 8 d i, j∈{1,...,d} α i α j E[ϕ (2) (S P T )S (2) i,T S (2) j,T ] (42)
Then [START_REF] Bompis | Stochastic Approximation Finite Element method for analytical approximation of multidimensional diffusion[END_REF] leads to the identity:

1 4 S (2) i,T S (2) j,T = S P T 2 T 0 t 0 dZ i,s dZ i,t × T 0 t 0 dZ j,s dZ j,t (43)
with dZ k,t defined in [START_REF] Bompis | Stochastic Approximation Finite Element method for analytical approximation of multidimensional diffusion[END_REF]. We decompose the last product of integrals using the Itô's Lemma:

T 0 t 0 dZ i,s dZ i,t × T 0 t 0 dZ j,s dZ j,t = T 0 t 0 s 0 dZ j,u dZ j,s × t 0 dZ i,s dZ i,t + T 0 t 0 s 0 dZ i,u dZ i,s × t 0 dZ j,s dZ j,t + T 0 t 0 dZ i,s × t 0 dZ j,s σ i (t) -σ(t) | σ (t) -σ(t) dt (44)
We begin with the last term in the r.h.s. of ( 44). An application of the Itô's Lemma yields for this term:

T 0 t 0 dZ i,s × t 0 dZ j,s σ i (t) -σ(t) | σ j (t) -σ(t) dt = T 0 t 0 s 0 dZ i,s dZ j,u σ i (t) -σ(t) | σ j (t) -σ(t) dt + T 0 t 0 s 0 dZ j,s dZ i,u σ i (t) -σ(t) | σ j (t) -σ(t) dt + T 0 t 0 σ i (s) -σ(s) | σ j (s) -σ(s) ds σ i (t) -σ(t) | σ j (t) -σ(t) dt
Then write ϕ (2) S P T S P T 2 = {L 2 -L}ϕ S P T to obtain the contribution, introducing the notation

ω i 3 ,i 4 ,(i 1 ,i 2 ) = ω σ i 3 -σ | σ , σ i 4 -σ | σ , σ i 1 -σ | σ i 2 -σ T 0 : E -L + 3L 2 -3L 3 + L 4 ϕ S P T × ω i, j,(i, j) + ω j,i,(i, j) + G ϕ 2 ω σ i -σ | σ j -σ , σ i -σ | σ j -σ T 0 = G ϕ 2 + 3G ϕ 3 + G ϕ 4 × ω i, j,(i, j) + ω j,i,(i, j) + G ϕ 2 ω σ i -σ | σ j -σ , σ i -σ | σ j -σ T 0 (45)
Then for the first term in the r.h.s. of (44), apply the Itô's Lemma twice to obtain:

T 0 t 0 s 0 dZ j,u dZ j,s × t 0 dZ i,s dZ i,t = T 0 t 0 s 0 dZ j,u × s 0 dZ i,u dZ j,s + t 0 s 0 u 0 dZ j,r dZ j,u dZ i,s + t 0 s 0 dZ j,u σ i (s) -σ(s) | σ j (s) -σ(s) ds dZ i,t = T 0 t 0 s 0 u 0 dZ i,r dZ j,u dZ j,s + t 0 s 0 u 0 dZ j,r dZ i,u dZ j,s + t 0 s 0 u 0 dZ j,r Z j,u dZ i,s + t 0 s 0 σ i,u -σ(u) | σ j (u) -σ(u) du dZ j,s + t 0 s 0 dZ j,u σ i (s) -σ(s) | σ j (s) -σ(s) ds dZ i,t .
This leads to the contribution, introducing the notations

ω i 1 ,i 2 ,i 3 ,i 4 = ω σ i 1 -σ | σ , σ i 2 -σ | σ , σ i 3 - σ | σ , σ i 4 -σ | σ T 0 , ω (i 1 ,i 2 ),i 3 ,i 4 = ω σ i 1 -σ | σ i 2 -σ , σ i 3 -σ | σ , σ i 4 -σ | σ T 0 and ω i 3 ,(i 1 ,i 2 ),i 4 = ω σ i 3 -σ | σ , σ i 1 -σ | σ i 2 -σ , σ i 4 -σ | σ 0 : E -L + 5L 2 -10L 3 + 10L 4 -5L 5 + L 6 ϕ S P T × ω i, j, j,i + ω j,i, j,i + ω j, j,i,i + E -L + 3L 2 -3L 3 + L 4 ϕ S P T × ω (i, j), j,i + ω j,(i, j),i = G ϕ 2 + 15G ϕ 3 + 25G ϕ 4 + 10G ϕ 5 + G ϕ 6 × ω i, j, j,i + ω j,i, j,i + ω j, j,i,i + G ϕ 2 + 3G ϕ 3 + G ϕ 4 × ω (i, j), j,i + ω j,(i, j),i (46)
Similarly the contribution of the second term in the r.h.s. of (44) is, inverting the indexes i and j:

(47) G ϕ 2 + 15G ϕ 3 + 25G ϕ 4 + 10G ϕ 5 + G ϕ 6 × ω j,i,i, j + ω i, j,i, j + ω i,i, j, j + G ϕ 2 + 3G ϕ 3 + G ϕ 4 × ω ( j,i),i, j + ω i,( j,i), j
Then use the mathematical reductions: ω i, j,(i, j) + ω j,i,(i, j) + ω (i, j), j,i + ω j,(i, j),i + ω ( j,i),i, j + ω i,( j,i), j =C i C j C i, j , ω i, j, j,i + ω j,i, j,i + ω j, j,i,i + ω j,i,i, j + ω i, j,i, j + ω i,i, j, j = 1 4

C 2 i C 2 j ,
and gather all the contributions (45)-( 46)-( 47) of (44) and use ( 43)-(42) to achieve the proof.

A.3. Malliavin integration by parts formula and proof of Proposition 12.

We consider the Malliavin calculus w.r.t the (q + 1)-dimensional BM (W, W ⊥ ), the Malliavin derivatives associated to W (respectively W ⊥ ) being denoted by D (respectively D ⊥ ) whereas D denotes the Malliavin operator w.r.t. the full BM (W, W ⊥ ). We refer to [START_REF] Nualart | Malliavin calculus and related topics[END_REF] for the related theory and notations for the n-th Malliavin derivatives D n , the Sobolev spaces D n,p and the associated norms || • || k,p in which we extend estimates provided in Lemma 10.

Lemma 17. Under (H σ ), for any η ∈ [0, 1], i ∈ {1, . . . , d}, k ≥ 1 and t ∈ [0, T ], S η i,t ∈ D ∞,∞ , S η,(k) i,t ∈ D ∞,∞ and S η,(k) t ∈ D ∞,∞
. In addition we have the following estimates for any p ≥ 1: i,t , a straightforward induction shows that ||D r t 0 . . .

(48)                                                  ||D r S η i,t || p + ||D r S η t || p ≤ c |σ| ∞ , ||D 2 r,s S η i,t || p + ||D 2 r,s S η t || p ≤ c |σ| 2 ∞ , ||D 3 r,s,u S η i,t || p + ||D 3 r,s,u S η t || p ≤ c |σ| 3 ∞ , ||D r S η,(k) i,t || p ≤ c M k σ,i T (k-1)+ 2 , ||D r S η,(k) t || p ≤ c µ[M k σ,. ]T (k-1)+ 2 , ||D 2 r,s S η,(k) i,t || p ≤ c M k σ,i T (k-2)+ 2 , ||D 2 r,s S η,(k) t || p ≤ c µ[M k σ,. ]T (k-2)+ 2 , ||D 3 r,s,u S η,(k) i,t || p ≤ c M k σ,i T (k-3)+ 2 , ||D 3 r,s,u S η,(k) t || p ≤ c µ[M k σ,. ]T (k-3)+ 2 , ||D r R S ,k t || p ≤ c µ[(M σ,. ) k+1 ]T k 2 , ||D 2 r,s R S ,k t || p ≤ c µ[(M σ,. ) k+1 ]T (k-1)+ 2 , ||D 3 r,s,u R S ,k t || p ≤ c µ[(M σ,. ) k+1 ]T (k-2)+ 2 , uniformly in r, s, t, u ∈ [0, T ], i ∈ {1, . . . ,
t 1 0 dZ η i,t 0 . . . dZ η i,t k-1 || p ≤ c M k σ,i T (k-1)+ 2
. Then using the representation (9), it comes ∀r, t ∈ [0, T ], ∀i ∈ {1, . . . , d} and for any k ≥ 1 the next calculus:

D r S η,(k) i,t =k!1 r≤t [S η,(k) i,t σ η i (r) * + S η i,t D r t 0 . . . t 1 0 dZ η i,t 0 . . . dZ η i,t k-1 ].
The proof is easily achieved using the previous estimate and Lemma 10.

We now state the crucial result related to integration by parts formula which is proved in Appendix A.4.

Proposition 18 (Integration by parts formulas). Assume (H σ ) and suppose that h ∈ H 1 P (R). For λ ∈ [0, 1], we define the random variable F λ δ = S P T + λR S ,1

T × exp(δW ⊥ T/2
). Let j ∈ {2, 3}, for any Y ∈ D j-1,∞ , there exist random variables Y λ 2 and Y λ 3 ∈ ∩ p≥1 L p such that for any j ∈ {2, 3}:

E[h ( j) δ/ √ 2 (F λ δ )Y] = E[h (1) δ/ √ 2 (F λ δ )Y λ j ],
where for any p ≥ 1 and any j ∈ {2, 3} the following estimates hold::

sup λ∈[0,1] ||Y λ j || p ≤ c ||Y|| j-1,p+ 1 2 (V P T ) -( j-1) 2 . (49)
To achieve the proof of Proposition 12, consider Error 3,h explicitly written in (26). The first term is easily handled using (29) with k = 3. For the second term of (26), use (31), apply for any i, j ∈ {1, . . . , d} the above Proposition 18 with Y = R S ,2 T (R S ,2 T + S (2) T ) and use Lemmas 10 and 17 to get the following estimate:

E h (2) δ (S P T )R S ,2 T (R S ,2 T + S (2) T ) = E h (2) δ/ √ 2 (F 0 δ )Y = E h (1) δ/ √ 2 (F 0 δ )Y 0 2 =O C h (1) ||Y|| 1,2 (V P T ) -1 2 = O C h (1) µ[(M σ,. ) 4 ]T 2
The last term of (26) is handled similarly; apply the Proposition 18 with Y = R S ,1 T 3 , j = 3 and use estimates of Lemmas 10 and 17. To summarize, the proof of Theorem 5 is now complete provided that we establish Proposition 18, which is done in the following subsection.

A.4. Proof of Proposition 18.

Owing to (H σ ), S P T is a non degenerate random variable, with Malliavin covariance matrix equal to (S P T ) 2 V P T > 0, but F λ := S P T +λR S ,1 T may be degenerate for λ > 0. Fortunately, F λ δ = F λ × e δW ⊥ T/2 as defined in Proposition 18 is in D ∞,∞9 and is non degenerate since its Malliavin covariance matrix γ F λ δ is strictly positive with γ F λ δ -1 ∈ ∩ p≥1 L p . We namely have the inequality:

γ F λ δ = T 0 |D t F λ δ | 2 dt + T 0 |D ⊥ t F λ δ | 2 dt ≥ T 0 |D ⊥ t F λ δ | 2 dt = F λ δ ) 2 δ 2 T 2 > 0, (50)
and the lower bound, bounding from below the arithmetic mean by the geometric mean:

F λ δ = e δW ⊥ T 2 (1 -λ)S P T + d i=1 λα i S i,T ≥ e δW ⊥ T 2 (S P T ) (1-λ) d i=1 (S i,T ) λα i ,
which is a GBM. Then in view of ( 48)-(32), the next bound readily comes for j ∈ {1, 2} and any p ≥ 1:

||DF λ δ || j,p ≤ c |σ| ∞ √ T . (51)
Then [START_REF] Nualart | Malliavin calculus and related topics[END_REF]Proposition 1.5.6 and Proposition 2.1.4] proves the existence of Y λ 2 and Y λ 3 such that for any j ∈ {2, 3} and any p ≥ 1, the following estimates hold:

||Y λ j || p ≤ c ||Y|| i-1,p+ 1 2 ||DF λ δ || i-1 i-1,2p(2p+1) ||γ -1 F λ δ || i-1 i-1,2p (2p+1)
(52)

≤ c ||Y|| i-1,p+ 1 2 (|σ| ∞ √ T ) i-1 ||γ -1 F λ δ || i-1 i-1,2p (2p+1)
.

We now provide accurate estimates of γ F λ δ -1 . First observe that for the covariance matrix of F λ :

γ F λ = T 0 |D t S P T | 2 dt + λ 2 T 0 |D t R S ,1 T | 2 dt + 2λ T 0 D t S P T | D t R S ,1 T dt =(S P T ) 2 V P T + λ 2 T 0 |D t R S ,1 T | 2 dt + 2λS P T T 0 σ(t) | D t R S ,1
T dt, leading to the bound, using estimates (48) and (H σ ):

sup λ∈[0,1] γ F λ -(S P T ) 2 V P T p ≤ c |σ| ∞ µ[(M σ,. ) 2 ]T 3 2 , (53)
for any p ≥ 1. This intermediate result allows to prove the next lemma.

||(γ F η δ ) -1 || p ≤ c (V P T ) -1 , (54) sup t∈[0,T ],λ∈[0,1] D t (γ F λ δ ) -1 p ≤ c |σ| ∞ (V P T ) -1 , (55) sup s,t∈[0,T ],λ∈[0,1] ||D 2 s,t (γ F λ δ ) -1 || p ≤ c |σ| 2 ∞ (V P T ) -1 . (56)
Proof. First notice the useful inequality:

γ F λ δ =e 2δW T 2 T 0 |D t F λ | 2 dt + T 0 |D ⊥ t F λ δ | 2 > e 2δW T 2 T 0 |D t F λ | 2 dt = e 2δW T 2 γ F λ . (57)
For (54), use (57), Hölder and Markov inequalities to get for any p ≥ 1 and q ≥ 1 :

E (γ F λ δ ) -p = E (γ F λ δ ) -p 1 γ F λ ≤ 1 2 (S P T ) 2 V P T + E (γ F λ δ ) -p 1 γ F λ >(S P T ) 2 1 2 V P T ≤ (γ F λ δ ) -1 p 2p P (S P T ) 2 V P T -γ F λ 1 2 (S P T ) 2 V P T ≥ 1 + 1 2 V P T -p e -2δW T 2 (S P T ) -2 p p ≤ c δ 2 T ) -p γ F λ -(S P T ) 2 V P T q 2 2q 2(S P T ) -2 (V P T ) -1 q 2 2q + V P T -p ≤ c V P T -p (δ 2 T ) -p γ F λ -(S P T ) 2 V P T q 2 2q V P T p-q 2 + 1 . + 6 H 2 (d 1 ) V P T + H 3 (d 1 ) V P T 3 2 , G 6 = N (d 1 )
V P T 24 + 50 H 1 (d 1 )

V P T + 35 H 2 (d 1 ) V P T + 10 H 3 (d 1 ) V P T 3 2 + H 4 (d 1 ) V P T 2 ,
where for any n ∈ N, H n represents the n-th Hermite polynomial defined by, for any x ∈ R :

H n (x) = (-1) n e x 2 /2 ∂ n x n (e -x 2 /2 ).

Thus one has the following identities:

G 2 = N (d 1 ) V P T , G 2 + G 3 = -N (d 1 ) H 1 (d 1 ) V P T , G 2 + 3G 3 + G 4 =N (d 1 )
H 2 (d 1 )

V P T 3 2 , G h 2 + 15G h 3 + 25G h 4 + 10G h 5 + G h 6 =N (d 1 )
H 4 (d 1 )

V P T 5 2
.

Then we provide the first and second sensitivities w.r.t. the volatility parameter σ(0, T ) = 1 T V P T :

Vega BS 1, σ(0, T ) 2 T, K =∂ σ(0,T ) Call BS 1, σ(0, T ) 2 T, K = √ T N (d 1 ), Vomma BS 1, σ(0, T ) 2 T, K =∂ σ(0,T ) Vega BS (1, σ(0, T ) 2 T, K) = Vega BS 1, σ(0, T ) 2 T, K σ(0, T ) d 1 d 2 = Vega BS 1, σ(0, T ) 2 T, K σ(0, T ) [ ln 2 (K) σ2 (0, T ) 2 T - σ2 (0, T )T 4 ].

B.2. Benchmark Price Approximations.

In this Subsection, we make explicit the price approximation formulas provided in [START_REF]Interest rate model calibration using semidefinite programming[END_REF]- [START_REF] Carmona | Generalizing the black-scholes formula to multivariate contingent claims[END_REF]- [START_REF] Gobet | Weak approximation of averaged diffusion processes[END_REF]- [START_REF] Landon | Stratégies de couverture presque optimale : théorie et applications[END_REF]- [START_REF] Ju | Pricing asian and basket options via taylor expansion[END_REF]- [START_REF] Shiraya | Pricing multi-asset cross-currency options[END_REF] using our notations to allow easy comparisons. In addition, when relevant, we propose modifications of some formulas combining the methodologies of the authors and our martingale point of view.

£ D'Aspremont approximation formula. The second order approximation of E (S T -K) + provided in [14, Propositions 10 and 12 p.13 and p.15] is given by :

App 2 (DAF) =Call BS (1, V P T , K) (58)
+ d i=1 α i T 0 ξ i (t) | σ(t) V P T e 2 t 0 ξ i (s)| σ(s) ds N -ln(K) + 1 2 V P T + t 0 ξ i (s) | σ(s) ds V P T dt,
introducing the residual volatilities of the assets w.r.t. the Proxy volatility :

ξ i (t) := σ i (t) -σ(t). ( 59
)
The corrective term is a sum of temporal integrals of Gaussian functions and is a little bit less tractable than our second order formula whose corrective term is a linear combination of Gaussian functions computed at the maturity T . There is a priori no closed formula and one needs to compute numerically the integrals.

£ Modification of the D'Aspremont approximation formula. Actually, according to us, it seems that the Thesis [START_REF]Interest Rate Model Calibration and Risk-Management Using Semidefinite Programming[END_REF] and the paper [START_REF]Interest rate model calibration using semidefinite programming[END_REF] contain an error which we reported to the author, but not verified by the latter. We propose here a modification of the formula provided in [START_REF]Interest rate model calibration using semidefinite programming[END_REF] by following the approach developed by author. The idea is to work directly on the dynamic of the basket and to consider the stochastic weights α i :

dS t =S t d i=1 α i (t) σ i (t) | dW t , α i (t) = α i S i,t S t , α i (0) =α i ,
The dynamic of the stochastic weights is given by an application of the Itô formula :

d α i (t) = α i S i,t S t σ i (t) | dW t - α i S i,t S t d j=1 α j (t) σ j (t) | dW t + α i S i,t S t d j,k=1 α j (t) α k (t) σ j (t) | σ k (t) dt - α i S i,t S t d j=1 α j (t) σ i (t) | σ j (t) dt = α i (t) σ i (t) - d j=1 α j (t)σ j (t) | dW t - d j=1 α j (t)σ j (t)dt = α i (t) ξ i (t) - d j=1 α j (t)ξ j (t) | dW t -σ(t)dt - d j=1 α j (t)ξ j (t)dt ,
with the residual volatilities ξ i (t) defined in (59) and using the fact that d j=1 α j (t) = 1. In the paper [START_REF]Interest rate model calibration using semidefinite programming[END_REF]Equation (18) page 13] and in the Thesis [13, top of the page ix], the dynamic of the stochastic weights is however indicated as:

d α i (t) = α i (t) ξ i (t) - d j=1 α j (t)ξ j (t) | dW t + σ(t)dt+ d j=1 α j (t)ξ j (t)dt
Then the idea is to parametrize the basket and the stochastic weights dynamics using, similarly to our approach, a linear interpolation between the volatility of the assets and the Proxy volatility σ, replacing the vector of small magnitude d j=1 α j (t)ξ j (t) (because d j=1 α j ξ j (t) is null) by ε d j=1 α j (t)ξ j (t):

dS ε t = S ε t σ(t) -ε d i=1 α ε i (t)ξ i (t) | dW t , d α ε i (t) = α ε i (t) ξ i (t) -ε d j=1 α ε j (t)ξ j (t) | dW t -σ(t)dt -ε d j=1 α ε j (t)ξ j (t)dt , (60)
The next step is to perform an expansion of the pricing equation around ε = 0 using a PDE approach :

P ε (t, S ε t , α ε (t)) = E (S ε T -K) + |(t, S ε t , α ε (t)) = P 0 + P (1) ε + O(ε 2 ).
P 0 is given by the Black-Scholes price whereas, starting from the dynamic (60) and following the proof of [START_REF]Interest rate model calibration using semidefinite programming[END_REF]Lemma 11], we find that P (1) solves :

             -d i=1 ξ i (t) | σ(t) y i x 2 ∂ 2 x 2 P 0 = ∂ t P (1) + | σ(t)| 2 x 2 2 ∂ 2 x 2 P (1) + d i=1 ξ i (t) | σ(t) xy i ∂ 2 xy i P (1) + d i=1 |ξ i (t)| 2 y 2 i 2 ∂ 2 y 2 i P (1) -d i=1 ξ i (t) | σ(t) y i ∂ y i P (1) , P (1) = 0, for t = T.
The difference with [START_REF]Interest rate model calibration using semidefinite programming[END_REF] is the minus sign on the last term of the r.h.s. . Then following the proof of [14, Lemma 12], the Feynman-Kac representations give :

P (1) (0, S 0 , α) = T 0 E d i=1 ξ i (t) | σ(t) α P i (t)S P t N ln( S P t K )+ 1 2 V P t,T V P t,T V P t,T dt with V P t,T = T t | σ(s)| 2 ds
and the Proxy diffusions α P j and S P , solutions of (60) with ε = 0. In particular, observe that the product α P i S P is a log-normal martingale with volatility σ + ξ i . The Cameron-Martin formula and elementary calculus on Gaussian functions give that the corrective terms writes:

P (1) (0, S 0 , α) = d i=1 α i T 0 ξ i (t) | σ(t) V P T N -ln(K) + 1 2 V P T + t 0 ξ i (s) | σ(s) ds V P T dt
It is simpler that the corrective term in (58) and in addition one can compute explicitly the integral to get the final approximation formula:

App 2 (MDAF) =Call BS (1, V P T , K) + d i=1 α i N d 1 + C i V P T -N(d 1 ) , (61)
with the coefficients C i defined in Theorem 5. Notice that a Taylor expansion yields the next approximation, using the fact that d i=1 α i C i = 0:

d i=1 α i N d 1 + C i V P T -N(d 1 ) = 1 2 d i=1 α i N (2) (d 1 ) C 2 i V P T + d i=1 α i C 3 i (V P T ) 3 2 1 0 N (3) d 1 + λ C i V P T (1 -λ) 2 2 dλ = - 1 2 N (d 1 ) H 1 (d 1 ) V P T d i=1 α i C 2 i + O µ[(M σ,. ) 3 ]T 3 2 .
Hence up to an error of order O µ[(M σ,. ) 3 ]T 3 2 , the approximation (61) coincides with our second order approximation formula (35). This illustrates the duality between the stochastic analysis approach and the PDE techniques.

£ Carmona and Durrleman approximation formulas. The paper [12, Propositions 3 and 6 p.5 and p.7] provides lower and upper bounds of E (S T -K) + . Though the formulas are not strictly analytical but solutions of maximization/minimization problems, the functions to be optimized are interestingly very close to the one-dimensional Black-Scholes formula. One denotes by Σ the covariance matrix of the assets and by Σ k i = Σ ii + Σ kk -2Σ ik for any i ∈ {1, . . . , d}, k ∈ {0, . . . , d} where the index k = 0 is related to the variance/covariances involving the strike which are set to 0. The bounds are solutions of:

App(CDF) inf =sup y∈R sup u Σu=1 d i=1 α i N y + (Σu) i √ T -KN(y) , (62) App(CDF) sup = min k∈{0,...,d} d i=1 α i N d k + Σ k i T -KN d k -Σ kk T , (63)
with d k solution of:

d i=1 α i N d k + Σ k i T -KN d k -Σ kk T = 0.
Consider the i.i.d. case with equal weights. A symmetry argument leads to u = 1/(σ √ d)(1, . . . , 1) * (with σ the common volatility coefficient of the assets) and the solution of the optimization problem is y

= -ln(K)/( σ √ T ) -1 2 σ √ T with the proxy volatility σ = σ/ √ d.
Hence App(CDF) inf is simply given by the Black-Scholes price Call BS (1, V P T , K) and coincides with our second order formula showing in particular that this is a lower bound in the i.i.d. case. On the other hand for the upper bound, one has the identity d i=1 α

i N d 0 + Σ k i T -KN(d 0 ) = Call BS (1, σ 2
T, K). £ Gobet and Miri approximation formula. We recall the notation

S GM T = e T 0 σ(s)|dW s -1 2 d i=1 T 0 α i |σ i (s)| 2 ds
for the geometric mean proxy used in [START_REF] Gobet | Weak approximation of averaged diffusion processes[END_REF]. The link with our proxy S P is given by the following equation:

S GM T =S P T e λ , with λ = - 1 2 T 0 d j=1 α j |σ j (t)| 2 - d j=1 α j σ j (t) 2 dt = - 1 4
i, j∈{1,...,d} α i α j C i, j ≤ 0, (64)

C i, j = T 0 |σ i (t) -σ j (t)| 2 dt, i, j ∈ {1, . . . , d}
Using the change of variables described in Remark 1 to take into account the drift term e λ , the second order approximation of E (S T -K) + provided in [17, Theorem 2.2 p.12] reads as, in the particular case of log-normal processes: App 2 (GMF) =Call BS (e λ , V P T , K)

+e λ N( d 1 ) d i=1 α i 1 8 T 0 d j=1 α j |σ j (t)| 2 -|σ i (t)| 2 dt 2 + 1 2 T 0 | σ(t) -σ i (t)| 2 dt +e λ N ( d 1 ) V P T 1 2 d i=1 α i C i T 0 d j=1 α j |σ j (t)| 2 -|σ i (t)| 2 dt - 1 2 1 + H 1 ( d 1 ) V P T d i=1 α i C 2 i , (65)
with the coefficients C i defined in Theorem 5 and

d 1 = (λ -ln(K))/ V P T + 1 2 V P T .
The formula is a little bit more complicated than our corresponding second order approximation (35).

£ Landon approximation formula. In the particular case of basket options, the second order approximation of E (S T -K) + provided in [20, Theorem 8.4.2 p.160] can be written:

App 2 (LF) = d i=1 α i P i -KP 0 , P i = N( d i ) + λ N ( d i ) V P T , i ∈ {0, . . . , d} (66)
d i = -ln(K) -1 2 V P T + T 0 σ(t) | σ i (t) dt V P T , i ∈ {1, . . . , d}, d 0 = d 2 = -ln(K) -1 2 V P T V P T ,
where λ is defined in (64). This expansion reads as a multi-dimensional generalisation of the Black-Scholes formula with approximated probabilities of exercise. In the i.i.d. case with equal weights, the approximation is simply given by Call BS (1, V P T , K) = N(d 1 ) -KN(d 2 ) and coincides with our second order formula. We namely have d i = d 1 , ∀i ∈ {0, . . . , d} and the well known identity N (d 1 ) -KN (d 2 ) = 0. Finally notice that the term λ (appearing also in the approximation provided in [START_REF] Gobet | Weak approximation of averaged diffusion processes[END_REF]) is certainly due to the fact that the parametrization

S η t = d i=1 α i e t 0 σ η i (s)|dW s -1 2 t 0 η|σ i (s)| 2 +(1-η)| σ(s)| 2 ds
used in [START_REF] Landon | Stratégies de couverture presque optimale : théorie et applications[END_REF] does not maintain the martingale property, what could lead to numerical arbitrage. £ Modification of the Landon approximation formula. We propose here a modification of the approximation provided in [START_REF] Landon | Stratégies de couverture presque optimale : théorie et applications[END_REF] using our martingale parametrization. Following the approach developed in [START_REF] Landon | Stratégies de couverture presque optimale : théorie et applications[END_REF], we linearise the payoff, use changes of probability measure and perform approximations of probabilities of exercise regions. In the following heuristic calculus, P i for i ∈ {1, . . . , d} is the probability defined by the Radon-Nikodym derivative dP i dP | F T = S i,T under which S i k,. := S k,.

S i,. is a log-normal martingale with volatility σ k (.) -σ i (.), for any k ∈ {1, . . . , d}. One has to compute, setting S i . := S . S i,. :

E (S T -K) + = d i=1 α i P i ln( S i T ) > ln K S i,T -KP ln(S T ) > ln(K)
For h ∈ C 2 P (R), denoting by E P i the expectation under P i , S i,P T the log-normal Proxy with volatility σ -σ i and S i, (1) . the first order correction process ∀i ∈ {1, . . . , d}, we use the approximation :

E P i h ln( S i T ≈ E P i h ln( S i,P T + E P i h ln( S i,P T S i, (1)
T S i,P T = E P i h ln( S i,P T , using the fact that S i,(1) T = 0. Similarly we use the approximation E ln h S T ≈ E ln h S P T and we transpose the approach to the non-smooth indicator functions to get the next formula :

App 2 (MLF) = d i=1 α i P i ln( S i,P T > ln K S i,T -KP ln(S P T ) > ln(K) = d i=1 α i N( d i ) -KN(d 2 ).
As a result, one obtains a simpler formula. Finally notice (see (61)) that this formula exactly coincides with our modification of the D'Aspremont formula, remarking that

d i = d 1 + C i √ V P T
. We use the following unified notation for the resulting approximation : £ Ju approximation formula. In [START_REF] Ju | Pricing asian and basket options via taylor expansion[END_REF], the weighted average of log-normal variables with time homogeneous volatilities is approximated by a scalar log-normal random variable with a matching of the two first moments at the maturity T . The approximation is improved with a Taylor expansion of the ratio of the characteristic function of the basket to that of the proxy around zero volatility and takes the form of a third order formula :

α i α j σ i | σ j 2 T 2 , U (3) 2 =
i, j∈{1,...,d}

α i α j σ i | σ j 3 T 3 , U 3 =2
i, j,k∈{1,...,d}

α i α j α k σ i | σ k σ j | σ k T 2 , U 4 = 6
i, j,k∈{1,...,d}

α i α j α k σ i | σ k σ j | σ k 2 T 3 , U 5 =8
i, j,k∈{1,...,d}

α i α j α k σ i | σ j σ i | σ k σ j | σ k T 3 , U 6 =8
i, j,k,l∈{1,...,d}

α i α j α k α l σ i | σ l σ j | σ k σ k | σ l T 3 + 2U 2 U 2 , U 7 =6
i, j,k,l∈{1,...,d}

α i α j α k α l σ i | σ l σ j | σ l σ k | σ l T 3 .
The main term of the formula is the log-normal approximation given in [START_REF] Levy | Pricing european average rate currency options[END_REF]. In the one-dimensional case, the reader can verify that the identity a 2 = a 3 = d 1 = d 2 = d 3 = d 4 = z 1 = z 2 = z 3 = 0 holds and hence the approximation is exact. Notice that the formula is a little bit more complicated than our third order formula.

£ Shiraya and Takahashi approximation formula. We introduce the notations Y = 1 -K, d = Y/ V P T , C i = This reads as an expansion around the Bachelier (Gaussian) model. Observe that the corrective terms are the same than our third order formula (36) replacing C i , C i, j and d 1 by C i , C i, j and d respectively. This enlightens the link between the Black-Scholes model and the Bachelier model. However, one can notice that the approximation is inexact in the one-dimensional case. B.3. Benchmark implied volatility shapes ATM for short maturity. In this Subsection, we make explicit the implied volatility shape ATM for short maturity provided in [START_REF] Piterbarg | Markovian Projection Method for Volatility Calibration[END_REF]- [START_REF] Durrleman | From implied to spot volatilities[END_REF]- [START_REF] Alos | On the short-maturity behaviour of the implied volatility skew for random strike options and applications to option pricing approximation[END_REF] using our notations to allow easy comparisons.

£ Piterbarg formula. In [START_REF] Piterbarg | Markovian Projection Method for Volatility Calibration[END_REF], the dynamic of the basket (with time-homogeneous volatilities for the assets) is approximated with the following SDE: S such that ψ (1) = φ (1) -φ(1), we retrieve our formulas presented in Corollary 15 using the well-known rules that ATM for short maturity, the local volatility equals the implied volatility and its slope is twice the slope of the implied volatility (see for instance [START_REF] Berestycki | Asymptotics and calibration of local volatility models[END_REF]). Next the skew for ATM/short maturity is given by v 0 /(4|σ(0)| 2 ) = v 0 /(4| σ(0)| 2 ) where the stochastic process v t steps in the following decomposition of the dynamic of the variance process |σ(t)| 2 : d|σ(t)| 2 = δ t dt + v t dW S t + v ⊥ t dW S ,⊥ , for an extra BM W S ,⊥ independent of W S . An application of the Itô formula and a straightforward calculus of correlation leads to (we don't write the full decomposition for the sake of brevity) : This is our formula.

d|σ(t)| 2 = i,

Figure 1 .

 1 Figure 1. Estimates of the implied volatility by MC and second order formulas in the 10-dimensional correlated case.

Figure 2 .

 2 Figure 2. Estimates of the implied volatility by MC and third order formulas in the 10-dimensional correlated case.

Proof.

  Observe that the Malliavin derivative S P T is D t S P T = 1 t≤T S P T σ * (t) and ϕ S P T ∈ D 1,∞ with D t ϕ(S P T ) = 1 t≤T ϕ (1) S P T S P T σ * (t) = 1 t≤T Lϕ S P T σ * (t). Then identify T 0 c(t) | dW t with the Skorohod operator and apply the duality relationship [22, Definition 1.3.1 and Proposition 1.3.11].

Lemma 19 .

 19 Assume (H σ ). Then (γ F λ δ ) -1 ∈ D2,∞ and we have the next estimates for any p ≥ 1: sup λ∈[0,1]

App 2 (

 2

App 3 ((y 2 ) + z 2 y 2 σ 2 T N (y 2 ) + z 3 H 2 (y 2 ) σ 3 T 3 2 N1 =d 2 -d 3 + d 4 , z 2 = d 3 -d 4 , z 3 = d 4 ,d 1 =0. 5 (6a 2 1 + a 2 -

 322232232234234341512 JU) =N(y 1 ) -KN(y 2 ) (y 2 )] with the proxy volatility σ allowing the variance matching:σ = 1 T ln i, j∈{1,...,d} α i α j e T σ i |σ j , 4b 1 + 2b 2 ) -1/6[120a 3 1a 3 + 6(24c 1 -6c 2 + 2c 3c 4 )], d 2 =0.5(10a 2 1 + a 2 -6b 1 + 2b 2 ) -[128a 3 1 /3a 3 /6 + 2a 1 b 1a 1 b 2 + 50c 1 -11c 2 + 3c 3c 4 ], d 3 =2a 2 1b 1 -1/3[88a 3 1 + 3a 1 (5b 1 -2b 2 ) + 3(35c 1 -6c 2 + c 3 )], d 4 = -20a 3 1 /3 + a 1 (-4b 1 + b 2 ) -10c 1 + c 2 , c 1 =a 1 b 1 , c 2 = 1/144(9U 6 + 4U 7 ), c 3 = 1/48(4U 4 + U 5 ), c 4 = a 1 a 2 -2/3a 3 1 -1/6a 3 , b 1 =0.25U 3 , b 2 = a 2 1 -0.5a 2 , a 1 = -0.5U 2 , a 2 = 2a 2 1 -0.5U 2 , a 3 = 6a 1 a 2 -4a 3 1 -0.5U (3) 2 , U 2 =i, j∈{1,...,d} α i α j σ i | σ j T, U 2 = i, j∈{1,...,d}

T0

  σ i (t) | σ(t) dt and C i, j = T 0 σ i (t) | σ j (t) dt for any i, j ∈ {1, . . . , d}. Then the third order approximation of E (S T -K) + provided in [25, Theorem 3.1 p.7] is given by: App 3 (STM) = YN( d) + V P T N

α

  dS t = φ(S t )dW S twhere W S t is a scalar BM and φ a R-valued function with the following properties:φ(S 0 ) =φ(1) = p := d i=1 α i α j σ i | σ j = | σ| φ (S 0 ) =φ (1) = q :i σ i | σ 2 = | σ| + 1 | σ| 3 d i=1 α i C i (0) 2 , with C i (0) = σ iσ | σ .Hence introducing the local volatility function ψ(S ) = φ(S )

£

  Durrleman formula.[START_REF] Durrleman | From implied to spot volatilities[END_REF] Theorem 3.1] says that the level of the implied volatility ATM for short maturity equals the current value of the spot volatility, that is| d i=1 α i S i,0 S 0 σ i (0)| = | σ(0)| what coincides with our results. Regarding the skew ATM for short maturity, we follow [16, Section 5] and adopt a stochastic volatility point of view by writing, for a scalar BM W S associated to the Basket:

αii

  i (t) σ i (t) -σ(t) | σ(t)2 and we finally obtain v 0 /(4| σ(0)| 2 ) = 1/(2| σ(0)| 3 ) d i=1 α i σ i (0)σ(0) | σ(0)2 what is exactly our formula.£ Alos and León formula.[START_REF] Alos | On the short-maturity behaviour of the implied volatility skew for random strike options and applications to option pricing approximation[END_REF] provides approximations of E[(e X T -K T ) + ] for a stock price e X T and a positive random strike K T . To recover our setting (with a put option instead of a call option), we set X T := k for the deterministic log-strike and K T := S T for the basket with volatilityσ(t) = d i=1 S i,t S t α i σ i (t). [2, Remark 6] gives that σ I (T, k)| (T =0,k=0) ≈ |σ(0)| = | σ(0)| whereas [2,Theorem 14] provides for the skew using the notations x = ln(S 0 ) for the log-spot and D for the Malliavin derivative operator :∂ k σ I (T, k)| (T =0,k=0) = -∂ x σ I (T, k)| (T =0,x=0) = 1 2|σ(0)| 2 σ(t) | (D t |σ(t)|) * | (t=0) = 1 4|σ(0)| 3 σ(t) | (D t |σ(t)| 2 ) * | (t=0) j∈{1,...,d} α i α j S i,t S j,t σ i (t) | σ j (t) σ i (t)α i σ i (0) | σ(0) σ(0) | σ i (0) -σ(0) α i σ(0) | σ i (0) -σ(0) 2 = 1 2| σ(0)| 3 d i=1 α i σ i (0)σ(0) | σ(0) 2 .

methodology. 2.2.1. Corrective processes. Approximations

  S i,T ≈ S P T are expected to be accurate if |σ| ∞ , M σ and T are globally small enough. Nevertheless, we can not reasonably expect E[h(S T )] ≈

	E[h(S P T )] to be solely accurate enough and we provide correction terms. To derive them, we lever-
	age the next interpolations:

Table 1

 1

	K (%)	80	86	89	92	95	98	100	102	105	108	112	116	124
	(MC)	9.66	9.66	9.67	9.67	9.67	9.68	9.68	9.68	9.69	9.69	9.69	9.70	9.71
	Vol(App2)	9.49	9.49	9.49	9.49	9.49	9.49	9.49	9.49	9.49	9.49	9.49	9.49	9.49
	(CDF) sup	30.00 30.00 30.00 30.00 30.00 30.00 30.00 30.00 30.00 30.00 30.00 30.00 30.00
	(GMF)	ND	4.99	7.66	8.13	8.37	8.52	8.59	8.65	8.72	8.78	8.84	8.88	8.91
	Vol(App3)	9.67	9.67	9.68	9.68	9.68	9.68	9.68	9.68	9.68	9.68	9.68	9.67	9.67
	App3Vol	9.68	9.68	9.68	9.68	9.68	9.68	9.68	9.68	9.68	9.68	9.68	9.68	9.68
	(JUF)	9.65	9.66	9.67	9.67	9.67	9.68	9.68	9.68	9.69	9.69	9.69	9.70	9.71
	(STF)	9.82	9.72	9.70	9.69	9.68	9.68	9.68	9.67	9.67	9.66	9.65	9.63	9.57

  Table 2 and we also provide two Figures 1 and 2 thereafter representing the comparison between the estimated implied volatility (MC) and the approximated implied volatilities of order two (excepting (CDF) sup because estimations are too large) and then of order three : 17.29 18.01 18.40 18.81 19.12 19.35 19.57 19.77 20.00 20.24 20.58 20.85 21.31 Vol(App2) 15.98 17.74 18.22 18.63 18.93 19.16 19.37 19.58 19.81 20.06 20.41 20.66 20.99 App2Vol 17.39 17.95 18.28 18.65 18.93 19.16 19.37 19.58 19.81 20.07 20.44 20.76 21.33 (DAF) 16.68 18.13 18.59 19.00 19.29 19.52 19.73 19.93 20.14 20.37 20.67 20.87 21.10 (CDF) inf 17.15 17.86 18.25 18.66 18.97 19.21 19.42 19.62 19.85 20.09 20.42 20.69 21.14 (CDF) sup 31.61 31.73 31.79 31.86 31.92 31.96 32.00 32.03 32.07 32.12 32.18 32.23 32.32

	K	66	75	81	88	94	99	104	109	115	122	133	143	163
	(MC)
	(LF)	ND	17.19 18.05 18.67 19.06 19.33 19.55 19.75 19.95 20.13 20.33 20.44 20.49
	(MDALF)	15.22 17.65 18.19 18.64 18.95 19.19 19.40 19.61 19.83 20.07 20.39 20.62 20.91
	(GMF)	ND	18.39 18.75 18.92 18.99 19.03 19.06 19.08 19.10 19.12 19.14 19.15 19.18
	Vol(App3)	17.44 18.04 18.42 18.81 19.12 19.36 19.58 19.78 20.01 20.24 20.57 20.83 21.26
	App3Vol	17.41 18.04 18.42 18.81 19.12 19.36 19.58 19.78 20.01 20.24 20.57 20.82 21.22
	(JUF)	17.02 18.13 18.51 18.85 19.12 19.32 19.52 19.72 19.96 20.22 20.61 20.92 21.37
	(STF)	17.41 17.85 18.38 18.82 19.13 19.36 19.58 19.78 20.00 20.24 20.61 20.88 20.54

Table 3

 3 Estimates of the implied volatility (%) in the 2-dimensional case with ρ = -90%. One has | σ| ≈ 6.71%, M σ = µ[M σ,. ] ≈ 28.50% and σ ≈ 9.03%.

	K (%)	78	83	87	91	94	96	100	102	105	109	114	119	127
	(MC)	11.29 11.37 11.44 11.52 11.60 11.67 11.76 11.82 11.92 12.08 12.30 12.55 13.01
	Vol(App2)	10.61 10.61 10.61 10.61 10.61 10.61 10.61 10.61 10.61 10.61 10.61 10.61 10.61
	(GMF)	ND	ND	9.00	9.46	9.63	9.75	9.83	9.87	9.93	9.98	10.04 10.08 10.13
	Vol(App3)	11.72 11.85 11.93 12.00 12.04 12.06 12.07 12.07 12.05 12.01 11.95 11.87 11.74
	App3Vol	12.07 12.07 12.07 12.07 12.07 12.07 12.07 12.07 12.07 12.07 12.07 12.07 12.07
	(JUF)	10.15 10.87 11.18 11.42 11.58 11.73 11.88 11.97 12.11 12.28 12.48 12.66 12.87
	(STF)	12.53 12.28 12.17 12.12 12.10 12.09 12.07 12.05 12.01 11.92 11.75 11.54 11.16

Table 4

 4 Estimates of the implied volatility (%) in the 2-dimensional case with ρ = -75%. One has | σ| ≈ 10.61%, M σ = µ[M σ,. ] ≈ 26.25% and σ ≈ 11.98%.

	K (%)	70	78	82	88	92	96	100	103	108	113	121	127	141
	(MC)	15.43 15.48 15.51 15.56 15.59 15.62 15.66 15.69 15.74 15.80 15.90 15.98 16.22
	Vol(App2)	15.00 15.00 15.00 15.00 15.00 15.00 15.00 15.00 15.00 15.00 15.00 15.00 15.00
	(GMF)	ND	13.80 14.18 14.40 14.48 14.53 14.56 14.58 14.61 14.63 14.66 14.68 14.70
	Vol(App3)	15.67 15.71 15.73 15.75 15.75 15.76 15.76 15.76 15.75 15.75 15.73 15.72 15.68
	App3Vol	15.76 15.76 15.76 15.76 15.76 15.76 15.76 15.76 15.76 15.76 15.76 15.76 15.76
	(JUF)	15.17 15.35 15.43 15.52 15.58 15.64 15.69 15.73 15.79 15.85 15.94 16.00 16.12
	(STF)	16.55 16.03 15.92 15.83 15.80 15.78 15.76 15.75 15.72 15.68 15.58 15.48 15.12

Table 5

 5 Estimates of the implied volatility (%) in the 2-dimensional case with ρ = -50%. One has | σ| ≈ 15.00%, M σ = µ[M σ,. ] ≈ 22.50% and σ ≈ 15.74%.

	K (%)	65	73	79	85	90	95	100	104	110	116	126	135	153
	(MC)	18.65 18.68 18.69 18.71 18.73 18.75 18.76 18.78 18.80 18.82 18.86 18.90 19.00
	Vol(App2)	18.37 18.37 18.37 18.37 18.37 18.37 18.37 18.37 18.37 18.37 18.37 18.37 18.37
	(GMF)	14.04 17.60 17.89 18.01 18.06 18.09 18.11 18.12 18.14 18.15 18.16 18.17 18.19
	Vol(App3)	18.78 18.79 18.79 18.80 18.80 18.80 18.80 18.80 18.80 18.80 18.79 18.79 18.78
	App3Vol	18.80 18.80 18.80 18.80 18.80 18.80 18.80 18.80 18.80 18.80 18.80 18.80 18.80
	(JUF)	18.57 18.63 18.67 18.70 18.73 18.75 18.77 18.79 18.82 18.84 18.88 18.91 18.96
	(STF)	19.40 19.01 18.91 18.86 18.83 18.81 18.80 18.80 18.78 18.75 18.70 18.64 18.21

Table 6

 6 Estimates of the implied volatility (%) in the 2-dimensional case with ρ = -25%. One has | σ| ≈ 18.37%, M σ = µ[M σ,. ] ≈ 18.75% and σ ≈ 18.80%.

	K (%)	62	71	77	84	89	95	100	105	111	118	129	139	159
	(MC)	20.33 20.35 20.36 20.38 20.39 20.40 20.41 20.42 20.43 20.44 20.47 20.49 20.54
	Vol(App2)

  17.48 19.56 19.77 19.86 19.90 19.92 19.94 19.95 19.95 19.96 19.97 19.98 19.99

	Vol(App3)	20.42 20.42 20.43 20.43 20.43 20.43 20.43 20.43 20.43 20.43 20.43 20.42 20.42
	App3Vol	20.43 20.43 20.43 20.43 20.43 20.43 20.43 20.43 20.43 20.43 20.43 20.43 20.43
	(JUF)	20.30 20.33 20.35 20.37 20.38 20.40 20.41 20.42 20.44 20.45 20.47 20.49 20.52
	(STF)	20.87 20.55 20.50 20.47 20.45 20.44 20.43 20.42 20.41 20.40 20.37 20.33 19.82

Table 7

 7 Estimates of the implied volatility (%) in the 2-dimensional case with ρ = -10%. One has | σ| ≈ 20.12%, M σ = µ[M σ,. ] ≈ 16.50% and σ ≈ 20.43%.

	K (%)	52	63	69	78	86	93	100	107	115	126	143	158	192
	(MC)	28.07 28.07 28.07 28.07 28.07 28.07 28.07 28.07 28.07 28.07 28.07 28.07 28.07
	Vol(App2)	28.06 28.06 28.06 28.06 28.06 28.06 28.06 28.06 28.06 28.06 28.06 28.06 28.06
	(GMF)	27.97 28.04 28.05 28.05 28.05 28.05 28.05 28.06 28.06 28.06 28.06 28.06 28.06
	Vol(App3)	28.07 28.07 28.07 28.07 28.07 28.07 28.07 28.07 28.07 28.07 28.07 28.07 28.07
	(JUF)	28.07 28.07 28.07 28.07 28.07 28.07 28.07 28.07 28.07 28.07 28.07 28.07 28.07
	(STF)	26.29 27.71 27.94 28.05 28.07 28.07 28.07 28.08 28.08 28.10 28.21 28.22 26.53

Table 12

 12 Estimates of the implied volatility (%) in the 2-dimensional case with ρ = 75%. One has | σ| ≈ 28.06%, M σ = µ[M σ,. ] ≈ 9.92% and σ ≈ 28.07%.

	K (%)	50	61	68	78	85	92	100	107	116	127	145	161	197
	(MC)	29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24
	Vol(App2)	29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24
	(GMF)	29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24
	Vol(App3)	29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24
	(JUF)	29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24
	(STF)	29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24 29.24

  d}, and η ∈ [0, 1]. Proof. Inclusions in D ∞,∞ are standard to justify under our assumptions. Assume w.l.o.g. that p ≥ 2. We only detail results for ||D r S η i,t || p and ||D r S || p . The other estimates, presenting no extra difficulty, are left to the reader. For S

		η,(k)
		i,t η i,t , one has D r S η i,t = 1 r≤t S	η i,t σ	η i (r) * and the estimate is obvious.
	For S	η,(k)

  j∈{1,...,d}α i (t) α j (t) σ i (t) | σ j (t) d α i (t) (t) σ i (t) | σ(t) σ i (t) -σ(t) | σ(t) = 2 |σ(t)|

	Hence the following identity comes using the definition of σ(t):
	v t =	2 |σ(t)|	d i=1	α i
						α i (t)	+	d α j (t) α j (t)	+ . . . dt ,
	d α i (t) α i (t)	= σ i (t) -σ(t) | dW t -σ(t)dt =	1 |σ(t)|	σ i (t) -σ(t) | σ(t) dW S t + . . . dW S ,⊥ t	+ . . . dt

The risk-free rate and the dividend yield are supposed null. Adaptation of the results for non-zero but deterministic interest rates and dividends is straightforward by considering the change of variable discussed in Remark 1 .

The presented approximations are not analytical, we use the heuristics that derivatives of the expansion should coincide with the expansions of the derivatives.

1 bp (basis point) equals 0.01%.

This is due to the fact that, at the second order, the main part of the approximation is captured by the main term, i.e. the Black-Scholes price as noticed in subsection 2.3.

We have generated a covariance matrix using the genPositiveDefMat function of the R package clusterGeneration with a random correlation matrix and uniform random variances in [0 : 16%] such that volatilities are in [0 : 40%] with mean close to 30%.

actually D 3,∞ is sufficient.

Acknowledgments. I am grateful to Emmanuel Gobet for fruitful discussions as well as his valuable feedback and for pointing out references [18], [16] and [2].

Table 11

Estimates of the implied volatility (%) in the 2-dimensional case with ρ = 50%. One has | σ| ≈ 25.98%, M σ = µ[M σ,. ] ≈ 12.99% and σ ≈ 26.02%.

Appendix A. Technical proofs.

A.1. Proof of Lemma 7. We proceed by induction and provide firstly a technical result:

Then choose q = 12p and use (32)-( 53) to obtain the final bound:

. This leads to the announced result, using (54):

Proof of (56) is very similar and is left to the reader.