Vincent Lemaire
email: vincent.lemaire@upmc.fr

Gilles Pagès
email: gilles.pages@upmc.fr

Multilevel Richardson-Romberg extrapolation

Keywords: Multilevel Monte Carlo estimator, Richardson-Romberg extrapolation, Multistep, Euler scheme, Nested Monte Carlo method, Stratification, Option pricing primary 65C05, secondary 65C30, 62P05

). Thus, in standard frameworks like discretization schemes of diffusion processes, the root mean squared error (RMSE) ε > 0 can be achieved with our ML2R estimator with a global complexity of ε -2 log(1/ε) instead of ε -2 (log(1/ε)) 2 with the standard MLMC method, at least when the weak error

). The ML2R estimator is then halfway between a regular MLMC and a virtual unbiased Monte Carlo. When the strong error

), β < 1, the gain of ML2R over MLMC becomes even more striking. We carry out numerical simulations to compare these estimators in two settings: vanilla and path-dependent option pricing by Monte Carlo simulation and the less classical Nested Monte Carlo simulation.

Introduction

The aim of this paper is to combine the multilevel Monte Carlo estimator introduced by S. Heinrich in [START_REF] Heinrich | Multilevel Monte Carlo methods[END_REF] and developed by M. Giles in [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF] (see also [START_REF] Kebaier | Statistical Romberg extrapolation: a new variance reduction method and applications to option pricing[END_REF] for the statistical Romberg approach) and the (consistent) Multistep Richardson-Romberg extrapolation (see [START_REF] Pagès | Multi-step Richardson-Romberg extrapolation: remarks on variance control and complexity[END_REF]) in order to minimize the simulation cost of a quantity of interest I 0 = E [Y 0 ] where the random variable Y 0 cannot be simulated at a reasonable cost (typically a generic multidimensional diffusion process or a conditional expectation). Both methods rely on the existence of a family of simulatable random variables Y h , h > 0, which strongly approximate Y 0 as h goes to 0 whose bias E [Y h ] -E [Y 0 ] can be expanded as a polynomial function of h (or h α , α > 0).

However the two methods suffer from opposite but significant drawbacks: the multilevel Monte Carlo estimator does not fully take advantage of the existence of such an expansion beyond the first order whereas the Multistep Richardson-Romberg extrapolation cannot prevent an increase of the variance of the resulting estimator. Let us be more precise.

Consider a probability space (Ω, A, P) and suppose that we have a family (Y h ) h 0 of real-valued random variables in L 2 (P) associated to Y 0 supposed to be non degenerate and satisfying lim h→0 Y h -Y 0 2 = 0 where h takes values in an admissible subset of parameters H ⊂ (0, h] such that 0 ∈ H, h ∈ H and H n ⊂ H for every integer n 1. We also assume that h ∈ H. Usually, the random variable Y h appears as a functional of a time discretization scheme of step h or from an inner approximation in a nested Monte Carlo simulation. The parameter h is called the bias parameter in what follows. Furthermore, we make the pseudo-assumption, that for every admissible h ∈ H, the random variable Y h is simulatable whereas Y 0 is not (at a reasonable cost). For this reason, the specification of h will be often made in connection with the complexity of the simulation Y h with in mind to make it inverse linear in h.

We aim at computing an as good as possible approximation of I 0 = E [Y 0 ] by a Monte Carlo type simulation. The starting point is of course to fix a parameter h ∈ H to consider a standard Monte Carlo estimator based on Y h to compute I 0 . So, let (Y (k) h ) k 1 be a sequence of independent copies of Y h and the estimator I

(h) N = 1 N N k=1 Y (k)
h . By the strong law of numbers and the central limit theorem we have a standard control of the renormalized statistical error √ N (I

N -E [Y h ]) which behaves as a centered Gaussian with variance var(Y h ). On the other hand, there is a bias error due to the approximation of I 0 by I h = E [Y h ]. This bias error is also known as the weak error when Y h is a functional of the time discretization scheme of a stochastic differential equation with step h. In many applications, the bias error can be expanded as

E [Y h ] -E [Y 0 ] = c 1 h α + • • • + c R h αR + o(h αR ) (1)
where α is a positive real parameter (usually α = 1 2 , 1 or 2). In this paper, we fully take into account this error expansion and provide a very efficient estimator which can be viewed as a coupling between an MLMC estimator and a Multistep Richardson-Romberg extrapolation.

We first present a brief description of the original MLMC estimator as described in [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF]. The main idea is to use the following telescopic summation with depth L 2

E [Y h L ] = E [Y h ] + L j=2 E Y h j -Y h j-1
where (h j ) j=0,...,L is a geometrically decreasing sequence of different bias parameters h j = M -(j-1) h. For each level j ∈ 1, . . . , L , the computation of E Y h j -Y h j-1 is performed by a standard Monte Carlo procedure. The key point is that, at each level, we consider a number N j = ⌈N q j ⌉ of scenarios where q = (q 1 , . . . , q L ) ∈ S + (L) = q ∈ (0, 1) L , L j=1 q j = 1 (L-dimensional simplex) and a random sample of Y h j and Y h j-1 supposed to be perfectly correlated. More precisely, we consider L copies of the biased family denoted Y (j) = (Y (j) h ) h∈H , j ∈ 1, . . . , L attached to independant random copies Y (j) 0 of Y 0 . The MLMC estimator then writes

I N h,L,q = 1 N 1 N 1 k=1 Y (1),k h + L j=2 1 N j N j k=1 Y (j),k h j -Y (j),k h j-1 (2)
where (Y (j),k ) k 1 , j = 1, . . . , L are independent sequences of independent copies of Y (j) and N 1 , . . . , N L are positive integers. The analysis of the computational complexity and the study of the bias-variance structure of this estimator will appear as a particular case of a generalized multilevel framework that we will introduce and analyze in Section 3. This framework, following the original MLMC, highly relies on the combination of a strong rate of approximation of Y 0 by Y h and a first order control of weak error

E [Y h ] -E [Y 0 ]
. This MLMC estimator has been extensively applied to various fields of numerical probability (jump diffusions [START_REF] Dereich | A multilevel Monte Carlo algorithm for Lévy-driven stochastic differential equations[END_REF][START_REF] Dereich | Multilevel Monte Carlo algorithms for Lévy-driven SDEs with Gaussian correction[END_REF], computational statistics and more general numerical analysis problems (high dimensional parabolic SPDEs, see [START_REF] Barth | Multilevel Monte Carlo method for parabolic stochastic partial differential equations[END_REF], etc). For more references, we refer to the webpage http://people.maths.ox.ac.uk/gilesm/mlmc_community.html and the references therein.

On the other hand, the multistep Richardson-Romberg extrapolation takes advantage of the full expansion (1). Let us first recall the one-step Richardson-Romberg Monte Carlo estimator. We consider one biased family denoted Y = (Y h ) h∈H attached to the random variable Y 0 . The one-step Richardson-Romberg Monte Carlo estimator then writes

I N h, h 2 = 1 N N k=1 2Y k h 2 -Y k h
where (Y k ) k 1 is a sequence of independent copies of Y . It is clear that this linear combination of Monte Carlo estimators satisfies the following bias error expansion (of order 2 in h)

E 2Y h 2 -Y h -E [Y 0 ] = - c 2 2 h 2 + o(h 2 ).
Moreover, the asymptotic variance of this estimator satisfies lim h→0 var(I N h, h 2 ) = var(Y 0 )/N which is the same as the crude Monte Carlo estimator. It is natural to reiterate this extrapolation to obtain a linear estimator with bias error of order 3 in h and so on. This extension called Multistep Richardson-Romberg extrapolation for Monte Carlo estimator has been introduced and extensively investigated in [START_REF] Pagès | Multi-step Richardson-Romberg extrapolation: remarks on variance control and complexity[END_REF] in the framework of discretization of diffusion processes. More details are given in Section 2.4. The aim of this paper is to show that an appropriate combination of the MLMC estimator and the Multistep Richardson-Romberg extrapolation outperforms the standard MLMC. More precisely, we will see in Section 3 that an implementation of the Multilevel Richardson Romberg estimator (ML2R) turns out to be a weighted version of MLMC and writes

I N h,R,q = 1 N 1 N 1 k=1 Y (1),k h + R j=2 W j N j N j k=1 Y (j),k h j -Y (j),k h j-1 (3)
where R 2 is the depth level -similar to L in (2) -and (Y (j),k ) k 1 are like in (2). We denote by n j = M j-1 the j-th refiner coefficient of the initial bias parameter h ∈ H and call M 2 the root of the refiners. A strong feature of our approach comes from the fact that the weights (W k ) k=2,...,R are explicit and only depend on α (given by (1)), M 2 and R 2. In practice these ML2R weights read W j = R i=j w i where w i are given by (12). These derivative weights (w i ) i∈{1,...,R} have been introduced in [START_REF] Pagès | Multi-step Richardson-Romberg extrapolation: remarks on variance control and complexity[END_REF] to kill the successive bias terms that appear in the expansion (1).

To compare the two methods MLMC and ML2R, we consider the following optimization problem: minimizing the global simulation cost (of one estimator) subject to the constraint that the resulting L 2 -error or root mean squared error (RMSE) must be lower than a prescribed ε > 0. We solve the problem step by step for both estimators (in fact for a more general unifying class of estimators). The first stage is a stratification procedure to minimize the effort estimator (product of the variance by the complexity) to optimally dispatch the N j across all level. Doing so we are able to specify the initial bias parameter h and the depth level parameter R as a function of ε and structural fixed parameters (α, β, V 1 , var(Y 0 )). A light preprocessing makes possible to optimize the choice of the root M 2 of the refiners n j , j = 1, . . . , R. Basically (see Theorem 3.11), the numerical cost of the ML2R estimator implemented with these optimal parameters (depending on ε) and denoted Cost(ML2R) satisfies lim sup

ε→0 v(β, ε) × Cost(ML2R) K(α, β, M )
where K(α, β, M ) is an explicit bound and

v(β, ε) =        ε 2 log(1/ε) -1 if β = 1, ε 2 if β > 1, ε 2 e -1-β √ α √ 2 log(1/ε) log(M ) if β < 1.
As first established in [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF], we prove in a similar way that the optimal numerical cost of the MLMC estimator denoted Cost(MLMC) satisfies a similar result with v(β, ε)

= ε 2 log(1/ε) -2 if β = 1 and v(β, ε) = ε 2+ 1-β α if β < 1.
In the case β = 1, the gain of log(1/ε) may look as a minor improvement but, beyond the fact that it is halfway to a virtual unbiased simulation, this improvement is obtained with respect to a tremendously efficient method to speed up crude Monte Carlo simulation. In fact, as emphasized in our numerical experiments (see Section 5), this may lead to a significant reduction factor for CPU time, e.g. when α < 1: pricing a Black-Scholes Lookback Call option with a prescribed quadratic error ε = 2 -8 , yields a reduction factor of 3.5 in favor of MLMC. When β < 1, ML2R the above theoretical reduction factor asymptotically goes to +∞ as ε goes to 0 in a very steep way. Thus, the reduction in CPU time factor reaches, mutatis mutandis, 22 for a Black-Scholes Up&Out Barrier call option for which β = 1 2 (still using a regular Euler scheme without Brownian bridge). When compared on the basis of the resulting empirical RMSE, these factors become even larger (approximately 48 and 61 respectively). In fact, it confirms that β < 1 is the setting where our Multilevel Richardson-Romberg estimator is the most powerful compared to regular MLMC.

The paper is organized as follows: in Section 2, we propose a general parametrized framework to formalize the optimization of a biased Monte Carlo simulation based on the L 2 -error minimization. The crude Monte Carlo estimator and the multistep Richardson-Romberg estimator appear as the first two examples, allowing us to make precise few notations as well as our main assumptions. In Section 3, we first introduce the extended family of multilevel estimators attached to allocations matrix T. Among them, we describe in more details our proposal: the new ML2R estimator, but also the standard MLMC estimator. Two typical fields of application are presented in Section 4: the time discretization of stochastic processes (Euler scheme) and the nested Monte Carlo method, for which a weak expansion of the error at any order is established in the regular case. In Section 5, we present and comment on numerical experiments carried out in the above two fields.

Notations: • Let N * = {1, 2, . . .} denote the set of positive integers. • If n = (n 1 , . . . , n R ) ∈ (N * ) R , |n| = n 1 + • • • + n R and n! = 1 i R n i .
• Let (e 1 , . . . , e R ) denote the canonical basis of R R (viewed as a vector space of column vectors). Thus e i = (δ ij ) 1 j r where δ ij stands for the classical Kronecker symbol.

• ., . denotes the canonical inner product on R R .

• For every x ∈ R + = [0, +∞), ⌈x⌉ denotes the unique n ∈ N * satisfying n -1 < x n and ⌊x⌋ denotes the unique n ∈ N satisfying n x < n + 1.

• For every x ∈ R + = [0, +∞), and y ∈ R * + = (0, +∞), ⌈x⌉ y denotes the unique n ∈ N * satisfying (n -1)y < x ny and ⌊x⌋ y denotes the unique n ∈ N satisfying ny x < (n + 1)y.

• If (a n ) n∈N and (b n ) n∈N are two sequences of real numbers, a n ∼ b n if a n = ε n b n with lim n ε n = 1, a n = O(b n ) if (ε n ) n∈N is bounded and a n = o(b n ) is lim n ε n = 0.

Preliminaries

Mixing variance and complexity (effort)

We first introduce some notations and recall basic facts on (possibly biased) linear estimators. We consider a family of linear statistical estimator (I N π ) N 1 of I 0 ∈ R where π lies in a parameter set Π ⊂ R d . By linear, we mean, on the one hand, that

E I N π = E I 1 π , N 1,
and, on the other hand, that the numerical cost Cost(I N π ) induced by the simulation of I N π is given by

Cost(I N π ) = N κ(π)
where κ(π) = Cost(I 1 π ) is the cost of a single simulation or unitary complexity. We also assume that our estimator is of Monte Carlo type in the sense that its variance is inverse linear in the size N of the simulation:

var(I N π ) = ν(π) N where ν(π) = var(I 1 π ) denotes the variance of one simulation. For example, in a crude biased Monte Carlo π = h ∈ H, in a Multilevel Monte Carlo π = (h, R, q) ∈ H ×N * × S + and in the Multistep Monte Carlo π = (h, R) ∈ H ×N * .
We are looking for the "best" estimator in this family (I N π ) N 1 , π ∈ Π i.e. the estimator minimizing the computational cost for a given error ε > 0. In the sequel, we will consider N as a continuous variable lying in R + . A natural choice for measuring the random error I N π -I 0 is to consider the L 2 -error or root mean squared error (RMSE) E (I N π -I 0 ) 2 = I N π -I 0 2 . Our aim is to minimize the cost of the simulation for a given target error, say ε > 0. This generic problem reads π(ε), N (ε) = argmin

I N π -I 0 2 ε Cost(I N π ). (4)
In order to solve with this minimization problem, we introduce the notion of effort φ π of a linear Monte Carlo type estimator I N π .

Definition 2.1. The effort of the estimator I N π is defined for every π ∈ Π by

φ π = ν(π) κ(π). (5)
By definition of a linear estimator I N π we have that

φ π = ν(π) κ(π) = var(I N π ) Cost(I N π ) = var(I 1 π ) Cost(I 1 π )
for every integer N 1, so that we obtain the fundamental relation

Cost(I N π ) = N φ π ν(π) . (6)
• If the estimators (

I N π ) N 1 are unbiased i.e. E I N π = I 0 for every π ∈ Π, then E (I N π -I 0 ) 2 = I N π -I 0 2 2 = var I N π = 1 N ν(π).
The solution of the generic problem (4) then reads

π(ε) = π * = argmin π∈Π φ π , N (ε) = ν(π * ) ε 2 = φ(π * ) κ(π * )ε 2 . (7)
Consequently, the most performing estimator I N π is characterized as a minimizer of the effort φ π as defined above (and the parameter π does not depend on ε).

• When the estimators (I N π ) N 1 , π ∈ Π, are biased, the mean squared error writes

E (I N π -I 0 ) 2 = µ 2 (π) + ν(π) N where µ(π) = E I N π -I 0 = E I 1 π -I 0
denotes the bias (which does not depend on N ). Using that ν(π) = N I N π -I 0 2 2µ(π) 2 , the solution of the generic problem (4) reads

π(ε) = argmin π∈Π, |µ(π)|<ε φ π ε 2 -µ 2 (π) , N (ε) = ν(π(ε)) ε 2 -µ 2 (π(ε)) = φ(π(ε)) κ(π(ε))(ε 2 -µ 2 (π(ε))) . ( 8
)

Assumptions on weak and strong approximation errors

We come back to the framework described in the introduction: let (Y h ) h∈H be a family of real-valued random variables associated to a random variable Y 0 ∈ L 2 . The index set H is a consistent set of step parameters in the sense that H ⊂ (0, h], h ∈ H and, for every integer n 1, H n ⊂ H (hence 0 is a limiting value of H). All random variables Y h are defined on the same probability space (Ω, A, P). The family satisfies two assumptions which formalize the strong and weak rates of approximation of Y 0 by Y h when h → 0 in H. These assumptions are the basement of multilevel simulation methods (see [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF][START_REF] Heinrich | Multilevel Monte Carlo methods[END_REF]):

Bias error expansion (weak error rate):

∃ α > 0, R 1, E [Y h ] = E [Y 0 ] + R k=1 c k h αk + h α Rη R (h), lim h→0 η R (h) = 0, (W E α, R)
where c k , k = 1, . . . , R are real coefficients and η R is a real valued function defined on H.

Strong approximation error assumption:

∃ β > 0, V 1 ∈ R + , Y h -Y 0 2 2 = E Y h -Y 0 2 V 1 h β . (SE β )
Note that the parameters α, β and R are structural parameters which depend on the family (Y h ) h∈H . When (Y h ) h∈H satisfies (W E α, R) for every integer R, we will say that (W E α,∞ ) is fulfilled. Such a family is said to be admissible (at level R with parameters β and α).

Note that when c 1 = 0, consistency of strong and weak errors implies in what follows that β 2α. In the sequel we will consider a free parameter R ∈ 2, . . . , R for which (W E α, R) is always satisfied (with the same coefficients c r up to r = R). This parameter R corresponds to the depth level L used in the multilevel literature and will also be referred to as the depth from now on.

In what follows, we will use the following ratio

θ = V 1 var(Y 0 ) (9)
which relates the quadratic rate of convergence of Y h to Y 0 and the variance of Y 0 All estimators considered in this work are based on independent copies (Y

(j)
h ) h∈H , (attached to random variables Y (j) 0 ) of (Y h ) h∈H , j = 1, . . . , R. All random variables are supposed to be defined on the same probability space. Note that, since the above properties (SE β ) and (W E α, R), R 1, only depend on the distribution of (Y h ) h∈H , all these copies will also satisfy these two properties.

We associate to the family (Y h ) h∈H and a given bias parameter h ∈ H, the R R -valued random vector

Y h,n = (Y h , Y h n 2 , . . . , Y h n R
)

where the R-tuple of integers n := (n 1 , n 2 , . . . , n R ) ∈ N R , called refiners in the sequel, satisfy

n 1 = 1 < n 2 < • • • < n R .
One defines likewise Y (j) h,n for the (independent) copies (Y

(j) h ) h∈H .
✄ Specification of the refiners: In most applications, we will choose refiners n i as n i = M i-1 where M 2. Indeed, this is the standard choice in the regular Multilevel Monte Carlo method as described in [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF]. Other choices like n i = i are possible (see below).

Crude Monte Carlo estimator

In our formalism a crude Monte Carlo simulation and its cost can be described as follows.

Proposition 2.2. Assume (W E α, R) with R 1 and c 1 = 0. The Monte Carlo estimator of

E [Y 0 ] defined by ∀N 1, h ∈ H, Ȳ N h = 1 N N k=1 Y k h where Y k h k 1 is an i.i.d. sequence of copies of Y h , satisfy µ(h) = c 1 h α (1 + η 1 (h)), κ(h) = 1 h , φ(h) = var(Y h ) h
and, for a prescribed L 2 -error ε > 0, the optimal parameters h * (ε) and N * (ε) solution to (4) are given by

h * (ε) = (1 + 2α) -1 2α ε |c 1 | 1 α , N * (ε) = 1 + 1 2α var(Y 0 )(1 + θ(h * (ε)) β 2 ) 2 ε 2
.

Furthermore, we have

lim sup ε→0 ε 2+ 1 α min h∈H, | µ(h)|<ε Cost( Ȳ N h ) |c 1 | 1 α 1 + 1 2α (1 + 2α) 1 2α var(Y 0 ).
Proof. The proof is postponed to Annex B.

We refer to the seminal paper [START_REF] Duffie | Efficient Monte Carlo simulation of security prices[END_REF] for more details on practical implementation of this estimator.

Remark 2.3. For crude Monte Carlo simulation, Assumption (SE β ) is not necessary. Note that, at order 1, one can always assume c 1 = 0 considering the first non-zero term h α in the expansion (if any).

Background on Multistep Richardson-Romberg extrapolation

The so-called Multistep Richardson-Romberg estimator has been introduced in [START_REF] Pagès | Multi-step Richardson-Romberg extrapolation: remarks on variance control and complexity[END_REF] in the framework of Brownian diffusions. It relies on R (refined) Euler schemes X( h n i

) , 1 i R, defined on a finite interval [0, T ] (T > 0) where the bias parameter h = T n , n 1. In that case, the refiners are set as n i = i, i = 1, . . . , R, (in order to produce a better control of both the variance and the complexity for the proposed estimator, see Remark 2.5 below). The main results are obtained when all the schemes are consistent i.e. all the Brownian increments are generated from the same underlying Brownian motion. As a consequence, under standard smoothness assumptions on the coefficients of the diffusion, the family Y h = X(h) , h ∈ H, makes up an admissible family in the above sense as will be seen further on in more details.

For a refiner vector (n 1 , n 2 , . . . , n R ) we define the weight vector w = (w 1 , . . . , w R ) as the unique solution to the Vandermonde system V w = e 1 where

V = V (1, n -α 2 , . . . , n -α R ) =      1 1 • • • 1 1 n -α 2 • • • n -α R . . . . . . • • • . . . 1 n -α(R-1) 2 • • • n -α(R-1) R      .
The solution w of the system has a closed form given by Cramer's rule (see Lemma A.1 in Appendix A):

∀i ∈ 1, . . . , R ,

w i = (-1) R-i n α(R-1) i 1 j<i (n α i -n α j ) i<j R (n α j -n α i ) . ( 10
)
We also derive the following identity of interest

w R+1 := R i=1 w i n αR i = (-1) R-1 n! α . (11)
Note that all coefficients (w i ) 1 i R depend on the depth level R of the combined extrapolations. For the standard choices n i = i or n i = M i-1 , i = 1, . . . , R, we obtain the following expressions:

w i =                (-1) R-i i αR i-1 j=1 (i α -j α ) R i+1 (j α -i α ) if n j = j, j ∈ 1, . . . , R , (-1) R-i M -α 2 (R-i)(R-i+1) i-1 j=1 (1 -M -jα ) R-i j=1 (1 -M -jα ) if n j = M j-1 , j ∈ 1, . . . , R . (12)
Note that when α = 1 and n j = j, then

w i = (-1) R-i i R i!(R -i)! , i = 1, . . . , R.
Assume now (W E α, R) and R ∈ 1, . . . , R . In order to design an estimator which kills the bias up to order R, we focus on the random variable resulting from the linear combination w,

Y h,n = R i=1 w i Y h n i . The first equation of the Vandermonde system V w = e 1 , namely R r=1 w r = 1, implies that lim h→0 E w, Y h,n = E [Y 0 ] .
Furthermore, when expanding the (weak) error, one checks that the other R -1 equations satisfied by the weight vector w make all terms in front of the c r , r = 1, . . . , R -1 vanish. Finally, we obtain

E w, Y h,n = E [Y 0 ] + c R w R+1 h αR 1 + η R,n (h) (13) where η R,n (h) = 1 c R w R+1 R r=1 w r n αR r η R h n r → 0 as h → 0. (14) Proposition 2.4. Assume (W E α, R) and R ∈ 2, . . . , R . The Multistep Richardson-Romberg esti- mator of E [Y 0 ] defined by ∀N 1, h ∈ H, Ȳ N h,n = 1 N N k=1 w, Y k h,n = w, 1 N N k=1 Y k h,n (15)
where

Y k h,n k 1 is an i.i.d. sequence of copies of Y h,n , satisfies µ(h) = (-1) R-1 c R h R n! α 1 + η R,n (h) , κ(h) = |n| h , φ(h) = |n| var( w, Y h,n ) h
and, for a prescribed L 2 -error ε > 0, the optimal parameters h * (ε) and N * (ε) solution of (4) are

h * (ε) = (1 + 2αR) -1 2αR ε |c R | 1 αR n! 1 R , N * (ε) = 1 + 1 2αR var(Y 0 )(1 + θ(h * (ε)) β 2 ) 2 ε 2 . Furthermore, inf h∈H | µ(h)|<ε Cost( Ȳ N h ) ∼ (1 + 2αR) 1+ 1 2αR 2αR |c R | 1 αR n var(Y 0 ) n! 1 R ε 2+ 1 αR as ε → 0. ( 16
)
Proof. The proof is postponed to Appendix B (but takes advantage of the formalism developed in the next Section).

Remark 2.5. In this approach the bias reduction suffers from an increase of the simulation cost by a |n| factor which appears in the numerator of (16). The choice of the refiners in [START_REF] Pagès | Multi-step Richardson-Romberg extrapolation: remarks on variance control and complexity[END_REF], namely

n i = i, i = 1, . . . , R, is justified by the control of the ratio n n! 1 R
: for such a choice it behaves linearly in R -like e 2 (R + 1) -for large values of R, whereas with

n i = M i-1 it goes to infinity like M R-1 2 .
3 A paradigm for Multilevel simulation methods

General framework

Multilevel decomposition

In spite of Proposition 2.4 which shows that the numerical cost of the Multistep method behaves like ε 2+ 1 αR , one observes in practice that the increase ot the ratio |n| n! (when R grows) in front of var(Y 0 ) in ( 16) reduces the impact of the bias reduction.

An idea is then to introduce independent linear combination of copies of Ȳh,n to reduce the variance taking advantage of the basic fact that if X and X ′ are independent with the same distribution then

E X+X ′ 2 = E [X] and var( X+X ′
2 ) = 1 2 var(X), combined with an appropriate stratification strategy to control the complexity of the resulting estimator. So, let us consider now R independant copies (Y (j) h,n ), j = 1, . . . , R, of the random vector Y h,n and the linear combination

R j=1 T j , Y (j) h,n = R i,j=1 T j i Y (j) h n i
where

T = [T 1 . . . T R ] is an R × R matrix with column vectors T j ∈ R R satisfying the constraint 1 i,j R T j i = 1. Under Assumption (W E α, R) and R ∈ 2, . . . , R , E R j=1 T j , Y (j) h,n -E [Y 0 ] = R j=1 T j , E Y (j) h,n -E [Y 0 ] = R i=1 R j=1 T j i E Y h n i -E [Y 0 ] = R i=1 R j=1 T j i E Y h n i -E [Y 0 ] = o(h α ).
The strong L 2 -convergence of the estimator also holds (without rate) as soon as Y h strongly converges toward Y 0 (in L 2 ).

As emphasized further on, we will also need that each column vector T j , j ∈ 2, . . . , R, has zero sum. In turn, this suggests to introduce the notion of Multilevel estimator as a family of stratified estimators of E [Y 0 ] attached to the random vectors T j , Y (j) h,n , j = 1, . . . , R. This leads to the following definitions.

Definition 3.1 (Allocation matrix). Let R 2. An R × R-matrix T = [T 1 . . . T R ] is an R-level allocation matrix if T j , 1 = R i=1 T j i = 0, j = 2, . . . , R. (17)
Note that such an allocation matrix always satisfies d i,j=1

T j i = 1. Definition 3.2 (General Multilevel estimator). Let R 2 and let Y (j),k h,n
k 1 be an i.i.d. sequence of copies of Y (j) h,n . A Multilevel estimator of depth R attached to a stratification strategy q = (q 1 , . . . , q R ) with q j > 0, j = 1, . . . , R, and j q j = 1 and an allocation matrix T, is defined for every integer N 1 and h ∈ H by

Ȳ N,q h,n = R j=1 1 N j N j k=1 T j , Y (j),k h,n (18)
where for all j ∈ 1, . . . , R ,

N j = ⌈q j N ⌉ (allocated budget to compute E T j , Y (j)
h,n ).

• If furthermore the R-level allocation matrix T satisfies

T 1 = e 1 and R j=1 T j = e R ,
the estimator is called a Multilevel Monte Carlo (MLMC) estimator of order R.

• If, furthermore, the R-level allocation matrix T satisfies

T 1 = e 1 and R j=1
T j = w, where w is the unique solution to (10), the estimator is called a Multilevel Richardson-Romberg (ML2R) estimator of order R.

Remark 3.3.

• Note that the assumption T 1 = e 1 is not really necessary. It simply allows for more concise formulas in what follows.

• In this framework, denoting by 0 the null column vector of R R , the crude Monte Carlo is associated to the allocation matrix T = (e 1 , 0, . . . , 0) an the multistep Richardson-Romberg estimator is associated to T = (w, 0, . . . , 0).

Within the abstract framework of a parametrized Monte Carlo simulation described in Section 2.1, the structure parameter π of the multilevel estimators ( Ȳ N,q h,n ) N 1 defined by ( 18) is π = (π 0 , q) where q = (q 1 , . . . , q R ) ∈ (0, 1) R ,

i q i = 1, π 0 = (h, n 1 , . . . , n R , R, T) ∈ Π 0 .

Cost, complexity and effort of a Multilevel estimator

In order to optimize (minimize) the effort φ(π) of the estimator (18), let us evaluate its unitary computational complexity. For a simulation size N , the numerical cost induced by the estimators

Y N,q h,n , N 1, reads Cost( Ȳ N,q h,n ) = R j=1 N j R i=1 1 h n i 1 {T j i =0} = N κ(π) (19)
where κ(π) is the unitary complexity given by

κ(π) = 1 h R j=1 q j R i=1 n i 1 {T j i =0} . (20)
However, it may happen, like for nested Monte Carlo (see Section 4.2), that the internal consistency of the family Y h leads to reduce the computational cost since the computational complexity entirely results from the most refined "scheme". If so,

κ(π) = 1 h R j=1 q j max 1 i R n i 1 {T j i =0} .
It follows that the effort of such a Multilevel estimator is given by

φ(π) = ν(π) κ(π) =   R j=1 1 q j var T j , Y (j) h,n   κ(π). (21)

Bias error of a Multilevel estimator

We now establish the bias error in this general framework. The following bias error result follows straightforwardly from the weak error decomposition (W E α, R) and the definition of an allocation matrix T.

Proposition 3.4. Assume (W E α, R).

(a) ML2R estimator: Let R ∈ 2, . . . , R be the depth of an ML2R estimator. For any admissible stratification strategy q = (q 1 , . . . , q R ), the bias error reads

µ(π 0 , q) = (-1) R-1 c R h R n! α 1 + η R,n (h) ( 22
)
where η R,n (h) = (-1) R-1 n! α R r=1 w r n αR r η R h n r (see (14)) with η R defined in (W E α, R).
(b) MLMC estimator: Let R 2 be the depth of an MLMC estimator. For any admissible stratification strategy q = (q 1 , . . . , q R ), the bias error reads

µ(π 0 , q) = c 1 h n R α 1 + η 1 h n R (23) with η 1 defined in (W E α, R).

Toward the optimal parameters

The optimization problem (8) is not attainable directly, so we decompose it into two successive steps:

Step 1: Minimization of the effort φ over all stratification strategies q = (q j ) 1 j R (as a function of a fixed bias parameter h). In practice, we will optimize an upper-bound φ of the true problem

q * = argmin q∈S + (R)
φ(π 0 , q), where φ(π) φ(π), and

φ * (π 0 ) = φ(π 0 , q * ). ( 24
)
This phase is solved in Theorem 3.6 below (an explicit expression for φ is provided in ( 27)). The quantity φ * (π 0 ) is called the optimal stratified effort (with a slight abuse of terminology since φ is only an upper bound of φ).

Step 2: Minimization of the resulting cost as a function of the remaining parameters π 0 for a prescribed L 2 -error ε > 0 (and specification of the resulting size of the simulation and its cost):

π 0 (ε) = argmin π 0 ∈Π 0 |µ(π 0 ,q * )|<ε φ * (π 0 ) ε 2 -µ 2 (π 0 , q * ) , N (π 0 (ε)) = φ * (π 0 (ε)) κ(π 0 (ε), q * )(ε 2 -µ 2 (π 0 , q * ))
.

This second phase is solved asymptotically when ε goes to 0 in Theorem 3.8 and Proposition 3.11, with some closed forms for some h * and R * as functions of ε and of the structural parameters coming from assumptions (W E α, R) and (SE β ).

Optimally stratified effort (Step 1)

Throughout our investigations on these estimators, we will make extensive use in what follows of the following lemma which is a straightforward consequence of Schwarz's Inequality.

Lemma 3.5. For all j ∈ 1, . . . , R , let a j > 0, b j > 0 and q j > 0 such that

R j=1 q j = 1. Then   R j=1 a j q j     R j=1 b j q j     R j=1 a j b j   2
and equality holds if and only if

q j = µ a j b -1 j , j = 1, . . . , R, with µ = R k=1 a k b -1 k -1
.

Theorem 3.6. Assume (SE β ) holds and let θ be defined by (9). Then, the optimally stratified effort φ * defined by (24) satisfies

φ * (π 0 ) φ(π 0 , q * ) = var(Y 0 ) h   1 + θh β 2 R j=1 R i=1 T j i n -β 2 i R i=1 n i 1 T j i =0 1 2   2
where q * = q * (π 0 ) is an optimal strategy given by

       q * 1 (π 0 ) = µ * (1 + θh β 2 ) q * j (π 0 ) = µ * θh β 2 R i=1 T j i n -β 2 i R i=1 n i 1 T j i =0 -1 2 , j = 2, . . . , R, (25)
and µ * is the normalizing constant such that R j=1 q * j = 1.

Proof. Under assumption (17), we have T 1 , Y

h,n = Y

(1) h

and, for every j ∈ {2, . . . , R}, T j , Y

h,n = T j , Y (j) h,n -Y (j) (j)
0 1 since T j , 1 = 0. Hence, using Minkowski's inequality and the strong error assumption, we obtain

∀j 2, var T j , Y (j) h,n R i=1 T j i Y (j) h n i -Y (j) 0 2 2 V 1 h β R i=1 T j i n -β 2 i 2 .
The variance of the Multilevel estimator is then

var Ȳ N,q h,n 1 N   var Y (1) h q 1 + V 1 h β R j=2 1 q j R i=1 T j i n -β 2 i 2   . ( 26
)
On the other hand we have, var

Y (1) h = var (Y h ) E [Y h -E [Y 0 ]] 2 Y h -Y 0 2 2 + 2E [(Y h -Y 0 )(Y 0 -E [Y 0 ])] + var (Y 0 ) var(Y 0 ) + V 1 h β + 2 V 1 h β/2 var Y 0 = var(Y 0 )(1 + θh β 2 ) 2 .
Combining (20), the above inequality (26) and the above upper-bound for var Y

(1) h , we derive the following upper bound for the effort φ(π) φ(π) with

φ(π) = var(Y 0 ) h   (1 + θh β 2 ) 2 q 1 + θ 2 h β R j=2 1 q j R i=1 T j i n -β 2 i 2     R i,j=1
q j n i 1

T j i =0   . (27)
Applying Lemma 3.5 with a 1 = (1+θh

β 2 ) 2 , b 1 = 1 and a j = θ 2 h β R i=1 T j i n -β 2 i 2 , b j = R i=1 n i 1 T j i =0
, j ∈ 2, . . . , R completes the proof.

Remark 3.7 (About variance minimization). We proved in the above proof that for every stratification strategy q = (q 1 , . . . , q R ),

var Ȳ N,q h,n var(Y 0 ) N   (1 + θh β 2 ) 2 q 1 + θ 2 h β R j=2 1 q j R i=1 T j i n -β 2 i 2   .
Then, applying Lemma 3.5 with

a 1 = (1 + θh β 2 ) 2 , b 1 = 1 and a j = θ 2 h β R i=1 T j i n -β 2 i 2 , b j = 1, j ∈ 2, . . . , R , we obtain (since R j=1 q j b j = 1) inf q∈S + (R) var Ȳ N,q h,n var(Y 0 )   1 + θh β 2 R j=1 R i=1 T j i n -β 2 i   2
with an optimal choice (to minimize the variance):

q † 1 = µ † (1 + θh β 2 ), q † j = µ † θh β 2 R i=1 T j i n -β 2 i (µ †
normalizing constant such that n j=1 q † j = 1). Note that this choice is non-optimal when dealing with the effort optimization approach.

Resulting cost optimization (

Step 2)

Bias parameter optimization (first approach)

In this first approach, we fix the depth R 2, the allocation matrix T and the refiners n 1 , . . . , n R and we only optimize the bias parameter h ∈ H with respect to ε > 0, so that

π 0 (ε) = h(ε, n 1 , . . . , n R , R, T).
We recall that φ * (h) φ(h, q * ) =: φ * (h) where

φ * (h) = var(Y 0 ) h   1 + θh β 2 R j=1 R i=1 T j i n -β 2 i R i=1 n i 1 T j i =0 1 2   2 . ( 28
)
Theorem 3.8 (Bias parameter optimization). Assume (W E α, R) and (SE β ). Let R 2 and let n i , i = 1, . . . , R ,be fixed refiners.

(a) ML2R estimator: Assume R ∈ 2, . . . , R is such that c R = 0. An ML2R estimator of depth R satisfies inf h∈H |µ(h,q * )|<ε Cost Ȳ N,q * h,n ∼ (1 + 2αR) 1+ 1 2αR 2αR |c R | 1 αR var(Y 0 ) n! 1 R ε 2+ 1 αR as ε → 0
with q * defined in (25). This asymptotically optimal bound is achieved with a bias parameter given by

h * (ε, R) = (1 + 2αR) -1 2αR ε |c R | 1 αR n! 1 R . ( 29
) (b) MLMC estimator: Assume c 1 = 0. An MLMC estimator of depth R satisfies inf h∈H |µ(h,q * )|<ε Cost Ȳ N,q * h,n ∼ (1 + 2α) 1+ 1 2α 2α |c 1 | 1 α var(Y 0 ) n R ε 2+ 1 α as ε → 0
with q * defined in (25). This done asymptotically optimal bound is achieved with a bias parameter given by

h * (ε, R) = (1 + 2α) -1 2α ε |c 1 | 1 α n R . (30)
Proof. (a) By definition of the effort φ and the bias µ of the estimator we have (see Section (2.1))

Cost Ȳ N,q * h,n = φ * (h) ε 2 -µ 2 (h, q * ) .
It follows from (28) that the cost minimization problem is upper-bounded by the more tractable problem

inf h∈H, |µ(h,q * )|<ε h φ * (h) h(ε 2 -µ 2 (h, q * ))
with a bias µ(h, q * ) satisfying (22). First note that lim h→0 h φ(h, q * ) = var(Y 0 ). We will consider now the denominator h(ε 2µ 2 (h, q * )). Elementary computations show that, for fixed real numbers a, R ′ > 0, the function g a,R ′ defined by g a,R

′ (ξ) = ξ(1 -a 2 ξ 2R ′ ), ξ > 0, satisfies ξ(a, R ′ ) := argmax ξ>0 g a,R ′ (ξ) = (2R ′ + 1) 1 2 a -1 R ′ and max (0,+∞) g a,R ′ = 2R ′ (2R ′ + 1) 1+ 1 2R ′ a -1 R ′ . Then, set R ′ = Rα, ã = | w R+1 c R | ε = |c R | |n! α ε .
Inspired by what precedes, we make the sub-optimal choice

h(ε) = h(ε, R, α) = ξ ã, αR = ε (2αR + 1) 1 2 |c R | 1 αR n! 1 R corresponding to the case η R,n ≡ 0. It is
clear that, for small enough ε, µ 2 (h, q * ) < ε 2 which makes this choice admissible. Hence

inf h∈H, |µ(h,q * )|<ε φ * (h) ε 2 -µ 2 (h, q * ) 1+ 1 2αR (2αR+1) 1 2αR |c R | 1 αR h(ε) φ * (h(ε)) n! 1 R ε 2+ 1 αR 1 1 - (η R,n (h(ε))+1) 2 -1 2αR . (31)
The "limsup" side of the result follows since lim h→0 η R,n (h) = 0. On the other hand, it follows from the definition (21) of the effort φ that

φ * (h) = 1 h   R j=1 1 q * j var T j , Y (j) h,n     R i,j=1
q * j n i 1

T j i =0   .
Then Schwarz's Inequality implies

φ * (h) 1 h   R j=1 var T j , Y (j) h,n R i=1 n i 1 {T j i =0}   2 1 h max 1 j R var T j , Y (j) h,n R i=1 n i 1 {T j i =0} 1 h max 1 j R var T j , Y (j) h,n since n i n 1 = 1, i = 1, . . . , R. Denoting g(h) = max 1 j R var T j , Y (j) h,n
one clearly has lim h→0 g(h) = var(Y 0 ) under the strong assumption (SE β ) and, as a consequence, lim h→0 h φ(h) = var(Y 0 ). Hence, the cost minimization problem is lower bounded by the more explicit problem

inf h∈H |µ(h,q * )|<ε g(h) h(ε 2 -µ 2 (h, q * ))
.

Let η ∈ 0, 1). There exists ε η > 0 such that, for every h ∈ (0, h(ε η )),

|g(h) -var(Y 0 )| η var(Y 0 ) and |η R,n (h)| η.
We derive from Equation (22) that

µ(h(ε η ), q * ) 2 ε 2 η (1 -η) 2 2αR + 1 . Consequently, if ε < εη(1-η) √ 2αR+1 , for every h ∈ (0, h(ε η )), µ(h, q * ) 2 µ(h(ε η ), q *
) 2 and, as soon as

µ 2 (h, q * ) < ε 2 , one has | µ(h, q * )| c R h R n! α (1 -η), so that g(h) h(ε 2 -µ(h, q * ) 2 ) var(Y 0 )(1 -η) h ε 2 -(1 -η) 2 (c R /n! α ) 2 h 2αR .
Taking advantage of what was done in the "lim sup" part, we get

inf h∈H µ(h,q * )<ε g(h) h(ε 2 -µ(h, q * ) 2 ) 1 + 1 2αR (2αR + 1) 1 2αR |c R | 1 αR var(Y 0 ) n! 1 R ε 2+ 1 αR (1 -η) 1+ 1 αR .
Letting ε and η successively go to zero, yields the "lim inf" side.

(ii) Owing to (23), the bias µ(h, q) is now given by

µ(h, q) = h n R st α c 1 + η 1 h n R with lim h→0 η 1 (h) = 0.
Following the lines of the proof of (i) with R ′ = α completes the proof.

Remark 3.9.

• The fact that the function lim h→0 hφ * (h) = var(Y 0 ) follows from the L 2 -strong convergence of Y h toward Y 0 . Its rate of convergence plays no explicit rôle in this asymptotic rate of the cost as ε → 0. However, this strong rate is important to design a practical stratification among the R independent Brownian motions, which is the key to avoid an explosion of this term.

• When c R = 0 the same reasoning can be carried out by considering any small parameter ǫ R 0 > 0. Anyway in practice c R is usual not known and the impact of this situation is briefly discussed further on in Section 3.3.3.

• When c 1 = 0, specific weights can be computed (see Practitioner's corner in Section 5.1).

Remark 3.10. The asymptotic number N of simulations given by (8) satisfies

N (ε) ∼ 1 + 1 2αR var(Y 0 ) ε 2   R j=1 q * j R i=1 n i 1 T j i =0   -1
as ε → 0

for an ML2R estimator and

N (ε) ∼ 1 + 1 2α var(Y 0 ) ε 2   R j=1 q * j R i=1 n i 1 T j i =0   -1
as ε → 0

for an MLMC estimator.

Templates for R-level allocation matrix T

We now specify the allocation matrix T.

MLMC estimator

The standard Multilevel Monte Carlo allocation matrix used by [START_REF] Heinrich | Multilevel Monte Carlo methods[END_REF][START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF] is derived from the telescopic summation

E Y h n R = E Y h + R j=2 E Y h n j -Y h n j-1 ,
This corresponds in our general framework to the allocation matrix T defined by T j = e je j-1 for j ∈ 2, . . . , R i.e.

T =           1 -1 0 • • • • • • 0 0 1 -1 0 • • • 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 • • • • • • 0 1 -1 0 • • • • • • • • • 0 1           . (MLMC)
In this particular case, the resulting upper-bound φ * of

φ * writes φ * (π 0 ) = var(Y 0 ) h   1 + θh β 2 R j=1 n -β 2 j-1 + n -β 2 j n j-1 + n j   2 (32)
with the convention n 0 = (n 0 ) -1 = 0.

ML2R estimator

The corresponding allocation matrix T for the ML2R estimator is defined by

T j = -W j e j-1 + W j e j for j ∈ 2, . . . , R with W j = R k=j
w k , w is given by (10) i.e.

T =           1 -W 2 0 • • • • • • 0 0 W 2 -W 3 0 • • • 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 • • • • • • 0 W R-1 -W R 0 • • • • • • • • • 0 W R           . ( ML2R
)
The resulting upper-bound φ * now reads (still with the convention

n 0 = (n 0 ) -1 = 0). φ * (π 0 ) = var(Y 0 ) h   1 + θh β 2 R j=1 W j n -β 2 j-1 + n -β 2 j n j-1 + n j   2 , (33)
In the sequel, we will focus on the above choice (ML2R) for the allocation matrix T which leads to the ML2R estimator (3) proposed in the introduction. With this allocation matrix (ML2R) the ML2R estimator writes as a weighted version of MLMC

Ȳ N,q h,n = 1 N 1 N 1 k=1 Y (1),k h + R j=2 W j N j N j k=1 Y (j),k h n j -Y (j),k h n j-1
where N j = ⌈q j N ⌉. Alternative choices for T are proposed in Section 5.1.

Bias parameter and depth R optimization (second approach) for geometric refiners

In this second approach, we consider geometric refiners with root M 2 of the form

n i = M i-1 , i = 1, . . . , R.
These are the refiners already considered in regular multilevel Monte Carlo framework.

π 0 (ε) = h(ε, M, R(ε), T), R(ε, M, T) . Theorem 3.11. Assume (SE β ) holds for β > 0. (a) ML2R estimator: Assume (W E α,∞ ), sup R∈N sup h ′ ∈(0,h) |η R (h ′ )| < +∞ for every h ∈ H and lim R→+∞ |c R | 1 R = c ∈ (0, +∞). The ML2R estimator with allocation matrix T in (ML2R) satisfies lim sup ε→0 v(β, ε) × inf h∈H,R 2 |µ(h,R,q * )|<ε Cost Ȳ N,q h,n K(α, β, M ) (34) with v(β, ε) =        ε 2 (log(1/ε)) -1 if β = 1, ε 2 if β > 1, ε 2 e -1-β √ α √ 2 log(1/ε) log(M ) if β < 1.
These bounds are achieved with an order

R * (ε) =     1 2 + log c 1 α h log(M ) + 1 2 + log c 1 α h log(M ) 2 + 2 log(A/ε) α log(M )     , A = √ 1 + 4α
satisfying lim ε→0 R * (ε) = +∞ and a bias parameter h * = h * (ε, R(ε)) given by (29). The finite real constant K(α, β, M ) depends on M and on the structural parameters α, β, V 1 , var(Y 0 ), h, namely

K(α, β, M ) =                2V 1 α Wα(M )M (1+M )(1+M -1 2 ) 2 log(M ) if β = 1, var(Y 0 )M h 1 + θ h β 2 Wα(M )M β-1 2 √ 1+M (1+M -β 2 ) 1-M 1-β 2 2 if β > 1, V 1 h 1-β c (1-β) α W 2 α (M )M (1+M )(1+M -β 2 ) 2 (M 1-β 2 -1) 2 if β < 1. (35)
(b) MLMC estimator: Assume (W E α,1 ) and c 1 = 0. The MLMC estimator (with allocation matrix T defined in (MLMC)) satisfies

lim sup ε→0 v(β, ε) × inf h∈H,R 2 |µ(h,R,q * )|<ε Cost Ȳ N,q h,n K(α, β, M ) (36) with v(β, ε) =      ε 2 (log(1/ε)) -2 if β = 1, ε 2 if β > 1, ε 2+ 1-β α if β < 1.
These bounds are achieved with an order

R * (ε) =     1 + log c 1 1 α h log(M ) + log(A/ε) α log(M )     , A = √ 1 + 2α
satisfying lim ε→0 R * (ε) = +∞ and a bias parameter h * = h * (ε, R(ε)) given by (30). The finite real constant K(α, β, M ) depends on M and the structural parameters α, β, V 1 , var(Y 0 ), h, namely

K(α, β, M ) =                1 + 1 2α V 1 α 2 M (1+M )(1+M -1 2 ) 2 log(M ) 2 if β = 1, 1 + 1 2α var(Y 0 )M h 1 + θ h β 2 M β-1 2 √ 1+M (1+M -β 2 ) 1-M 1-β 2 2 if β > 1, (1+2α) 1+ 1-β 2α 2α V 1 h 1-β |c 1 | (1-β) α M (1+M )(1+M -β 2 ) 2 (M 1-β 2 -1) 2 if β < 1.
Remark 3.12.

• It is proved in Appendix B that lim M →+∞ W α (M ) = 1 and, to be more precise, that W α (M ) -1 ∼ M -α as M → +∞.

• The assumption on the functions η R and the sequence (c R ) R 2 in (a) of the above proposition are reasonable though probably impossible to check in practice. In particular, note that as soon as the sequence (c R ) R 2 has at most a polynomial growth as a function of R, it satisfies the assumption since c = 1.

• When c = 0, one can replace c R in the proof below by ǫ R 0 (see also Remark 3.9) and carry on the computations (with c = ǫ 0 ). This constant has an impact when β < 1: when ǫ 0 → 0, K(α, β, M ) goes to 0 which emphasizes that we are not in the right asymptotics.

• If β = 1, ML2R is asymptotically more efficient than MLMC by a factor log 1/ε) → +∞ as ε → 0. When β < 1, ML2R (with M = 2) is asymptotically infinitely more efficient than MLMC by a factor ε

-1-β √ α e -1-β α √
2 log(M ) log(1/ε) which goes to +∞ as ε → 0 in a very steep way. To be precise the ratio is greater than 1 as soon as

ε 2 -2 α .
It is clear that it is for this setting that ML2R is the most powerful compared to regular MLMC. When β > 1, both estimators achieve the same rate ε -2 as a virtual unbiased Monte Carlo method based on the direct simulation of Y 0 .

Proof. We provide a detailed poof of claim (a), that of (b) following the same lines.

Step 1: We start from Equation (31) in the proof of Theorem 3.8 which reads

inf h∈H |µ(h,q * )|<ε Cost Ȳ N,q * h,n 1 + 1 2αR φ * (h * (ε)) ε 2 1 1 - (η R,n (h * (ε))+1) 2 -1 2αR with φ * (h * (ε)) = 1 h * (ε) var(Y 0 )   1 + θh * (ε) β 2 R j=1 W j n -β 2 j-1 + n -β 2 j n j-1 + n j   2 ,
(convention n 0 = (n 0 ) -1 = 0). The idea is to choose R = R * (ε) as large as possible provided the optimal bias parameter h * ∈ H. The form of the refiners

n i = M i-1 implies that n! = M R(R-1) 2 so that h * (ε, R) = (1 + 2αR) -1 2αR |c R | -1 αR ε 1 αR M R-1 2 .
Note that, by assumption we have lim R→+∞ |c R | -1 αR = c -1 α with c ∈ (0, +∞). We choose to try to saturate the constraint h * h which leads to impose formally (for R 2)

h * (ε, R) (1 + 4α) -1 2αR c -1 α ε 1 αR M R-1 2 = h,
(where we temporarily forget that R should be an integer). Let A = c 1 α h. As a consequence, we search for the positive zero R + (ε) of the polynomial

P (R) = R(R -1) 2 log(M ) -R log(A) - 1 α log √ 1 + 4α/ε , that is R + (ε) = 1 2 + log(A) log(M ) + 1 2 + log(A) log(M ) 2 + 2 log( √ 1+4α/ε) α log M and denoting R * (ε) = ⌊R + (ε)⌋ we obtain P (R * (ε)) 0. Hence, h * (ε, R * (ε)) = he P (R * ) R * - P (R + ) R +
h. Let us show that our choice h * (ε, R * (ε)) for the bias parameter (see (29)) is admissible -i.e. µ(h * (ε, R * (ε)), R * (ε), q * ) 2 < ε 2 -at least for small enough ε. Elementary computations show that

µ h * (ε, R * (ε)) 2 = (c R * (ε) w R * (ε)+1 ) 2 (h * (ε, R * (ε))) 2αR * (ε) = ε 2 e -αR * (ε) log(1+2αR * (ε)) 1 + η R * (ε),n h * (ε, R * (ε)) 2 .
Our choice for R * (ε) implies that h * (ε, R * (ε)) is upper-bounded by h. Claim 6 of Proposition A.2 in Appendix A and the assumption on η R imply that, sup

0<h ′ <h |η R * (ε),n (h ′ )| B α (M ) sup h ′ ∈(0,h) |η R * (ε) (h ′ )| B α (M ) sup R 1 sup h ′ ∈(0,h) |η R (h ′ )| < +∞.
As a consequence of the assumption made on the functions η R , it is clear that

µ h * (ε, R * (ε)), R * (ε), q * 2 = o(ε 2 ) since R * (ε) → +∞ as ε → 0.
Hence, our choice for the bias parameter is admissible at least for small enough ε.

Likewise, the assumption on the functions η R implies that lim

ε→0 η R * (ε),n (h(ε,R * (ε)))+1 2 -1 2αR * (ε) = 0. We have then proved that lim sup ε→0   l(ε, R * (ε)) inf h∈H |µ(h,R,q * )|<ε × Cost Ȳ N,q * h,n   M var(Y 0 ) h with l(ε, R) = ε 2   1 + θ h * (ε, R) β 2 R j=1 W j n -β 2 j-1 + n -β 2 j n j-1 + n j   -2
.

It follows from Claim 5 of Proposition A.2 in Appendix A that max j=1,...,R W i W α (M ). On the other hand, standard computations show that, for every j = 2, . . . , R, n

-β 2 j-1 + n -β 2 j n j-1 + n j = M β-1 M j 1-β 2 1 + M -β 2 1 + M 1 2 .
(37)

Moreover, with our convention on n 0 , it still holds true as an inequality ( ) for j = 1. So

l(ε, R) ε 2   1 + θh * (ε, R) β 2 W α (M )M β-1 √ 1 + M (1 + M -β 2 ) R j=1 M j 1-β 2   -2
.

Step 2: Now we will inspect successively the three cases depending on the strong rate convergence parameter β > 0.

Case β = 1. In that case,

l(ε, R * (ε)) ε 2 1 + θh * (ε, R * (ε)) β 2 W α (M ) √ 1 + M (1 + M -1 2 )R * (ε) -2 , ε 2 1 + θ h β 2 W α (M ) √ 1 + M (1 + M -1 2 )R + (ε) -2
, and, as

R 2 + (ε) ∼ 2 α log(M ) log(1/ε) as ε → 0, we get (34) with K(α, 1, M ) given by (35) keeping in mind that V 1 = var(Y 0 )θ 2 . Case β > 1. Noting that R j=1 M j 1-β 2 M 1-β 2 1-M 1-β 2
, we get

l(ε, R * (ε)) ε 2 1 + θ h β 2 W α (M )M β-1 2 √ 1 + M (1 + M -β 2 ) 1 -M 1-β 2 -2
, which yields (34) with K(α, β, M ) given by (35).

Case β < 1. In that setting, we note this time that R j=1 M j 1-β 2

M (R+1) 1-β 2 M 1-β 2 -1 so that l(ε, R * (ε)) ε 2 1 + θh β 2 W α (M ) √ 1 + M (1 + M -β 2 ) M 1-β 2 -1 M (R + (ε)-1) 1-β 2 -2 . As R + (ε) satisfies h * (ε, R + (ε)) = h, we obtain M R + (ε)-1 2 = h c 1 α ε - 1 αR + (ε) . We have ε - 1 αR + (ε) ∼ e log(M ) 2α
log(1/ε) as ε → 0. Elementary, although tedious computations yield (34) with K(α, β, M ) given by ( 35).

(b) The choice for R * (ε) follows from the formal constraint lim sup

ε→0 h * (ε, R * (ε)) = (1 + 2α) -1 2α |c 1 | -1 α ε 1 α M R * (ε)-1 = h.
Then, the proof follows the same lines as that of (a).

Remark 3.13 (On the constraint h). In the proof we choose to saturate the constraint h * h. If we consider h * = χ where χ is a free parameter in (0, h], then the asymptotic constants K(α, β, M ) for the renormalized optimized cost in Theorem 3.11 depends on χ and one verifies the following facts:

• When β < 1, one can write K(α, β, M, χ) = χ 1-β K(α, β, M, 1) which this time suggests to start the simulation with a small upper bias parameter χ < h.

• When β > 1, the asymptotic cost of the simulation increases in ε 2 like a (virtual) unbiased one.

In that very case, it appears that the asymptotic constant K(α, β, M, χ) can itself be optimized as a function of χ. Namely, if we set

κ 1 = var(Y 0 )M χ and κ 2 = θ 2 W α (M ) 2 M β-1 (1 + M )(1 + M -β ) (1 -M 1-β 2 ) 2 , then χ opt = β -2 β+1 κ -1 β+1 2 and K(α, β, M, χ opt ) = (β + 1) 2 β -2 β+1 κ 1 κ 1 β+1 2 .
• When β = 1, the asymptotic constant K(α, β, M, χ) does not depend on χ. This suggests that the choice of the upper bias parameter is not decisive, at least for high accuracy computations (ε close to 0). The choice χ = h remains the most natural.

Examples of applications 4.1 Brownian diffusion approximation

Euler scheme In fact, the (one-step) Richardson-Romberg extrapolation is well-known as an efficient mean to reduce the time discretization error induced by the use of an Euler scheme to simulate a Brownian diffusion. In this field of Numerical Probability, its introduction goes back to Talay and Tubaro in their seminal paper [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF] on weak error expansion, followed by the case of non smooth functions in [START_REF] Bally | The law of the Euler scheme for stochastic differential equations. I. Convergence rate of the distribution function[END_REF] under an Hörmander hypo-ellipticity assumption. It relies on the following theorem.

Theorem 4.1. Let b : R d → R d , σ : R d → M(d, q) and let (W t ) t 0 be a q-dimensional standard Brownian motion defined on a probability space (Ω, A, P). Let X = (X t ) t∈[0,T ] be a diffusion process, strong solution to the Stochastic Differential Equation (SDE)

dX t = b(X t )dt + σ(X t )dW t , t ∈ [0, T ] , X 0 = x 0 ∈ R d , ( 38
)
and its continuous Euler scheme Xh = ( Xh t ) t∈[0,T ] with bias (step) parameter h = T /n defined by

Xh t = X 0 + t 0 b Xh s ds + t 0 σ Xh s dW s ,
where s = kh on [kh, (k + 1)h) .

(a) Smooth setting (Talay-Tubaro [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF]): If b and σ are infinitely differentiable with bounded partial derivatives and if f : R d → R is an infinitely differentiable function, whose all partial derivatives have polynomial growth, then for a fixed T > 0 and every integer

R ∈ N * E f ( Xh T ) -E [f (X T )] = R k=1 c k h k + O h R+1 , ( 39
)
where the coefficients c k depend on b, σ, f , T (but not on h).

(b) (Hypo-)Elliptic setting (Bally-Talay [START_REF] Bally | The law of the Euler scheme for stochastic differential equations. I. Convergence rate of the distribution function[END_REF]): If b and σ are infinitely differentiable with bounded partial derivatives and if σ is uniformly elliptic in the sense that

∀ x ∈ R d , σσ * (x) ε 0 I q , ε 0 > 0
or, more generally, if (b, σ) satisfies the strong Hörmander hypo-ellipticity assumption, then (39) holds true for every bounded Borel function f : R d → R.

Other results based on the direct expansion of the density of the Euler scheme allow to deal with a drift b with linear growth (see [START_REF] Konakov | Edgeworth type expansions for Euler schemes for stochastic differential equations[END_REF], in a uniformly elliptic setting, see also [START_REF] Guyon | Euler scheme and tempered distributions[END_REF] at order 1 in a tempered distribution framework). It is commonly shared by the "weak error community", relying on an analogy with recent results on the existence of smooth density from the diffusion, that if the hypo-ellipticity assumption is satisfied except at finitely many points that are never visited by the diffusion, then the claim (b) remains true. The boundedness assumption on σ is more technical than a key assumption. For a recent review on weak error, we refer to [START_REF] Jourdain | A review of recent results on approximation of solutions of stochastic differential equations[END_REF].

To deal with our abstract multilevel framework we consider for a fixed horizon T > 0, the family of Euler schemes Xh with step h ∈ H = { T n , n 1}. We set Y h = f ( Xh T ) and Y 0 = f (X T ) for a smooth enough function f with polynomial growth. The above theorem says that condition (W E α, R) is satisfied with R = +∞ and α = 1. However, for a fixed R, the differentiability assumption on b, σ and f can be relaxed by simply assuming that these three functions are

C R+5 b on [0, T ] × R d .
On the other hand, as soon as f : R d → R is Lipschitz continuous, it is classical background that (SE β ) is satisfied with β = 1 as an easy consequence of the fact that the (continuous) Euler scheme Xh converges for the sup-norm toward X in L 2 (in fact in every L p -space) at rate √ h as the step h goes to 0.

In such a setting, we can implement multilevel estimators with α = β = 1.

Milstein scheme

The Milstein scheme is a second order scheme which satisfies (SE β ) with β = 2 and (W E α, R) still with α = 1 (like the Euler scheme). Consequently, provided it can be implemented, the resulting multilevel estimators should be designed with these parameters. However, the main drawback of the Milstein scheme when the SDE is driven by a multidimensional Brownian motion (q 2), it requires the simulation of Lévy areas, for which there is no known efficient method (except in dimension 2). In a recent work [START_REF] Giles | Antithetic multilevel Monte Carlo estimation for multi-dimensional SDEs without Lévy area simulation[END_REF], Giles and Szpruch introduce a suitable antithetic multilevel correction estimator which avoids the simulation of these Lévy areas. This approach could be easily combined with our weighted version of MLMC.

Note that in the β > 1 case, Rhee and Glynn introduce in [START_REF] Rhee | A new approach to unbiased estimation for SDE's[END_REF] a class of finite-variance optimally randomized multilevel estimators which are unbiased with a square root convergence rate.

Path-dependent functionals When a functional

F : C([0, T ], R d ) → R is Lipschitz continuous for the sup-norm, it is straightforward that F ( Xh ) and F (X) satisfy (SE β ), β = 1, with H = { T n , n 1}
, (but this is no longer true if one considers the stepwise constant Euler scheme since the rate of convergence is then log n/n ≍ √ -h log h). More generally, if F is β-Hölder, β ∈ (0, 1], then this family satisfies (SE β ). High order expansions of the weak error are not available in the general case, however first order expansion have been established for specific functionals like F (w) = f T 0 w(s)ds or F (w) = f (w(T ))1 {τ D (w)>T } where τ D (w) is the exit time of a domain D of R d which show that they satisfy (W E α, R) with α = 1 and R = 1 (see e.g. [START_REF] Lapeyre | Competitive Monte Carlo methods for the pricing of Asian Options[END_REF][START_REF] Gobet | Weak approximation of killed diffusion using Euler schemes[END_REF]). More recently, new results on first order weak error expansions have been obtained for functionals of the form F (w) = f w(T ), sup t∈[0,T ] w(t) (see [START_REF] Giles | Analysing multi-level Monte Carlo for options with nonglobally Lipschitz payoff[END_REF] and [START_REF] Alfonsi | Pathwise optimal transport bounds between a onedimensional diffusion and its Euler scheme[END_REF]). Thus, for the weak error expansion, it is shown in [START_REF] Alfonsi | Pathwise optimal transport bounds between a onedimensional diffusion and its Euler scheme[END_REF] that, for every η > 0, there exists a real constant C η > 0 such that

E f X T , sup t∈[0,T ] X t -E f Xn T , sup t∈[0,T ] Xn t C η N 2 3 -η .
For a review of recent results on approximation of solutions of SDEs, we again refer to [START_REF] Jourdain | A review of recent results on approximation of solutions of stochastic differential equations[END_REF].

Remark 4.2. Note that, as concerned the MLMC estimator, in the general setting of the discretization of a Brownian diffusion by an Euler scheme, a Central Limit Theorem (with stable weak convergence) has been obtained in [START_REF] Ben Alaya | Central Limit Theorem for the Multilevel Monte Carlo Euler Method and Applications to Asian Options[END_REF]. A similar approach applied to the ML2R estimator with allocation matrix (ML2R) (weighted version of MLMC), should yield a similar Central Limit Theorem.

Nested Monte Carlo

The purpose of the so-called nested Monte Carlo method is to compute by simulation quantities of the form

E f E [X | Y ]
where (X, Y ) is a couple of R × R q Y -valued random variable defined on a probability space (Ω, A, P) with X ∈ L 2 (P) and f : R → R is a Lipschitz continuous function with Lipschitz coefficient [f ] Lip .

Such quantities often appear in financial application, like compound option pricing or risk estimation (see [START_REF] Broadie | Efficient risk estimation via nested sequential simulation[END_REF]) and in actuarial sciences (see [START_REF] Devineau | Construction d'un algorithme d'accélération de la méthode des "simulations dans les simulations" pour le calcul du capital économique Solvabilité II[END_REF]) where nested Monte Carlo is widely implemented.

The idea of replacing conditional expectations by Monte Carlo estimates also appears in [START_REF] Belomestny | Multilevel dual approach for pricing American style derivatives[END_REF] where authors devise a multilevel dual Monte Carlo algorithm for pricing American style derivatives. We make the following more stringent assumption: there exists a Borel function F : R q Z × R q Y → R and a random variable Z : (Ω, A) → R q Z independent of Y such that

X = F (Z, Y ).
Then, if X ∈ L 2 , one has the following representation

E [X | Y ] (ω) = E [F (Z, y)] |y=Y (ω) = R q Z F (z, Y (ω))P Z (dz).
To comply with the multilevel framework, we set

H = {1/K, K 1}, Y 0 = f E [X | Y ] , Y 1 K = f 1 K K k=1 F (Z k , Y )
where (Z k ) k 1 is an i.i.d. sequence of copies of Z defined on (Ω, A, P) and independent of Y (up to an enlargement of the probability space if necessary).

The following proposition shows that the nested Monte Carlo method is eligible for multilevel simulation when f is regular enough with the same parameters as the Euler scheme for Brownian diffusions.

Proposition 4.3. Assume X ∈ L 2R . If f is Lipschitz continuous and 2R times differentiable with f (k) bounded, k = R, . . . , 2R, the nested Monte Carlo satisfies (SE β ) with β = 1 and (W E α, R) with α = 1 and R = R -1.
Remark 4.4. When f is no longer smooth, typically if it is the indicator function of an interval, it is still possible to show that nested Monte Carlo is eligible for multilevel Richardson-Romberg approach e.g. in the more constrained framework developed in [START_REF] Hong | Estimating the mean of a non-linear function of conditional expectation[END_REF][START_REF] Gordy | Nested simulation in portfolio risk measurement[END_REF] where X can be viewed as an additive perturbation of Y . Assuming enough regularity in y on the joint density g N (y, z) of Y and the renormalized perturbation, yields an expansion of the weak error (but seems in a different scale). However, in this work we focus on the regular case.

The proof follows form the two lemmas below.

Lemma 4.5 (Strong approximation). Assume f is Lipschitz continuous. For every h ∈ H, Y 0 -Y h 2 2 [f ] 2 Lip X 2 2 -E [X | Y ] 2 2 h (40) so that (Y h ) h∈H satisfies (SE β ) with β = 1. Proof. Let h = 1/K and set E Y [X] = E [X | Y ] for convenience. By Lipschitz continuity, Y 0 -Y h 2 2 = f E Y [X] -f 1 K K k=1 F (Z k , Y ) 2 2 [f ] 2 Lip E Y [X] - 1 K K k=1 F (Z k , Y ) 2 2 and E Y [X] - 1 K K k=1 F (Z k , Y ) 2 2 = E Y [X] 2 2 + 1 K 2 K k=1 F (Z k , Y ) 2 2 - 2 K K k=1 E E Y [X] F (Z k , Y ) + 2 K 2 1 j<k K E F (Z j , Y )F (Z k , Y ) .
Conditioning on Y gives

E E Y [X] F (Z k , Y ) = E E Y [X] E Y [F (Z k , Y )] = E Y [X] 2 2 and E F (Z j , Y )F (Z k , Y ) = E Y [X]
2 2 for j = k. Plugging these two identities in the first expansion finally yields

Y 0 -Y h 2 2 [f ] 2 Lip K F (Z, Y ) 2 2 -E Y [X] 2 2 .
Lemma 4.6 (Weak error). Let f : R → R be a 2R times differentiable function with f (k) , k = R, . . . , 2R, bounded over the real line. Assume X ∈ L 2R (P). Then there exists c 1 , . . . , c R-1 such that

∀h ∈ H, E [Y h ] = E [Y 0 ] + R-1 r=1 c r h r + O h R . (41)
Consequently

(Y h ) h∈H satisfies (W E α, R) with α = 1 and R = R -1. Proof. Let K 1 and X k = F (Z k , Y ) -E Y [F (Z k , Y )] = X k -E Y (X k ), k = 1, . . . , K.
By the multinomial formula we get

( X 1 + • • • + X K ) k = k 1 +•••+k K =k k! k 1 ! • • • k K ! X k 1 1 • • • X k K K so that, taking conditional expectation given Y , yields E Y ( X 1 + • • • + X K ) k = k! k 1 +•••+k K =k K i=1 E Y X k i k i ! since E Y X k i i = E Y X k i . As E Y X i = 0, we obtain E Y ( X 1 + • • • + X K ) k = k! k 1 +•••+k K =k, k i =1 K i=1 E Y X k i k i ! . Let I = I(k 1 , . . . , k K ) = {i | k i = 0, 1}. It is clear that 1 |I| k/2.
By symmetry, we have now that

k 1 +•••+k K =k, k i =1 K i=1 E Y X k i k i ! = 1 ℓ (k/2)∧K I⊂{1,...,K},|I|=ℓ, i∈I k i =k,k i 2 K i=1 E Y X k i k i ! = 1 ℓ k/2 K ℓ 1 i ℓ k i =k-2ℓ ℓ i=1 E Y X 2+k i (2 + k i )! .
As a consequence, for every integer R 1,

E Y Y h = E Y Y 0 + 2R-1 k=1 f (k) E Y [X] k!K k E Y ( X 1 + • • • + X K ) k + R 2R-1 (Y ) = E Y [Y 0 ] + 2R-1 k=1 f (k) E Y [X] k!K k 1 ℓ (k/2)∧K K ℓ a k,ℓ + R 2R-1 (Y )
where

a k,ℓ = k 1 +•••+k ℓ =k-2ℓ ℓ i=1 E Y X 2+k i (2 + k i )! and |R 2R-1 (Y )| f (2R) sup (2R)! 1 K 2R E Y X 1 + • • • + X K 2R .
By the Marcinkiewicz-Zygmund Inequality we get

|R 2R-1 (Y )| (B M Z 2R ) 2R f (2R) sup (2R)! 1 K 2R E Y | X 2 1 + • • • + X 2 K R f (2R) sup (B M Z 2R ) 2R (2R)! 1 K R E Y X 2R
where

B M Z p = 18 p 3 2 (p-1) 1 2
, p > 1 (see [START_REF] Shiryaev | Probability[END_REF] p.499). Now we write the polynomial x(x -1)

• • • (x -ℓ + 1)
on the canonical basis 1, x, . . . , x n ,. . . as follows

x(x -1) • • • (x -ℓ + 1) = ℓ m=0 b ℓ,m x m (b ℓ,ℓ = 1 and b ℓ,0 = 0).
Hence,

E Y [Y h ] = E Y [Y 0 ] + 2R-1 k=1 k 2 ℓ=1 ℓ m=1 f (k) E Y [X] k! 1 K k-m a k,ℓ b ℓ,m + O K -R where K R O(K -R ) is bounded by a deterministic constant. For every r ∈ {1, . . . , R -1}, set J R,r = (k, l, m) ∈ N 3 , 1 k 2R -1, 1 ℓ k/2, 1 m ℓ, k = m + r
(note that one always has k (2m) ∨ 1 so that km 1 when k, l, m vary in the admissible index set). We finally get

E Y [Y h ] = E Y [Y 0 ] + 2R-1 r=1 (k,ℓ,m)∈J R,r f (k) E Y [X] k! a k,ℓ b ℓ,m 1 K r + O K -R . = E Y [Y 0 ] + R-1 r=1 c r K r + O K -R .
Taking the expectation in the above equality yields the announced result.

Remark 4.7. Though it is not the only term included in the final O(K -R ), it is worth noticing that

(B M Z 2R ) 2R (2R)! 1 R ∼ (36R) 2 2R e -2
∼ 18 e 2 as R → +∞ owing to Stirling's formula. This suggests that, e.g. if all the derivatives of f are uniformly bounded, lim sup

R→+∞ |c R | 1 R < +∞.
5 Numerical experiments

Practitioner's corner

We summarize here the study of the Section 3. We have proved in Theorems 3.6, 3.8 and 3.11 that the asymptotic optimal parameters (as ε goes to 0) R, h, q and N depend on structural parameters α, β, V 1 , c 1 , var(Y 0 ) and h (recall that θ = V 1 / var(Y 0 )). Note we do not have optimized the design of the multilevel estimators, namely the allocation matrix T and the refiners n i , i = 2, . . . , R. We propose in this Section a numerical procedure to choose a good value of M in the case n i = M i-1 .

About structural parameters

Implementing MLMC or ML2R estimator needs to know both the weak and strong rates of convergence of the biased estimator Y h toward Y 0 . The exponents α and β are generally known by a mathematical study of the approximation (see Section 4.1 for Brownian diffusion discretization and Section 4.2 for nested Monte Carlo). The parameter V 1 comes from the strong approximation rate assumption (SE β ) and a natural approximation for

V 1 is V 1 ≃ lim sup h→0 h -β Y h -Y 0 2 2
Since Y 0 cannot be simulated at a reasonable cost, one may proceed as follows to get a good empirical estimator of V 1 . First, assume that, in fact, Y h -Y 0 2 2 ∼ V 1 h β as h → 0 but that this equivalence still holds for not too small step h. Then, one derives from Minkowski's Inequality that, for every integer

M 1, Y h -Y h M 2 Y h -Y 0 2 + Y 0 -Y h M 2 so that V 1 (1 + M -β 2 ) -2 h -β Y h -Y h M 2 2 .
As a consequence, if we choose M = M max large enough (see (46) below), we are led to consider the following estimator

V 1 (h) = 1 + M -β 2 max -2 h -β Y h -Y h Mmax 2 2 . ( 42
)
The estimation of the c i (c 1 for crude Monte Carlo and an MLMC estimators and c = lim

R→∞ |c R | 1 R
for the ML2R estimator) is much more challenging. So these methods are usually implemented in a blind way by considering the coefficient |c R | 1 R of interest equal to 1. Note that, even in a crude Monte Carlo method, these structural parameters are useful (and sometimes necessary) to deal with the bias error (see Proposition 2.2).

Design of the Multilevel

The allocation matrix is fixed by the template (ML2R) for the multilevel Richardson-Romberg estimator and by the template (MLMC) for the multilevel Monte Carlo estimator. Alternative choices could be to consider for the ML2R estimator another allocation matrix T satisfying (17) like T j = w j e 1 + w j e j for j ∈ 2, . . . , R which reads

T =         1 -w 2 -w 3 • • • -w R 0 w 2 0 • • • 0 0 0 w 3 . . . 0 . . . . . . . . . . . . . . . 0 0 • • • 0 w R         . ( 43
)
We could also consider a lower triangular allocation matrix (through it does not satisfy the conventional assumption

T 1 = e 1 ) T =            W 1 0 • • • • • • • • • 0 -W 1 W 2 0 • • • • • • 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 • • • • • • -W R-2 W R-1 0 0 • • • • • • • • • -W R-1 1            where W j = j k=1 w k . (44)
The refiners can be specified by users but it turns out that the parametrized family n i = M i-1 , i = 1, . . . , R (M ∈ N, M 2) seems the best compromise between variance control and implementability. The parameter α being settled, all the related quantities like (W i (R, M )) 1 i M can be tabulated for various values of M and R and can be stored offline.

Taking advantage of c 1 = 0 When c 1 = 0, only R -1 weights are needed to cancel the (remaining) coefficients up to order R i.e. c r , r = 2, . . . , R -1 (instead of R). One easily shows that, if w (R-1) r r=1,...,R-1 denotes the weight vector at order R -1 associated to refiners n 1 = 1 < n 2 , . . . , n R-1 (for a given α), then the weight vector w (R) at order R (with size R -1) reads

w (R) r = n α r w (R-1) r 1 s R-1 n α s w (R-1) s , r = 1, . . . , R -1.

Asymptotic optimal parameters

In the case n i = M i-1 (with the convention n 0 = n -1 0 = 0), we can summarize the asymptotic optimal value of the parameters q, R, h and N in the table 5.1 for the (ML2R) estimator and in the table 5.2 for the (MLMC) estimator.

R     1 2 + log c 1 α h log(M ) + 1 2 + log c 1 α h log(M ) 2 + 2 log (A/ε) α log(M )     , A = √ 1 + 4α h -1 (1 + 2αR) 1 2αR ε -1 αR M -R-1 2 h q q 1 = µ * (1 + θh β 2 ) q j = µ * θh β 2   W j (R, M ) n -β 2 j-1 + n -β 2 j √ n j-1 + n j   , j = 2, . . . , R; 1 j R q j = 1 N 1 + 1 2αR var(Y 0 )   1 + θh β 2 R j=1 W j (R, M ) n -β 2 j-1 + n -β 2 j n j-1 + n j   2 ε 2 R j=1
q j (n j-1 + n j ) Table 5.1: Optimal parameters for the ML2R estimator.

R     1 + log c 1 1 α h log(M ) + log(A/ε) α log(M )     , A = √ 1 + 2α h -1 (1 + 2α) 1 2α ε -1 α c 1 1 α M -(R-1) h q q 1 = µ * (1 + θh β 2 ) q j = µ * θh β 2   n -β 2 j-1 + n -β 2 j √ n j-1 + n j   , j = 2, . . . , R; 1 j R q j = 1 N 1 + 1 2α var(Y 0 )   1 + θh β 2 R j=1 n -β 2 j-1 + n -β 2 j n j-1 + n j   2 ε 2 R j=1 q j (n j-1 + n j )
Table 5.2: Optimal parameters for the MLMC estimator.

as an estimator of var( Ȳ N,q h,n ). The empirical L 2 -error or empirical root mean squared error (RMSE) εL of the estimator used in our numerical experiments is then given by εL = (μ L ) 2 + νL .

(47)

The computations were performed on a computer with 4 multithreaded(16) octo-core processors (Intel(R) Xeon(R) CPU E5-4620 @ 2.20GHz). The code of one estimator runs on a single thread (program in C++11 available on request).

Euler scheme of a geometric Brownian motion

We consider a geometric Brownian motion (S t ) t∈[0,T ] , representative in a Black-Scholes model of the dynamics of a risky asset price between time t = 0 and time t = T :

S t = s 0 e (r-σ 2 2 )t+σWt , t ∈ [0, T ], S 0 = s 0 > 0,
where r denotes the (constant) "riskless" interest rate, σ denotes the volatility and W = (W t ) t∈[0,T ] is a standard Brownian motion defined on a probability space (Ω, A, P). The price or premium of a so-called vanilla option with payoff ϕ is given by e -rT E [ϕ(S T )] and the price of a path dependent option with functional payoff ϕ is given by e -rT E ϕ((S t ) t∈[0,T ] ) . Since (S t ) t∈[0,T ] is solution to the diffusion SDE dS t = S t (rdt + σdW t ), S 0 = s 0 > 0, one can compute the price of an option by a Monte Carlo simulation in which the true process (S t ) t∈[0,T ] is replaced by its Euler scheme ( Skh ) 0 k n , h = T n (even if we are aware that S T is simulatable). The bias parameter set H is then defined by H = T /n, n 1 and h = T .

Although nobody would adopt any kind of Monte Carlo simulation to compute option price in this model in practice since a standard difference method on the Black-Scholes parabolic PDE is much more efficient to evaluate a vanilla option and many path-dependent ones, it turns out that the Black-Scholes model and its Euler scheme is a very demanding benchmark to test and evaluate the performances of Monte Carlo method(s). As a consequence, it is quite appropriate to carry out numerical tests with ML2R (and MLMC).

Vanilla Call option

α = β = 1
The Black-Scholes parameters considered here are s 0 = 100, r = 0.06 and σ = 0.4. The payoff is a European Call with maturity T = 1 year and strike K = 80.

In such a regular diffusion setting (both drift and diffusion coefficients are C ∞ b and the payoff function is Lipschitz continuous), one has α = β = 1. The parameters θ = V 1 / var(Y 0 ) and var(Y 0 ) have been roughly estimated following the procedure (42) on a sample of size 100 000 described in Section 5.1 leading to the values V 1 ≃ 56 and var(Y 0 ) ≃ 876 (so that θ ≃ 0.25). The empirical L 2 -error ǫL is estimated using L = 256 runs of the algorithm and the bias is computed using the true value of the price I 0 = 29.4987.

The results are summarized in Table 5.3 for the ML2R estimator and in Table 5.4 for the MLMC estimator.

As an example, the third line of the Table 5.3 reads as follows: for a prescribed RMSE error ε = 2 -3 = 0.125, the ML2R estimator Ȳ N,q h,n (with allocation matrix (ML2R)) is implemented with the parameters R = 3, h = 1 and refiners n i = 4 i-1 (then n 1 = 1, n 2 = 4 and n 3 = 16) and the sample size N ≃ 319 000. The stratification weights q i (not reported in this Table) are such that the numerical cost Cost( Ȳ N,q h,n ) ≃ 709 800. For such parameters, the empirical RMSE εL ≃ 0.0928 and the computational time of Y N,q h,n ≃ 0.559 seconds. The empirical bias error μL is reported in the 5th column (bias) and the empirical unitary variance νL is reported in the 6th column (var ). Recall that εL = (μ L ) 2 + νL .

Note first that, as expected, the depth parameter R 2 and the numerical cost Cost( Ȳ N,q h,n ) grow slower for ML2R than for MLMC as ε goes to 0. Consequently, regarding the CPU-time for a 0.01 0.1 1 10 100 1000

2 -9 2 -8 2 -7 2 -6 2 -5 2 -4 2 -3 2 -2 2 -1 2 0 CPU-time in seconds Empirical RMSE εL MLMC ML2R 2 -1 2 -2 2 -3 2 -4 2 -5 2 -6 2 -7 2 -8 2 -1 2 -2 2 -3 2 -4 2 -5 2 -6 2 -7 2 -8
Figure 1: Call option in a Black-Scholes model. CPU-time (y-axis, log scale) as a function of εL (x-axis, log 2 scale).

One observes that the L 2 -error of ML2R has a very small bias μL (5th column) due to the particular choice of the weights (W i ) 1 i R .

Figure 2 provides a graphical representation of the performance of both estimators, now as a function of the empirical RMSE ε. It shows that ML2R is faster then MLMC by a factor that goes from 18 up to 48 within the range of our simulations. Additional Figures 6(a that (W E α, R) holds with α = 0.5 and R = +∞. A first computational stage gives us var(Y 0 ) ≃ 303, V 1 ≃ 5.30 and θ ≃ 0.41. The results are summarized in Table 5.7 for ML2R and in Table 5.8 for MLMC. See Figure 3 for a graphical representation. First, note that since β = 0.5, we observe that the function (CPU-time)×ε 2 increases much faster for MLMC than ML2R as ε goes to 0 (see Figure 7 in Appendix C) which agrees with the theoretical asymptotic rates from Theorem 3.11. In fact, in this highly biased example with slow strong convergence, the ratio Cost(MLMC)/ Cost(ML2R) as a function of the prescribed ε = 2 -k goes from 1.1 (k = 2) up to 22 (k = 8), idem for the ratio between CPU-times. When looking at this ratio as a function of the empirical RMSE, it even goes from 3 up to 61 which is huge having in mind that MLMC provides similar gains with respect to a crude Monte Carlo simulation.

k ε = 2 -k L 2 -error time (s) bias var R M h -1 N Cost 1 5.
2 -10 2 -9 2 -8 2 -7 2 -6 2 -5 2 -4 2 -3 2 -2 2 -1 2 0 2 1 CPU-time in seconds Empirical RMSE εL MLMC ML2R 2 -1 2 -2 2 -3 2 -4 2 -5 2 -6 2 -7 2 -8 2 -9 2 -1 2 -2 2 -3 2 -4 2 -5 2 -6 2 -7 2 -8 2 -9
k ε = 2 -k L 2 -error time (s) bias var R M h -1 N Cost 1 5.
k ε = 2 -k L 2 -error time (s) bias var R M h -1 N Cost 1 5.

Nested Monte Carlo for compound option pricing

A compound option is simply an option on an option. The exercise payoff of a compound option involves the value of another option. A compound option has then two expiration dates T 1 < T 2 and two strike prices K 1 and K 2 . We consider here the example of a European style Put on a Call where the underlying risky asset S is still given by a Black-Scholes process with parameters (r, σ). On the first expiration date T 1 , the holder has the right to sell a new Call option using the strike price K 1 .

The new Call has expiration date T 2 and strike price K is given by w i = (-1) R-i n α(R-1) i

1 j<i (n α i -n α j ) i<j R
(n α jn α i )

.

(b) Furthermore

w R+1 = R i=1 w i n αR i = (-1) R-1 1 i R n α i .
Proof. (a) Let a i = n -α i . Note that by Cramer's rule the solution of this linear system is given by w i = det(V i ) det(V ) where V i is the matrix formed by replacing the ith column of V by the column vector e 1 . The first point is that V i is again a Vandermonde matrix of type V i = V (1, . . . , a i-1 , 0, a i+1 , . . . , a R ). On the other hand, the determinant of a square Vandermonde matrix can be expressed as det(V ) = 1 j<k n (a ka j ). So we have for every i ∈ 1, . . . , R

w i = 1 j<k R;j,k =i (a k -a j ) 1 j<i (-a j ) i<k R a k 1 j<k R (a k -a j ) = 1 j<i (-a j ) i<k R a k 1 j<i (a i -a j ) i<k R (a k -a i )
Using that a i = n -α i , i = 1, . . . , R, we have 1 j<i (-a j )

1 j<i (a ia j )

= n α(i-1) i 1 j<i (n α in α j ) and i<k R a k i<k R (a ka i ) = (-1) R-i n α(R-i) i i<k R (n α kn α i ) which completes the proof. (b) follows by setting x = 0 in the decomposition 1

1 i R (x -n α i ) = R i=1 1 (x -n α i ) j =i (n α i -n α j )
.

Proposition A.2. When n i = M i-1 , i = 1, . . . , R, the following holds true for the coefficients w i = w i (R, M ).

1. Closed form for w i , i = 1, . . . , R:

w i = w i (R, M ) = (-1) R-i M -α 2 (R-i)(R-i+1) 1 j i-1 (1 -M -jα ) 1 j R-i (1 -M -jα )
, i = 1, . . . , R.

2. Closed form for w R+1 :

w R+1 = (-1) R M -R(R-1) 2 α .

A useful upper bound:

sup

R∈N * R-1 i=1 | w i (R, M )| M -α π 2 α,M k 0 M -α k(k+3)
and max

1 i R | W i (R, M ) -1| W α (M ) -1 ∼ M -α → 0 as M → +∞.
In particular, the matrix T = T(R, M ) in (ML2R) converges toward the matrix of the standard Multilevel Monte Carlo (MLMC) at level M when M → +∞.

6. One more useful inequality B Appendix: sketch of proof of Propositions 2.2 and 2.4

∀ R ∈ N, 1 | w R+1 | R r=1 | w r (R, M )| n αR r B α (M ) 1 π 2 α,M k 0 M -α 2 k(k+1) .
The multistep Richardson-Romberg estimator with the formal framework of Section 3, is characterized by the allocation matrix T = (w, 0, . . . , 0). Note that the first column is not e 1 but this has no influence on what follows. The expansion of E Ȳ N h,n follows from Proposition 2.4. No stratification is needed here since only one Brownian motion is involved. The proof of Theorem 3.8 applies here with q = (1, 0, . . . , 0). Furthermore

φ( Ȳ N h,n ) =

  ) and 6(b) are given in Appendix C.

Figure 2 :

 2 Figure 2: Lookback option in a Black-Scholes model. CPU-time (y-axis, log scale) as a function of εL (x-axis, log 2 scale).

Figure 3 :

 3 Figure 3: Barrier option in a Black-Scholes model. CPU-time (y-axis, log scale) as a function of εL (x-axis, log 2 scale).

Figure 4 :

 4 Figure 4: Nested compound option in a Black-Scholes model. CPU-time (y-axis, log scale) as a function of εL (x-axis, log 2 scale).

2

  and1 w R (R, M ) 1 π α,Mwhere π α,M = k 1 (1 -M -αk ).

|

  4. Asymptotics of the coefficients wi when M → +∞: w i (R, M )| = 0 and lim M →+∞ sup R∈N * | w R (R, M ) -1| = 0.5. Asymptotics of the coefficientsW i = W i (R, M ) when M → +∞: the coefficients W i are defined in(ML2R). It follows from what precedes that they satisfy W 1 = 1, max 1 i R | W i (R, M )| W α (M ) := M

Proof. Claim 6 :

 6 For every r ∈ {1, . . . , R},| w r (R, Noting that ((Rr)(Rr + 1) + 2(r -1)R) = R(R -1) + r(r -1), we derive that R r=1 | w r (R, M )| n αR r 1 π 2 α,M M -α R(R-1) 2 R r=1 M -α r(r-1)2 which yields the announced inequality since M -α R(R-1) 2 = | w R+1 |.

Figure 5 :

 5 Figure 5: Call option in a Black-Scholes model (log 2 scale for the x-axis)

Figure 6 :Figure 7 :

 67 Figure 6: Lookback option in a Black-Scholes model (log 2 scale for the x-axis)

Figure 8 :

 8 Figure 8: Nested compound option in a Black-Scholes model (log 2 scale for the x-axis)

  00•10 -01 1.35•10 +00 1.47•10 -03 -1.32•10 +00 6.60•10 -02 2 8 1 1.17•10 +03 2.05•10 +03 2 2.50•10 -01 6.86•10 -01 1.13•10 -02 -6.72•10 -01 1.87•10 -02 3 6 1 6.80•10 +03 1.61•10 +04 3 1.25•10 -01 3.00•10 -01 6.27•10 -02 -2.91•10 -01 5.37•10 -03 4 6 1 3.59•10 +04 1.11•10 +05 4 6.25•10 -02 1.96•10 -01 2.73•10 -01 -1.92•10 -01 1.57•10 -03 4 8 1 1.49•10 +05 5.04•10 +05 5 3.12•10 -02 9.25•10 -02 1.46•10 +00 -9.03•10 -02 4.04•10

	-04 5	7	1	7.26•10 +05 2.93•10 +06
	6 1.56•10 -02 4.38•10 -02 6.80•10 +00 -4.25•10 -02 1.20•10 -04 5 10	1	3.10•10 +06 1.40•10 +07
	7 7.81•10 -03 2.47•10 -02 3.26•10 +01 -2.42•10 -02 2.87•10 -05 6	8	1	1.42•10 +07 7.17•10 +07
	8 3.91•10 -03 9.06•10 -03 1.72•10 +02 -8.64•10 -03 7.49•10 -06 7	8	1	6.62•10 +07 3.89•10 +08
	9 1.95•10 -03 6.16•10 -03 7.34•10 +02 -6.00•10 -03 1.97•10 -06 7	9	1	2.71•10 +08 1.66•10 +09

Table 5

 5

	.6: Lookback option (α = 0.5, β = 1): Parameters and results of MLMC estimator.

Table 5

 5

.7: Barrier option (α = 0.5, β = 0.5): Parameters and results of ML2R estimator.

  00•10 -01 7.83•10 -01 2.26•10 -03 7.25•10 -01 8.73•10 -02 2 8 1 1.36•10 +03 2.83•10 +03 2 2.50•10 -01 4.03•10 -01 2.05•10 -02 3.67•10 -01 2.75•10 -02 3 6 1 1.03•10 +04 3.57•10 +04 3 1.25•10 -01 1.81•10 -01 1.83•10 -01 1.56•10 -01 8.30•10 -03 4 6 1 7.18•10 +04 4.28•10 +05 4 6.25•10 -02 1.09•10 -01 9.52•10 -01 9.71•10 -02 2.47•10 -03 4 8 1 3.27•10 +05 2.40•10 +06 5 3.12•10 -02 5.33•10 -02 8.38•10 +00 4.70•10 -02 6.27•10

	-04 5	7	1	2.11•10 +06 2.40•10 +07
	6 1.56•10 -02 2.61•10 -02 6.16•10 +01 2.22•10 -02 1.88•10 -04 5 10	1	1.09•10 +07 1.74•10 +08
	7 7.81•10 -03 1.41•10 -02 4.90•10 +02 1.23•10 -02 4.51•10 -05 6	8	1	6.40•10 +07 1.43•10 +09
	8 3.91•10 -03 5.58•10 -03 6.05•10 +03 4.43•10 -03 1.15•10 -05 7	8	1	4.37•10 +08 1.67•10 +10

Table 5

 5

	.8: Barrier option (α = 0.5, β = 0.5): Parameters and results of MLMC estimator.

Note that these optimal parameters depend only on the structural parameters and on the user's choice of the root M 2 for the refiners. For a fixed ε > 0, if we emphasize the dependance in M = M (ε) i.e. R(M ), h(M ), q(M ) and N (M ) the global cost C ε as a function of M is given by

where κ(h, R, q) = 1 h R j=1 q j R i=1 n i 1

(in the framework of Section 4.1) and κ(h, R, q) = 1 h R j=1 q j max 1 i R n i 1

(in the framework of Section 4.2). This function can be optimized for likely values of M . In numerical experiments we consider

Simulating consistent Brownian increments In many situations (like e.g. the numerical experiments carried out below), discretization schemes of Brownian diffusions need to be simulated with various steps (say T nn i and T nn i+1 in our case). This requires to simulate consistent Brownian increments over [0, T n ], then [ (k-1)T n , kT n ], k = 2, . . . , n. This can be performed by simulating recursively the Brownian increments over all successive sub-intervals of interest, having in mind that the "quantum" size for the simulation is given by T nm where m = gcd(n 1 , . . . , n R ). One can also produce once and for all an abacus of coefficients to compute by induction the needed increments from small subintervals up to the root interval of length T n . This is done e.g. in [START_REF] Pagès | Multi-step Richardson-Romberg extrapolation: remarks on variance control and complexity[END_REF] up to R = 5 for α = 1 and up to R = 3 for α = 1 2 .

Methodology

We compare the two MLMC and ML2R estimators for different biased problems. In the sequel, we consider the allocation matrix (ML2R) for the ML2R estimator. After a crude evaluation of var(Y 0 ) and V 1 (using (42)) we compute the "optimal" parameter M solution to (46). The others parameters are specified according to Tables (5.1) and (5.2) with c = c 1 = 1. In the numerical simulations we do not round down R, we round to the nearest integer.

The empirical bias error of the estimator Ȳ N,q h,n is estimated by

using L = 256 independent replications, where

is the true value. To compute this bias error, it is necessary to use many replications of the estimator.

Let νN be the empirical (unitary) variance of the estimator Ȳ N,q h,n defined by

. Note that νN is an estimator of var( Ȳ N,q h,n ) that we compute online to obtain a (biased) confidence interval for one run of Ȳ N,q h,n . Since we use L independent replications of Ȳ N,q h,n (to obtain the bias error), then we consider the empirical means over the L = 256 runs prescribed error ε = 2 -k , ML2R is about 10% to 100% (twice) faster than MLMC when k goes from 2 to 8. On the other hand, both estimators ML2R and MLMC provide an empirical RMSE close to the prescribed RMSE i.e. εL ε (see also Figure 5(a) in Appendix C). We can conclude that the automatic tuning of the algorithm parameters is satisfactory for both estimators.

In Figure 1 is depicted the CPU-time (4th column) as a function of the empirical L 2 -error (3rd column). It provides a direct comparison of the performance of both estimators. Each point is labeled by the prescribed RMSE ε = 2 -k , k = 1, . . . , 8 for easy reading. The plot is in log 2 -log scale. The ML2R estimator (blue solid line) is below the MLMC estimator (red dashed line). The ratio of CPU-times for a given εL shows that ML2R goes from 1.28 up to 2.8 faster, within the range of our simulations. In Appendix C, Figure 5(b) represents the product (CPU-time)×ε 2 as a function of ε. Note that this function is almost constant (as expected) for both estimators.

Lookback option

We consider a partial Lookback Call option defined by its functional payoff

where λ

1. The parameters of the Black-Scholes model are s 0 = 100, r = 0.15, σ = 0.1 and T = 1 and the coefficient λ is set at λ = 1.1. For these parameters, the price given by a closed-form expression is I 0 = 8.89343.

For such functional Lipschitz continuous payoff, (SE β ) holds with β = 1 and (W E α, R) holds with α = 0.5. Note that the full expansion R = +∞ is not yet proved to our knowledge. An estimation of others structural parameters yields var(Y 0 ) ≃ 41 and V 1 ≃ 3.58 (and then θ ≃ 0.29). Both estimators are implemented using the automatic tuning previously exposed.

The results are summarized in Table 5.5 for the ML2R and in Table 5.6 for the MLMC. Note first that as a function of the prescribed ε = 2 -k the ratio between CPU-times goes from 1.1 (k = 2) up to 3.5 (k = 9), idem for the ratio Cost(MLMC)/ Cost(ML2R). However, the empirical RMSE of MLMC is greater than ε (certainly because c 1 = 1) unlike that of ML2R (see 2(a) in Appendix C). Table 5.5: Lookback option (α = 0.5, β = 1): Parameters and results of ML2R estimator.

5.3.3

Barrier option α = 0.5, β = 0.5

We consider now an up-and-out call option to illustrate the case β = 0.5 < 1 and α = 0.5. This path-dependent option with strike K and barrier B > K is defined by its functional payoff

The parameters of the Black-Scholes model are s 0 = 100, r = 0, σ = 0.15 and T = 1. With K = 100 and B = 120, the price obtained by closed-form solution is I 0 = 1.855225. We consider here a simple (and highly biased) approximation of max To comply with the multilevel framework, we set

Note that the underlying process (S t ) t∈[0,T 2 ] is not discretized in time. The bias error is then due to the inner Monte Carlo estimator of the conditional expectation.

The parameters used for the underlying process (S t ) t∈[0,T 2 ] are S 0 = 100, r = 0.03 and σ = 0.3. The parameters of the Put-on-Call payoff are T 1 = 1/12, T 2 = 1/2 and K 1 = 6.5, K 2 = 100. Section 4.2 strongly suggest that (SE β ) and (W E α, R) are satisfied with β = α = 1. A crude computation of others structural parameters yields var(Y 0 ) ≃ 9.09, V 1 ≃ 7.20 and θ ≃ 0.89.

The results are summarized in Table 5.9 for ML2R and in Table 5.10 for MLMC.

A Appendix

Lemma A.1. (a) The solution of the system V w = e 1 where V is a Vandermonde matrix