Pierre Henry-Labordère

Xiaolu Tan
email: tan@ceremade.dauphine.fr

Nizar Touzi
email: nizar.touzi@polytechnique.edu

Unbiased simulation of stochastic differential equations *

Keywords: Primary 65C05, 60J60; secondary 60J85, 35K10 Unbiased simulation of SDEs, regime switching diffusion, linear parabolic PDEs

We propose an unbiased Monte-Carlo estimator for E[g(X t 1 , • • • , X tn )], where X is a diffusion process defined by a multi-dimensional stochastic differential equation (SDE). The main idea is to start instead from a well-chosen simulatable SDE whose coefficients are updated at independent exponential times. Such a simulatable process can be viewed as a regime-switching SDE, or as a branching diffusion process with one single living particle at all times. In order to compensate for the change of the coefficients of the SDE, our main representation result relies on the automatic differentiation technique induced by Bismu-Elworthy-Li formula from Malliavin calculus, as exploited by Fournié et al. [14] for the simulation of the Greeks in financial applications. In particular, this algorithm can be considered as a variation of the (infinite variance) estimator obtained in Bally and Kohatsu-Higa [3, Section 6.1] as an application of the parametrix method.

MSC2010.

Introduction

Let d ≥ 1, T > 0 and W be a d-dimensional Brownian motion, µ : [0, T ] × R d → R d and σ : [0, T ] × R d → M d be the drift and diffusion coefficients, where M d denotes the collection of all d × d dimensional matrices. Under standard assumptions on these coefficients, we consider the process X defined as the unique strong solution of the multi-dimensional SDE, X 0 = x 0 , and dX t = µ t, X t dt + σ t, X t dW t , (

Our main focus in this paper is on the Monte-Carlo approximation of the expectation

V 0 := E g X t 1 , • • • , X tn , (1.2)
for some function g : R d×n → R and discrete time grid 0 < t 1 < • • • < t n = T .

When n = 1, the analytic formulation of the problem is obtained by the well-known representation V 0 = u(0, X 0 ), where u is the solution of the linear PDE

∂ t u + b(t, x) • Du + 1 2 Tr σσ ⊤ (t, x)D 2 u = 0, u T = g, (1.3)
In practice, the Monte-Carlo method consists in simulating N independent copies of a discrete-time approximation of X, and then estimating V 0 by the empirical mean value of the simulations. The corresponding error analysis consists of a statistical error induced by the central limit theorem, and a discretization error which induces a biased Monte Carlo approximation. Under some smoothness conditions, Talay and Tubaro [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF] proved that the discretization error for the Euler scheme is controlled by a rate C∆t, where ∆t denotes the time step discretization. Since then, many works focused on the analysis of the discretization error under various discretization techniques, see e.g. Kloeden and Platen [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF], and Graham and Talay [START_REF] Graham | Stochastic Simulation and Monte Carlo Methods[END_REF] for an overview. However, the statistical error estimate N -1 2 is lost in all cases, as its combination with the discretization error leads to an overall error estimate of the order N -1 2 +ε for some ε > 0.

In the context of one-dimensional homogeneous SDEs with constant volatility coefficient, Beskos and Roberts [START_REF] Beskos | Exact simulation of diffusions[END_REF] developed an exact simulation technique for X by using the Girsanov change measure together with a rejection algorithm, see also Beskos, Papaspiliopoulos and Roberts [START_REF] Beskos | Retrospective exact simulation of diffusion sample paths[END_REF], Jourdain and Sbai [START_REF] Jourdain | Exact retrospective Monte Carlo computation of arithmetic average Asian options[END_REF], etc... This technique also applies to more general SDEs by using of the so-called Lamperti transformation which reduces the approximation problem to the unit diffusion case. We also refer to the subsequent active literature of exact simulation of an L ∞ -approximation of X, see. e.g Blanchet, Chen and Dong [START_REF] Blanchet | ǫ-Strong Simulation for Multidimensional Stochastic Differential Equations via Rough Path Analysis[END_REF].

An alternative approximation method for V 0 was induced by the multilevel Monte-Carlo (MLMC) algorithm introduced by Giles [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF], which generalizes the statistical Romberg method of Kebaier [START_REF] Kebaier | Statistical Romberg extrapolation: A new variance reduction method and applications to option pricing[END_REF]. One of the main advantages of the MLMC algorithm is to control the global error (sum of discretization error and statistical error) with a much better rate w.r.t. the computation complexity. We refer to Giles and Szpruch [START_REF] Giles | Antithetic multilevel Monte Carlo estimation for multi-dimensional SDEs without Lévy area simulation[END_REF], Alaya and Kebaier [START_REF] Alaya | Central limit theorem for the multilevel Monte Carlo Euler method[END_REF], Rhee and Glynn [START_REF] Rhee | Unbiased estimation with square root convergence for SDE models[END_REF] for further developments. In particular, Rhee and Glynn [START_REF] Rhee | Unbiased estimation with square root convergence for SDE models[END_REF] obtained an unbiased simulation method for SDEs by using a random level in the estimator.

More recently, Bally and Kohatsu-Higa [START_REF] Bally | A probabilistic interpretation of the parametrix method[END_REF] provided a probabilistic interpretation of the parametrix method for PDEs. In particular, when n = 1, they obtained a representation formula for V 0 in form

E g X T W T , (1.4)
where X is defined by a Euler scheme of X on a random discrete-time grid (the time step follows an independent exponential distribution), and W T is a corrective weight function depending on X. The above representation is formally similar to the stochastic finite element method proposed by Bompis and Gobet [START_REF] Bompis | Stochastic Approximation Finite Element method for analytical approximation of multidimensional diffusion[END_REF], where one replaces X by its Euler scheme solution in (1.2) and then corrects partially the error by some well-chosen weight functions. Notice that in the above representation of Bally and Kohatsu-Higa [START_REF] Bally | A probabilistic interpretation of the parametrix method[END_REF], the process X can be exactly simulated and hence it may provide an unbiased estimator for V 0 . Nevertheless, the obtained weight function W T is integrable but has infinite variance, and hence the corresponding Monte-Carlo estimator looses the standard central limit error estimate.

In this paper, we provide a representation of V 0 in the spirit of (1.4), but with an alternative weight function for the representation. Our results follow from completely different arguments. More importantly, our unbiased approximation of V 0 has finite variance, and applies for a large class of SDEs.

Our main idea is to consider the Euler scheme solution X as solution to a regimeswitching SDE with some well-chosen coefficients. In order to compensate for the change of the coefficients of the SDE, we introduce some weight functions obtained by the automatic differentiation technique induced by Bismut-Elworthy-Li formula from Malliavin calculus, as exploited by Fournié et al. [START_REF] Fournié | Applications of Malliavin calculus to Monte Carlo methods in finance[END_REF] for the simulation of the Greeks in financial applications.

The technique introduced in the present paper is inspired by the numerical algorithm introduced in [START_REF] Henry-Labordère | Counterparty Risk Valuation: A Marked Branching Diffusion Approach[END_REF][START_REF] Henry-Labordère | A numerical algorithm for a class of BSDE via branching process[END_REF], for semilinear PDEs of the form ∂ t u + 1 2 ∆u + F 0 (t, x, u) = 0, u T = g, for some nonlinearity F 0 . The main idea in [START_REF] Henry-Labordère | Counterparty Risk Valuation: A Marked Branching Diffusion Approach[END_REF][START_REF] Henry-Labordère | A numerical algorithm for a class of BSDE via branching process[END_REF] is to use an approximation by a branching diffusion representation induced by approximating the nonlinearity F 0 by a polynomial in u. Namely, given the nature of the linear operator, the representation is obtained by means of a Brownian motion with branching driven by the polynomial approximation of F 0 . Loosely speaking, the method developed in the present paper follows by reading the PDE part of (1.3) in the following equivalent form:

∂ t u + 1 2 ∆u + F 1 (t, x, Du, D 2 u) = 0,
where

F 1 (t, x, z, γ) := b(t, x) • z + 1 2 Tr (σσ ⊤ (t, x) -I)γ .
However, in contrast with the nonlinearity F 0 , the above function F 1 involves the gradient and the Hessian of the solution u. Consequently the last PDE can not be handled by the existing literature on branching diffusion representation of PDEs. The automatic differentiation technique introduced in the present paper is an important new idea which allows to convert Du and D 2 u in F 1 into u. Since no powers of u are involved in the equation, this leads to a representation by means of a Brownian motion with exactly one descendent with two different possible types revealed by the weight function corresponding to the order of differentiation. We believe that the automatic differentiation trick introduced here has very important consequences, beyond the particular application of the present paper. Indeed, in our paper [START_REF] Henry-Labordère | Branching diffusion representation of semilinear PDEs and Monte Carlo approximation[END_REF], we provide a significant extension of the branching diffusion representation to a general class of semilinear PDEs.

The paper is organized as follows. In Section 2, we consider the SDE with constant diffusion coefficient, and propose an unbiased estimator for V 0 for both Markovian case and path-dependent case. Then in Section 3, we consider the SDE with general diffusion coefficient function, and obtain a similar representation formula for V 0 , which is integrable but of infinite variance. Next, in Section 4, we provide some numerical examples. Finally, we complete some technical proofs in Section 5. In particular, an easy example is studied in Section 5.1 to illustrate the main idea of the technical proofs.

Unbiased simulation of SDE with constant diffusion coefficient

In this section, we will restrict to the constant diffusion coefficient case, and propose an unbiased estimator for V 0 having finite variance.

The Markovian case

Let us start by the Markovian case, where the diffusion process X is defined by

X 0 = x 0 , dX t = µ(t, X t ) dt + σ 0 dW t , (2.1)
for some matrix σ 0 ∈ M d , and our objective is to compute

V 0 = E[g(X T )]. (2.2)
for some function g : R d → R. We impose the following conditions on µ and σ 0 :

Assumption 2.1. The diffusion coefficient σ 0 is non-degenerate, the drift function µ(t, x) is bounded continuous in (t, x), uniformly 1 2 -Hölder in t and uniformly Lipschitz in x, i.e. for some constant L > 0,

µ(t, x) -µ(s, y) ≤ L |t -s| + x -y , ∀(s, x), (t, y) ∈ [0, T ] × R d . (2.3)

The unbiased simulation algorithm

To introduce our unbiased simulation algorithm, let us first introduce a random discrete time grid. Let β > 0 be a fixed positive constant, (τ i ) i>0 be a sequence of i.i.d. E(β)exponential random variables. We define

T k := k i=1 τ i ∧ T, k ≥ 0, and N t := max k : T k < t . (2.4)
Then (N t ) 0≤t≤T is a Poisson process with intensity β and arrival times (T k ) k>0 . We denote also T 0 := 0 and ∆T k+1 := T k+1 -T k .

Let W be a d-dimensional Brownian motion independent of (τ i ) i>0 , we introduce

∆W T k := W T k -W T k-1 , k > 0.
and a process X as the Euler scheme of X on the random discrete grid (T k ) k≥0 , i.e. X 0 = x 0 and

X T k+1 := X T k + µ T k , X T k ∆T k+1 + σ 0 ∆W T k+1 , k = 0, 1, • • • , N T . (2.5)
Then our estimator is given by

ψ := e βT g X T -g X T N T 1 {N T >0} β -N T N T k=1 W 1 k , (2.6)
with

W 1 k := µ(T k , X T k ) -µ(T k-1 , X T k-1 ) • (σ ⊤ 0 ) -1 ∆W T k+1 ∆T k+1 . (2.7)
Theorem 2.2. Suppose that Assumption 2.1 holds true, and g is Lipschitz. Then

E ψ 2 < ∞ and V 0 = E ψ . Proof. (i) We first show that E ψ 2 < ∞. For simplicity, we denote ∆ X k := X T k -X T k-1 for k > 0.
Let L g be the Lipschitz constant of the function g, and set L 0 := σ 0 σ ⊤ 0 -1 > 0 by the non-degeneracy of σ 0 . Then using Assumption 2.1, it follows by direct computation that

e -βT ψ ≤ L g |g(x 0 )| + ∆T 1 + |∆ X T 1 | N T k=1 L( ∆T k+1 + |∆ X T k+1 |) β∆T k+1 (σ ⊤ 0 ) -1 ∆W T k+1 .
Then denoting

E T k := E • X T k , ∆T k+1 , we have E T k ∆T k+1 + |∆ X T k +1 | ∆T k+1 (σ ⊤ 0 ) -1 ∆W T k+1 2 ≤ E 1 + |µ| ∞ √ T + σ 0 Z 2 (σ ⊤ 0 ) -1 Z 2 ,
where

|µ| ∞ := d i=1 |µ i | 2 0 , |µ i | 0 := sup t,x |µ i (t, x)|, and Z is a standard centered normal distribution in R d . This provides E T k ∆T k+1 + |∆ X T k +1 | ∆T k+1 (σ ⊤ 0 ) -1 ∆W T k+1 2 ≤ 2 1 + |µ| ∞ √ T ) 2 E (σ ⊤ 0 ) -1 Z 2 + 2E σ 0 Z 2 (σ ⊤ 0 ) -1 Z 2 = 2 1 + |µ| ∞ √ T ) 2 Tr((σ 0 σ ⊤ 0 ) -1 ) + 2 3d + d(d -1) =: γ.
We therefore get the following upper bound:

E ψ 2 ≤ Ce 2βT e -βT + γL 2 T β , where C := L 2 g E |g(0)| + ∆T 1 + |∆X 1 | 2 .
(2.8)

(ii) The equality V 0 = E[ ψ] will be proved in Section 5, with illustration of the main idea in Section 5.1.

On the choice of β

Notice that the random variable ψ in (2.6) can be exactly simulated from a sequence of Gaussian N (0, 1) and exponential E(β) random variables. Then the integrability and representation results in Theorem 2.2 induce an unbiased simulation Monte-Carlo method to approximate V 0 , with error induced by the standard central limit theorem. We next observe that the constant β > 0 may be chosen so as to minimize the approximation error relative to the computational effort:

• By the central limit theorem, the error induced by the Monte Carlo estimator based on the representation ψ is characterized by the variance of ψ. For tractability reasons, we shall instead replace it by the bound (2.8).

• The computation effort is proportional to the number N T of arrivals of the Poisson process before the maturity T , and is thus given by

C ′ E[N T ] = C ′ βT .
In view of this, we shall choose β by minimizing the ratio of the variance bound (2.8) to the mean computational effort. This minimization problem is obviously independent of the constants C, C ′ , and reduces to:

min β>0 f (β), where f (β) := 1 βT exp T β + γL 2 β .
Direct computation shows that the equation f ′ (β) = 0 has a unique solution on (0, ∞) given by

β * := γL 2 + T 2 /4 + T 2 .
As lim βց0 f (β) = lim β→∞ f (β) = ∞, this shows that β * is the minimizer of the above defined criterion, and will be taken as our "best sub-optimal" choice of β for the unbiased estimator ψ.

The path-dependent case

In this part, we would like to provide an extension of the above estimator ψ in (2.6) to the path-dependent case. Let n > 0, 0 = t 0 < t 1 < • • • < t n = T , σ 0 ∈ M d be a non-degenerate matrix, and µ : [0, T ] × R d×n → R d be a continuous function, Lipschitz in the space variable. Let X be the unique solution of SDE, with initial condition

X 0 = x 0 , dX t = µ(t, X t 1 ∧t , • • • , X tn∧t ) dt + σ 0 dW t ; (2.9)
and the objective is to compute the value,

V 0 := E g X t 1 , • • • , X tn , (2.10)
for some Lipschitz function g : R d×n → R.

Remark 2.3. It is clear that the value V 0 defined above can be characterized by a parabolic PDE system. Namely, for every k

= 1, • • • , n and (x 1 , • • • , x k-1 ) ∈ R d×(k-1) , we define µ k (t, x) := µ(t, x 1 , • • • , x k-1 , x, • • • , x), ∀(t, x) ∈ [t k-1 , t k ] × R d . (2.11) Suppose that (u k ) k=1,••• ,n is a family of functions such that u k is defined on [t k-1 , t k ] × R d×k and x → u k (t, x 1 , • • • , x k-1 , x
) is a solution (at least in the viscosity sense) of

∂ t u k + 1 2 σ 0 σ ⊤ 0 : D 2 u k + µ k • Du k = 0, (2.12)
with terminal conditions

u k (t k , x 1 , • • • , x k ) = u k+1 (t k , x 1 , • • • , x k , x k ), for k = 1, • • • , n -1,
and u n (t n , x 1 , • • • , x n ) = g(x 1 , • • • , x n ).
Then we have V 0 = u 1 (0, x 0 ).

The algorithm

The unbiased simulation algorithm of V 0 can be obtained by an iteration of the estimator (2.6) on each time interval [t k , t k+1 ]. One should just be careful on the integrability issue. Let us first introduce the algorithm.

Recall that W be a standard d-dimensional Brownian motion, (τ i ) i>0 is a sequence of i.i.d. E(β)-exponential random variables independent of W . Then N = (N s ) 0≤s≤t and (T i ) i>0 are defined in (2.4). Define further for every k

= 1, • • • , n, Ñ k := N t k - N t k-1 the number of jump arrivals on [t k-1 , t k ), and T k 0 := t k-1 and T k j := T Nt k-1 +j ∧t k , ∆ T k j := T k j -T k j-1 , W k j := W T k j , ∆ W k j := W k j -W k j-1 , ∀j = 1, • • • , Ñ k + 1.
Example 2.4. We give below an example for the case n = 2. In the following example, the number of jump arrivals on [0,

t 1 ) is Ñ 1 = 2, that on [t 1 , t 2 ) is Ñ 2 = 1
, and total number of jump arrivals is N T = 3.

For k = 1, we have

T 1 0 = 0, T 1 1 = T 1 , T 1 2 = T 2 and T 1 3 = t 1 ; W 1 0 = 0, W 1 1 = W T 1 , W 1 2 = W T 2 and W 1 3 = W t 1 . For k = 2, we have T 2 0 = t 1 , T 2 1 = T 3 , T 2 2 = t 2 , and W 2 0 = W t 1 , W 2 1 = W T 3 and W 2 2 = W t 2 . ✲ 0 T 1 T 2 t 1 T 3 t 2
We next introduce a process X k,x j

, ∀j = 0, 1, • • • , N k + 1, for each k = 1, • • • , n and initial condition x = (x 0 , x 1 , • • • , x k-1 ) ∈ R d×k by X k,x 0 := x k-1 and X k,x j+1 := X k,x j + µ k T k j , X k,x j ∆ T k j+1 + σ 0 ∆ W k j+1 .
Similarly, for every j = 1, • • • , N k , we define a automatic differentiation weight, with µ k defined by (2.11),

W k j := µ k T k j , X k,x j -µ k T k j-1 , X k,x j-1 • σ ⊤ 0 -1 ∆ W k j+1 ∆ T k j+1 .
7

We now introduce the algorithm for the path-dependent case, in a recursive way. First, for x = (x 0 , x 1 ,

• • • , x n ) ∈ R d×(n+1) , set ψ x n+1 := g(x 1 , • • • , x n ). Next, for k = 1, • • • , n, denote X k,x := (x 0 , x 1 , • • • , x k-1 , X k,x Ñ k +1 ) and X k,x,0 := (x 0 , x 1 , • • • , x k-1 , X k,x Ñ k 1 { Ñ k >0} ).
Then given ψ • k+1 , we define

ψ x k := e β(t k -t k-1 ) ψ X k,x k+1 -ψ X k,x,0 k+1 1 { Ñ k >0} β -Ñ k Ñ k j=1 W k j .
(2.13)

We finally obtain the numerical algorithm of the path-dependent case:

ψ := ψ x 0 1 . (2.14)

The integrability and representation result

We notice that the algorithm in the path-dependent case is nothing else than an iterative algorithm of the Markovian case, as suggested by the PDEs (2.12) in Remark (2.3). When the random variable ψ in (2.14) is integrable, it is not surprising to obtain the representation V 0 = E ψ as a consequence of Theorem 2.2. However, because of the renormalization term (i.e.

ψ X k,x k+1 -ψ X k,x,0 k+1 1 { Ñ k >0} in (2.13
)), the variance analysis becomes less obvious. We provide here a sufficient condition to ensure that ψ has finite variance.

Theorem 2.5. Suppose that µ : [0, T ] × R d×n → R d and g : R d×n → R are differentiable up to the order n, with bounded derivatives. Then

E ψ 2 < ∞ and V 0 := E ψ .
We will prove the integrability result here, and leave the proof of the representation result V 0 := E ψ in Section 5. As preparation, let us first provide two technical lemmas

. Let π = (0 = s 0 < s 1 < • • • < s m = T ) be an arbitrary partition of the interval [0, T ], μ : [0, T ] × R d → R d a R d -valued function. We define X π,x by X π,x 0 := x and X π,x k+1 := X π,x k + μ s k , X π,x k ∆s k+1 + W s k+1 -W s k . (2.15) Further, let ϕ : R d → R be a smooth function, ℓ > 0 and i = (i 1 , • • • , i ℓ ) ∈ {1, • • • , d} ℓ , we denote ∂ ℓ x,i ϕ(x) := ∂ ℓ x i 1 •••x i ℓ ϕ(x).
Lemma 2.6. Suppose that x → μ(t, x) is differentiable up to order n with uniformly bounded derivatives, and X π,x is defined by (2.15) with initial condition X π,x 0 = x. Then x → X π,x k is differentiable up to order n and there is a constant C independent of the partition π such that

max 1≤ℓ≤n max i∈{1,••• ,d} ℓ max 0≤k≤m ∂ ℓ x,i X π,x k ≤ C.
Proof. For simplicity, we consider the one dimensional d = 1 case, while the multidimensional can be deduced by almost the same arguments. First, let ℓ = 1, we have

∂ x X π,x k+1 = ∂ x X π,x k + ∂ x μ s k , X π,x k ∂ x X π,x k ∆s k+1 , which implies that ∂ x X π,x k+1 = Π k+1 j=1 1 + ∂ x μ s k , X π,x k ∆s k+1 . Since ∂ x μ(t, x) is uniformly bounded, it follows that ∂ x X π,x k is bounded by some con- stant C 1 independent of 1 ≤ k ≤ m and the partition π. By induction, it is easy to deduce that for ℓ = 2, • • • , n, ∂ ℓ x ℓ X π,x k+1 = ∂ ℓ x ℓ X π,x k + P ℓ ∂ i x i μ(s k , X π,x k ), ∂ i x i X π,x k , i = 1, • • • , ℓ -1 ∆s k+1 ,
where

P ℓ is a Polynomial on ∂ i x i μ(s k , X π,x k ) and ∂ i x i X π,x k for i = 1, • • • , ℓ -1, which is uniformly bounded by some constant independent of k = 1, • • • , m and the partition π. Hence ∂ ℓ x ℓ X π,x k is also bounded by some constant C ℓ independent of k = 1, • • • , m and the partition π. Lemma 2.7. Let ( ψ x k ) 1≤k≤n+1 be defined by (2.13). Then for every k = 2, • • • , n + 1, and every x = (x 0 , x 1 , • • • , x k-1 ) ∈ R d×k , the map x k-1 → ψ x k has derivatives up to order k -1 and max 1≤ℓ≤k-1 ∂ ℓ x k-1 ,i ψ x k ≤ C n j=k ( Ñ j + 1) j-1 .
(2.16)

Proof. We will prove it by induction. First, let k = n+1, then ψ x n+1 := g(x, x 1 , • • • , x n ) and hence |∂ ℓ xn ψ x | ≤ C for some constant C and for every ℓ = 1, • • • , n. Next, suppose that (2.16) holds true for ψ x k+1 , we know from (2.13) that

ψ x k := ψ X k,x k+1 -ψ X k,x,0 k+1 1 { Ñ k >0} Ñ k j=1 µ k ( T k j , X k,x j ) -µ k ( T k j-1 , X k,x j-1 ) β∆ T k j+1 • (σ ⊤ 0 ) -1 ∆ W k j+1 .
Then using the estimation in Lemma 2.6, we see that (2.16) is also true for ψ x k , and we hence conclude the proof.

Proof of Theorem 2.5 (i). By Lemma 2.7, we know that x → ψ x,x 2 is differentiable and in particular uniformly Lipschitz with coefficient bounded by 2CΠ n j=2 ( Ñ j + 1) j-1 . Then the definition of ψ x 0 1 falls into the Markovian case n = 1, but with terminal condition x → ψ x,x 2 . Notice that Ñ k ≤ N T has a Poisson distribution:

P(N T = m) = e -βT (βT ) m m! . It follows that, for some constant C > 0, E ψ x 0 1 2 ≤ E C Ñ k 4C 2 n j=2 ( Ñ j + 1) 2(j-1) ≤ E 4C 2 C N T (N T + 1) n(n-1) < ∞,
which implies that ψ has finite variance.

Unbiased simulation of general SDEs

Let us now consider the SDE (1.1) with general diffusion coefficient function, i.e. with drift and diffusion coefficients µ :

[0, T ] × R d → R d and σ : [0, T ] × R d → M d :
X 0 = x 0 , and dX t = µ t, X t dt + σ t, X t dW t .

Our objective of study in this section is

V 0 = E g(X T
) , for some function g : R d → R.

We will provide a representation result of V 0 in the same spirit of that in Section 2.

Remark 3.1 (Lamperti's transformation). We also notice that in some cases, the above SDE (1.1) may be reduced to the constant diffusion coefficient case (2.1), by the so-called the Lamperti transformation.

(i) When d = 1 and σ(t, x) > 0, let us define a function h : [0, T ] × R → R by h(t, x) := x 0 1 σ(t, y) dy.
Notice that for fixed t ∈ [0, T ], x → h(t, x) is strictly increasing, we denote h -1 (t, •) its inverse function. Then by Itô's formula, it is easy to obtain that Y t := h(t, X t ) satisfies the SDE

dY t = ∂ t h t, h -1 (t, Y t ) + µ(t, h -1 (t, Y t )) σ(t, h -1 (t, Y t )) - 1 2 ∂ x σ t, h -1 (t, Y t ) dt + dW t ,
whose diffusion coefficient is a constant as in SDE (2.1).

(ii)When d > 1, σ is non-degenerate and satisfies some further compatibility conditions, one can also obtain a similar transformation to reduce SDE (5.8) to the constant diffusion coefficient case.

An estimator of infinite variance for general SDEs

Let us impose the following conditions on coefficient functions µ and σ.

Assumption 3.2. The function (µ, σ) : [0, T ] × R d → R d × M d and a := 1 2 σσ ⊤ : [0, T ]×R d → M d are
uniformly bounded, and are uniformly Hölder in the time variable, uniformly Lipschitz in the space variable, i.e. for some constant L,

µ, σ, a (t, x) -µ, σ, a (s, y) ≤ L |t -s| + x -y , (3.1)
for all (t, x), (s, y) ∈ [0, T ] × R d ; and σ(t, x) is non-degenerate such that, for some constant ε 0 > 0,

a(t, x) := 1 2 σσ ⊤ (t, x) ≥ ε 0 I d , ∀(t, x) ∈ [0, T ] × R d .
Recall that (T k ) k≥0 are defined by (2.4) with a sequence of i.i.d. E(β)-exponential random variables, and W is a Brownian motion; the increment of the Brownian motion are defined by ∆W t k := W t k -W t k-1 , and ∆T k := T k -T k-1 . As in (2.5), we introduce X as solution of the Euler scheme on discrete grid by X 0 := x 0 and

X T k+1 := X T k + µ T k , X T k ∆T k+1 + σ T k , X T k ∆W T k+1 k = 0, • • • , N T . (3.2)
We then introduce a representation formula by

ψ := e βT g X T -g X T N T 1 {N T >0} β -N T N T k=1 W 1 k + W 2 k , (3.3)
where, for each

k = 1, • • • , N T , W 1 k := µ(T k , X T k ) -µ(T k-1 , X T k-1 ) • σ ⊤ (T k , X T k ) -1 ∆W T k+1 ∆T k+1 ,
and

W 2 k := a T k , X T k -a T k-1 , X T k-1 : σ ⊤ (T k , X T k ) -1 ∆W T k+1 ∆W ⊤ T k+1 -∆T k+1 I d ∆T 2 k+1 σ(T k , X T k ) -1 . (3.4)
Theorem 3.3. Suppose that Assumption 3.2 holds true, and g is Lipschitz. Then

E ψ < ∞ and V 0 = E[ ψ]. Proof. (i) Consider the random vectors ξ 1 k := ∆W T k √ ∆T k and ξ 2 k := ∆W T k ∆W ⊤ T k -∆T k I d ∆T k , for all k = 1, • • • , N T + 1, which are independent of ∆T k conditional on {∆T k > 0} = {N T ≥ k -1},
and which have finite second order moment. Notice that µ(t, x) and a(t, x) are uniformly bounded, and 1/2-Hölder-continuous in t and Lipschitz in x, and σ is uniformly bounded from below above zero. Then, for each k

= 1, • • • , N T , W 1 k ≤ C ∆T k + X T k -X T k-1 ∆W T k+1 ∆T k+1 ≤ C 1 + ∆T k + |ξ 1 k | |ξ 1 k+1 | ∆T k ∆T k+1 ,
where the constant C > 0 may vary from term by term but is uniformly bounded for all k. Similarly, one obtains that

W 2 k ≤ C 1 + ∆T k + |ξ 1 k | |ξ 2 k+1 | ∆T k ∆T k+1 1 ∆T k+1 .
As ∆T k ≤ T , it follows that

W 1 k + W 2 k ≤ C 1 + |ξ 1 k | |ξ 1 k+1 | + |ξ 2 k+1 | ∆T k ∆T k+1 1 ∆T k+1 ,
for some constant C > 0 independent of k. In addition, we have by the Lipschitz condition on g that

E g( X T ) -g( X T N T ) ∆T N T +1 < C ∆T N T +1 .
Then, it follows from the expression of ψ in (3.3) that

E ψ ≤ CE N T k=1 C ∆T k+1 1 {N T ≥1} + CE g( X T ) 1 {N T =0} ≤ CE N T k=1 C ∆T k+1 + CE g(x 0 + µ(0, x 0 )T + σ(0, x 0 )W T )
for some constant C > 0, where we have also used the independence of the ξ i k 's and their the boundedness of their second order moments. The integrability of ψ is now a direct consequence of Lemma A.2.

(ii) The proof of the equality V 0 = E[ ψ] will be completed in Section 5.

To conclude, we notice that the variable ψ is of order Π N T k=1 1/ ∆T k+1 in general cases, and the latter is integrable but of infinite variance. Therefore, ψ is not a good estimator for Monte-Carlo method. Nevertheless, it should still have some theoretical value as an alternative representation formula obtained by Bally and Kohatsu-Higa [3, Section 6.1].

An estimator for one-dimensional driftless SDE

To overcome the problem of variance explosion of the estimator (3.3), we will consider the higher order approximation X of X, and obtain an estimator of finite variance for the one dimensional (d = 1) driftless SDE of form

X 0 = x 0 , dX t = σ(t, X t ) dW t , (3.5)
Our objective is to compute

V 0 := E g(X T
) , for some function g : R → R.

Recall (T k ) k≥0 has been introduced in (2.4) from a sequence of i.i.d. exponential random variables, independent of the Brownian motion W . We next define X by

X 0 = x 0 , d X t = σ(T k , X T k ) + ∂ x σ(T k , X T k ) X t -X T k dW t , on [T k , T k+1 ], (3.6
)

for k = 0, 1, • • • , N T . By denoting c k 1 := σ(T k , X T k ) -∂ x σ(T k , X T k ) X T k and c k 2 := ∂ x σ(T k , X T k ), (3.7)
then the above linear SDE (3.6) has an explicit solution which is given by

X T k+1 = X T k + σ(T k , X T k )∆W T k+1 , if c k 2 = 0, (3.8)
and

X T k+1 = - c k 1 c k 2 + c k 1 c k 2 exp - (c k 2 ) 2 2 ∆T k+1 + c k 2 ∆W T k+1 + X T k exp - (c k 2 ) 2 2 ∆T k+1 + c k 2 ∆W T k+1 , if c k 2 = 0. (3.9)
We then define ψ by

ψ := e βT g( X T ) -g( X T N T )1 {N T >0} β -N T N T k=1 W 2 k , (3.10)
where the automatic differentiation weight is given by (see Lemma 5.9 below)

W 2 k := a(T k , X T k ) -ãk 2a(T k , X T k ) -∂ x σ(T k , X T k ) ∆W T k+1 ∆T k+1 + ∆W 2 T k+1 -∆T k+1 ∆T 2 k+1 , (3.11)
with a(•)

:= 1 2 σ 2 (•), ãk := 1 2 σ2 k and σk := σ(T k-1 , X T k-1 ) + ∂ x σ(T k-1 , X T k-1 )( X T k - X T k-1 ).
Similarly to the discussion at the end of Section 3.1 (see also Remark 5.8 below), the variable ψ in (3.10) is integrable but of infinite variance in general. To make the variance finite, we introduce an alternative estimator using an antithetic variable. Let X - T be an antithetic variable of X T defined by

X - T := X T N T -σ(T N T , X T N T )∆W T N T , if c N T 2 = 0,
and

X - T = - c N T 1 c N T 2 + c N T 1 c N T 2 exp - (c N T 2 ) 2 2 ∆T N T +1 -c N T 2 ∆W T N T +1 + X T N T exp - (c N T 2 ) 2 2 ∆T N T +1 -c N T 2 ∆W T N T +1 , if c N T 2 = 0. Denote W - k := W 2 k for k = 1, • • • , N T -1 and W - N T := a(T N T , X N T ) -ãN T 2a(T N T , X N T ) ∂ x σ(T N T , X N T ) ∆W T N T +1 ∆T N T +1 + ∆W 2 T N T +1 -∆T N T +1 ∆T 2 N T +1
.

We then introduce

ψ := ψ + ψ - 2 with ψ -:= e βT g( X - T ) -g( X T N T )1 {N T >0} β -N T N T k=1 W - k . (3.12)
Notice that the Brownian motion is symmetric, thus ψ -has exactly the same distribution as ψ, and it serves as an antithetic variable.

Assumption 3.4. The diffusion coefficient σ(•) satisfies σ(t, x) ≥ ε > 0 for all (t, x) ∈ [0, T ] × R, σ(t,
(•) ∈ C 2 b (R).
Theorem 3.5. Suppose that Assumption 3.4 holds true. Then

E ψ + E ψ 2 < ∞; and V 0 = E ψ = E ψ . (3.13)
We will complete the proof in Section 5.4.

Remark 3.6. (i) As ψ has finite variance, we may use the representation of Theorem 3.5 to built an unbiased Monte-Carlo estimator of V 0 . However, given the assumed regularity conditions, and the restriction to the one-dimensional setting, such a Monte-Carlo approximation is not competitive with the corresponding PDE based approximation methods. However, we believe that the present methodology is open to potential improvements, and we hope to improve our results in some future work so as to address the higher dimensions.

(ii) For a general SDE with drift function and/or d ≥ 1, we can also consider a similar choice of (μ, σ), which leads to μ(t, x) = c 1 + c 2 x and σ(t, x) = c 3 + c 4 x and a linear SDE

d X t = c 1 + c 2 X t dt + c 3 + c 4 X t dW t , (3.14)
where

c 1 ∈ R d , c 2 , c 3 ∈ M d and

Numerical examples

Notice that our estimator ψ given by (2.6) (resp. ψ given by (2.13) and (2.14)) is an unbiased estimator for V 0 in (2.2) (resp. V 0 in (2.10)). Then the error analysis of the Monte-Carlo approximation reduces to the statistical error. Hence the computation cost to achieve the accuracy O(ε) for the approximation of V 0 (resp. V 0 ) is of order O(ε -2 ), thus avoiding of the dependence on the discretization error. By combining different level of simulations, the MultiLevel Monte Carlo (MLMC) method proposed by Giles [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF] achieves a computation cost of order O(ε -2 (log ε) 2 ) or O(ε -2 ) depending on the strong discretization error rate. In particular, by considering a randomization of the level, Rhee and Glynn [START_REF] Rhee | Unbiased estimation with square root convergence for SDE models[END_REF] obtained an unbiased estimator. In the following, we provide some numerical results and comparisons between our unbiased simulation method with the Euler based MLMC method proposed by [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF].

Two one-dimensional SDEs

Let W be a one-dimensional standard Brownian motion, we consider the SDE given by

S 0 = 1, dS t = 0.1 M ∧ S t -1 S t dt + 1 2 S t dW t ,
where M is a large constant introduced in order to guarantee the Lipschitz property of the drift coefficient (in our numerical implementation, we have observed that the value of M is not relevant for large M , and that the numerical finding are not changed by taking M = ∞; this hints that our results may be extended beyond the case of Lipschitz coefficients). Applying Lemperti's transformation X t := log(S t ), we reduce the above SDE to the constant diffusion coefficient case, in form of (2.1), We implement our unbiased simulation method for the two following expectations:

X 0 = 0, dX t = 0.1 √ M ∧ e Xt -
V 0 := E (S T -K) + , and
V 0 := E 1 n n k=1 S t k -K + , (4.2)
where we choose K = 1, T = 1, n = 10 and t k := k n T . Notice that the pathdependent example does not satisfy the differentiability sufficient condition in Theorem 2.5. However, our numerical findings do not show any numerical difficulty in the present setting.

Using different numbers N of simulations, we obtain the standard deviation as (statistical) error of our estimator. Next, using the errors obtained by our unbiased simulation method, we implement the MLMC algorithm in Section 5 of Giles [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF], and we compare the computation time (in second) of the two methods. More precisely, the statistical error of the unbiased simulation method is given by Var[ ψ]/N , where Var[ ψ] denotes the estimated variance of ψ. For the implementation of MLMC, we choose M = 4, N L = 10 4 and use equation [START_REF] Dupire | Functional Itô calculus[END_REF] in [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF] as criteria to stop the loop in MLMC (see more details in Section 5 of [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF] for the meaning of M and N L ).

The numerical results are given in Tables 1 and2. We observe that with the same Monte-Carlo error, both methods have very close performance. In the present particular example, the computational time of our methods is slightly smaller. However, the conclusion may change depending on the nature of the example. Let us consider the problem

V 0 := E sin(X T ) , (4.3)
where X is defined by SDE, for some constant µ 0 ∈ R,

X 0 = 0, dX t = µ 0 cos(X t )dt + 1 2 dW t .
We implement the MLMC algorithm and our unbiased simulation method with different value of β, but with a given fixed error ε = 0.0002. The two methods provide very close estimation of value V 0 , so we give a comparison on the computation time in Figure 1. We can observe that β in the unbiased simulation method should not be too big nor too small, to minimize the computation effort. When µ 0 = 0.2, the computation time of MLMC method is slightly longer than the US method with β ≈ 0.05. However, when µ 0 = 0.5, the computation time MLMC method is always smaller than the US method for any choice of β > 0. This shows that, in the context of the present example, the performance of our unbiased simulation method is of the order of that of the multilevel Monte Carlo method.

A multi-dimensional SDE

We next consider a d-dimensional SDE with d = 4. Let W = (W 1 , • • • , W 4 ) ⊤ be a 4-dimensional standard Brownian motion, and σ 0 the 4 × 4 be the lower triangular matrix such that

σ 0 σ ⊤ 0 =      1 1/2 1/2 1/2 1/2 1 1/2 1/2 1/2 1/2 1 1/2 1/2 1/2 1/2 1      .
We consider the SDE

dX t = µ(t, X t )dt + σ 0 dW t , X i 0 = 0, i = 1, • • • , 4, with drift function µ(t, x) = (µ i (t, x), i = 1, • • • , 4) be given by µ i (t, x 1 , • • • , x 4 ) = 0.1 3 4 exp(x i ) + 1 4 exp(x) -1 -1 8
, where exp(x) := (e x 1 + • • • + e x 4 )/4. We then consider two problems:

V 0 := E 1 4 4 i=1 e M ∧X i T -K + , and
V 0 := E 1 4n n k=1 4 i=1 e M ∧X i t k -K + , (4.4)
where we choose K = 1, T = 1, n = 10 and t k := k n T and M is a large number so as to ensure that the terminal condition is Lipschitz. As in the one-dimensional case, we implement our unbiased simulation method using different sample sizes N . Then, we use the errors, obtained from our unbiased simulation method, in the MLMC algorithm in Section 5 of Giles [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF], and we compare the computation time (in second) of the two methods.

The numerical results are given in Tables 3 and4. We observe that both methods have very similar performance, with a slightly small advantage for our method. However, similar to the one-dimensional case, the MLMC algorithm could be better in other examples.

A one-dimensional driftless SDE

Finally, we provide an example of a one-dimensional driftless SDE. We recall that under the assumed regularity in Theorem 3.5, we are not expecting our method to be competitive with the PDE based approximations. Instead, our objective is to study numerically the performance of the estimator (3.12).

Let us consider the SDE

X 0 = 1, dX t = 2σ 1 + X 2 t dW t , (4.5)
and we aim to compute

V 0 = E X T -K + . (4.6)
We implemented the simple Euler scheme with time step ∆t = 1/10 and simulation number N = 10 6 , and next the unbiased simulation method (3.12) with β = 0. 5: Numerical results for V 0 in (4.6) (case d = 1), US denotes the unbiased simulation algorithm (3.12) with β = 0.1, the computation times are expressed in second. Notice also that the unbiased algorithm (3.12) contains implicitly an antithetic variance reduction, which makes its statistical even error smaller than that of the Eurler scheme.

simulation number N = 10 6 , and then the MLMC scheme using the statistical error obtained from the unbiased simulation method. The results are given in Table 5, and we can observe that all three methods provides very similar estimation of V 0 . In particular, the unbiased simulation method has a significant advantage.

Proofs

A toy example

Before completing the technical part of the proofs for Theorems 2.2, 2.5 and 3.3, we would like to illustrate the main idea by studying a simplified example in the one dimensional case with unit diffusion:

X 0 = x 0 , dX t = µ(t, X t )dt + dW t .
Let b ∈ R and β > 0, we define a sequence of i.i.d. random variable (τ k ) k≥1 of distribution E(β). Then let (T k ) k≥1 , ∆T k+1 := T k+1 -T k and ∆W k+1 := W T k+1 -W T k be defined in and below (2.4), we introduce X by . Then for all constants β > 0 and β > 0, one has

X t := x 0 + bt + W t . and then define ψ = e βT g X T N T k=1 µ(T k , X T k ) -b ∆W T k+1 β∆T k+1 . ( 5
E ψ < ∞ and E g(X T ) = E ψ .
Proof. Since µ and g are uniformly bounded, and for some constant C > 0, the conditional expectation E ∆T k+1 [|∆W T k+1 |] ≤ C/ ∆T k+1 , then by Lemma A.2, it is obvious that ψ is integral. Then it is enough to prove that E g(X T ) = E ψ . In preparation, let us introduce

ψ n = e βT n+1 N T ∧n k=1 (µ(T k , X T k ) -b)∆W T k+1 β∆T k+1 g X T 1 N T ≤n + µ -b β ∂ x u T n+1 , X T n+1 1 N T >n .
for all n ≥ 0, with the convention 0 k=1 ≡ 1. It is clear that (ψ n ) n≥0 are all integrable by Lemma A.2.

(i) Notice that µ and g are both smooth functions, then by Feynmann-Kac formula, we know

E[g(X T )] = u(0, x 0 ), where u ∈ C ∞ b ([0, T ] × R) is a smooth function of PDE ∂ t u(t, x) + 1 2 ∂ 2 xx u(t, x) + µ(t, x)∂ x u(t, x) = 0, for all (t, x) ∈ [0, T ) × R,
with terminal condition u(T, x) = g(x). Rewriting the above PDE in the following equivalent way:

∂ t u(t, x) + b∂ x u(t, x) + 1 2 ∂ 2 xx u(t, x) + µ(t, x) -b ∂ x u(t, x) = 0,
it follows from the Feynmann-Kac formula that

u(0, x 0 ) = E g( X T ) + T 0 µ t, X t -b ∂ x u t, X t dt = E e βT g X T 1 {T 1 ≥T } + e βT 1 β µ T 1 , X T 1 -b ∂ x u T 1 , X T 1 1 {T 1 <T } (5.2) = E ψ 0 ,
where the second equality follows from the fact that T 1 = T ∧ τ 1 , and τ 1 is a random variable independent of X, with density function βe -βt 1 {t≥0} .

(ii) Next, notice that for any t and bounded continuous function φ 0 , one has by integration by parts that

∂ x E φ 0 (x + bt + W t ) = E φ 0 (x + bt + W t ) W t t . (5.3) Notice also that ∆W T 1 = W T 1 , ∆T 1 = T 1 and X T 1 := x 0 + bT 1 + W T 1 . It follows by Lemma A.3 that ∂ x u(0, x 0 ) = E e β∆T 1 ∆W T 1 ∆T 1 g X T 1 {T ≤T 1 } + µ -b β ∂ x u T 1 , X T 1 1 {T 1 <T } . (5.4)
Changing the initial condition (0, x 0 ) to (T 1 , X T 1 ), one obtains that, whenever T 1 < T ,

∂ x u(T 1 , X T 1 ) = E e β∆T 2 ∆W T 2 ∆T 2 g X T 1 {T ≤T 2 } + µ -b β ∂ x u T 2 , X T 2 1 {T 2 <T } T 1 , X T 1 .
Plugging the above expression of ∂ x u(T 1 , X T 1 ) into the r.h.s. of (5.2), and using the fact that T ≤ T 2 is equivalent to N T -≤ 1, and

P {N T -≤ 1} \ {N T ≤ 1} = 0, it follows that u(0, x 0 ) = E[ψ 1 ].
(iii) Next, changing the initial condition in (5.4) from (0, x 0 ) to (T 2 , X T 2 ) when T 2 < T , and then plugging the corresponding expression of

∂ x u(T 2 , X T 2 ) into ψ 1 , it follows that u(0, x 0 ) = E[ψ 2 ].
Repeating the procedure, we have for all n ≥ 0,

E g(X T ) = u(0, x 0 ) = E[ψ n ].
Finally sending n → ∞, and using Lemma A.2 together with the dominated convergence theorem, it follows that E g(X

T ) = E[lim n→∞ ψ n ] = E[ψ].
Remark 5.2. We can also interpret formally the representation ψ in (5.1) as the expansion of the diffusion process X around a Brownian motion. Let b = 0 and µ(t, x) ≡ µ 0 for some constant µ 0 ∈ R, so that X t = W t and X t = µ 0 t + W t . Using the fact that P(N

T = k) = e -βT (βT ) k k! , ∀k ≥ 0, it follows formally that E g(X T ) = E ∞ k=0 (µ 0 T ) k k! g (k) (W T ) = E e βT g(W T ) Π N T k=1 µ 0 ∆W k+1 β∆T k+1 , (5.5)
where the second equality follows by the fact that

P[N T = k] = e βT (βT ) k k! , and E g (k) k i=0 ∆W i+1 = E g k i=0 ∆W i+1 Π k i=1 ∆W i+1 ∆T i+1 .
In particular, the r.h.s. of (5.5) is exactly E ψ defined by (5.1) in this case.

To conclude this part, we notice that ψ in (5.1) is integrable, but has an infinite variance in general. In the next subsection, we exploit the arbitrariness of the constant b which is involved in the definition of ψ. More precisely, we shall choose different constants b at each time T k , in an adaptive way. This will lead to the estimator ψ in Theorems 2.2 and 3.3.

A regime switching diffusion representation

For d ≥ 1, T > 0, let (µ, σ) : [0, T ] × R d → R d × M d be bounded continuous functions satisfying µ(t, x) -µ(t, y) + σ(t, x) -σ(t, y) ≤ L|x -y|; (t, x, y) ∈ [0, T ] × R d × R d , (5.6)
for some constant L > 0. We start by considering a linear parabolic PDE ) s∈[0,T ] defined as unique strong solution of the SDE X 0 = x 0 , and dX s = µ s, X s ds + σ s, X s dW s , s ∈ [0, T ].

∂ t u + µ • Du + a : D 2 u = 0, on [0, T ) × R d , ( 5
(5.8)

When PDE (5.7) admits a classical solution in C 1,3 b ([0, T ]×R d ), i.e. the collection of all functions φ(t, x) such that φ, ∂ t φ, Dφ, D 2 φ and D 3 φ all exit and bounded continuous, it follows by Feynmann-Kac formula that V 0 := E[g(X 0,x 0 T )] = u(0, x 0 ). Remark 5.3. For technical reason, we will assume that u ∈

C 1,3 b ([0, T ] × R d ) rather than in C 1,2 b ([0, T ] × R d ).
But by approximating the coefficient µ, σ with smooth functions, one can relax this regularity condition in more concrete context.

Recall that for β > 0, (τ i ) i>0 is a sequence of i.i.d. E(β)-exponential random variables, which is independent of the Brownian motion W . We define

T k := k i=1 τ i ∧ T, k ≥ 0, and N t := max k : T k < t .
Then (N t ) 0≤t≤T is a Poisson process with intensity β and arrival times (T k ) k>0 , and T 0 = 0. We also introduce, for all k > 0, ∆W k

t := W (T k-1 +t)∧T k -W T k-1 . It is clear that the sequence of processes (∆W k • ) k>0 are mutually independent. Let (μ, σ) : (s, y, t, x) ∈ [0, T ] × R d × [0, T ] × R d -→ R d × M d be
uniformly bounded, and continuous in t, Lipschitz in x, we define X by

X 0 := x 0 and d X t = μ(Θ t , t, X t )dt + σ(Θ t , t, X t )dW t , (5.9)
with Θ t := (T Nt , X T N t ). In other words, the process X is defined recursively by, X 0 = x 0 and for all k ≥ 0,

X T k+1 = X T k + T k+1 T k μ T k , X T k , s, X s ds + T k+1 T k σ T k , X T k , s, X s dW s .
Example 5.4. (i) Let (μ, σ)(s, y, t, x) = (µ, σ)(s, y), then X is defined as a Euler scheme as in (3.2), i.e. X 0 = x 0 , and

X T k+1 = X T k + µ(T k , X T k )∆T k+1 + σ(T k , X T k )∆W T k+1 .
(ii) When μ(•) ≡ 0 and σ(s, y, t, x) = σ(s, y) + ∂ x σ(s, y)(xy), then SDE (5.9) turns to be a linear SDE, whose solution is given explicitly in (3.8).

We first formulate an assumption on the existence of automatic differentiation weights associated to SDE (5.9)

. Let θ ∈ [0, T ) × R d and (t, x) ∈ [0, T ] × R d , the process ( X t,x,θ s ) s∈[t,T ] is defined by SDE X t,x,θ t := x, d X t,x,θ s = μ θ, s, X t,x,θ s ds + σ θ, s, X t,x,θ s dW s , (5.10)
Assumption 5.5. There is a pair of measurable functions

W 1 θ (•), W 2 θ (•) , called auto- matic differentiation weights, taking values in R d ×M d , such that, for all θ ∈ [0, T )×R d , (t, x) ∈ [0, T ) × R d , s > t, one has W i θ t, x, s -t, (W r -W t ) r∈[t,s]
), i = 1, 2 are both integrable. Moreover, for all bounded continuous function φ : R d → R,

D i E φ X t,x,θ s = E φ X t,x,θ s W i θ t, x, s -t, (W r -W t ) r∈[t,s] , i = 1, 2,
where D, D 2 denote the gradient and Hessian operators with respect to the variable x.

Let a(•) := 1 2 σσ ⊤ (•) and â(•) := 1 2 σσ ⊤ (•), we denote Θ 0 = (t, x) and then Θ k = (T t k , X T k ), for all k > 0.

and then for k > 0,

∆f k := (µ, a) T k , X T k -(μ, â) Θ k-1 , T k , X T k ∈ R d × M d ,
and for k ≥ 0,

W k := W 1 Θ k , W 2 Θ k T k , X T k , T k+1 , ∆W k+1 • ∈ R d × M d ,
with the weight functions W 1 θ (•), W 2 θ (•) given in Assumption 5.5. We then define

ψ := e βT g( X T ) -g( X T N T )1 {N T >0} β -N T N T k=1 ∆f k • W k , (5.11)
where (p, P )

• (q, Q) := p • q + P : Q for all p, q ∈ R d , P, Q ∈ M d .
Here we use the convention Π 0 k=1 = 1. Finally, for all n ≥ 1, we also introduce

ψ n = e βT n+1 N T ∧n k=1 β -1 ∆f k • W k g X T -g X T N T 1 {N T >0} 1 {N T ≤n} + β -1 ∆f n+1 • Du, D 2 u T n+1 , X T n+1 1 {N T >n} . (5.12) Assumption 5.6.
(i) The sequence (ψ n ) n≥0 is uniformly integrable. (ii) Let (e i ) i=1,••• ,d denote the canonical basis of R d .
There is some ε 0 > 0, such that for all (t, x) ∈ [0, T ) × R d and θ ∈ [0, T ) × R d , n ≥ 0 and i = 1, • • • , d, all the following random vectors is integrable:

W 1 θ (t, x, τ 1 ∧ (T -t), (W r -W t ) r∈[t,(t+τ 1 )∧T ] ), sup ε∈(0,ε 0 ] 1 ε W 1 θ t, x + εe i , τ 1 ∧ (T -t), (W • -W t ) -W 1 θ t, x, τ 1 ∧ (T -t), (W • -W t ) and ∆f n+1 • Du, D 2 u T n+1 , X T n+1 W n .
Theorem 5.7. Suppose that the PDE (5.7)

has a classical solution u ∈ C 1,3 b [0, T ] × R d ,
suppose in addition that Assumptions 5.5 and 5.6 hold true. Then ψ is integrable and u(0, x 0 ) = E ψ . (ii) By definition, the automatic differentiation weight satisfies E W k = 0, then ψ in (5.11) has the same mean than the estimator

e βT g( X T ) β -N T N T k=1 ∆f k • W k .
However, in practice, the weight function W k is typically of infinity variance, or even not integrable, in general. Indeed, as we will see in the following, W k is generally of order

1 ∆T k+1 = 1 T k+1 -T k , where conditioning on N T = n, (T 1 , • • • , T N T )
follows the law of statistic order of uniform distribution on [0, T ]. Then by direct computation, one knows E 1/∆T N T +1 = ∞. In the definition of ψ in (5.11), the additional term -g X T N T 1 {N T >0} can be seen as a control variate so as to guarantee the integrability of ψ. (iii) As a consequence of the integrability problems raised in (ii), Assumption 5.6 is in fact implicitly a restriction on the choice of the coefficients μ and σ, and we cannot expect a representation for u(t, x) with arbitrary μ and σ, see Section 5.3 below.

Proof of Theorem 5.7. (5.7). Denote (μ θ , âθ )(•) = (μ, â)(θ, •), one can rewrite (5.7) in the following equivalent way:

(i) Recall that u ∈ C 1,3 b [0, T ] × R d is a classical solution of PDE
-∂ t u -μθ • Du -âθ : D 2 u -(µ -μθ ) • Du + (a -âθ ) : D 2 u = 0 .
(5.13)

Using Feynmann-Kac formula, it follows that

u(0, x 0 ) = E g X 0,x 0 ,θ T + T 0 (µ -μθ ) • Du + (a -âθ ) : D 2 u s, X 0,x 0 ,θ s ds , (5.14)
where X 0,x 0 ,θ is defined by (5.10), which coincides with X in (5.9) on [0, T 1 ] whenever θ = (0, x 0 ). Recall that T 1 = τ 1 ∧ T , where τ 1 is a random variable of density βe -βs 1 {s≥0} independent of the Brownian motion W . Fixing θ = (0, x 0 ), it follows that

u(0, x 0 ) = E e βT 1 g X T 1 {N T =0} + β -1 ∆f 1 • (Du, D 2 u) T 1 , X T 1 1 {N T >0} = E[ ψ 0 ].
(ii) Let us now go back to the expression (5.14), and derive an expression for the derivatives Du(0, x 0 ) and D 2 u(0, x 0 ). First, for Du(0, x 0 ), we use the integrability condition in Assumption 5.6 with Lemma A.3, and also the fact the Du(•) is continuous, it follows that

Du(0, x 0 ) = E g X 0,x 0 ,θ T W 1 θ (x 0 , T ) + T 0 (µ -μθ ) • Du + (a -âθ ) : D 2 u s, X 0,x 0 ,θ s W 1 θ (x 0 , s)ds ,
where we simplify the notation W 1 θ (0, x 0 , s, (W r -W t ) r∈[0,s] ) to W 1 θ (x 0 , s). Then by the independence of τ 1 to the Brownian motion W , and setting θ = (0, x 0 ), it follows that

Du(0, x 0 ) = E ψ 0 W 1 (0,x 0 ) (0, x 0 , T 1 , ∆W 1 • ) . (5.15)
Next, for D 2 u(0, x 0 ), we use again Lemma A.3 together with the integrability condition in Assumption 5.6 and Lipschitz property of x → ((µμ) • Du + (aâθ ) : D 2 u)(s, x), and the continuity of D 2 u(•) that

D 2 u(0, x 0 ) = D 2 x 0 E g X 0,x 0 ,θ T + T 0 D 2 x 0 E (µ-μ θ )•Du+(a-â θ ) : D 2 u (s, X 0,x 0 ,θ s ) ds.
Setting θ = (0, x 0 ) and using Assumption 5.5, it leads to

D 2 u(0, x 0 ) = E ψ 0 W 2 (0,x 0 ) (0, x 0 , T 1 , ∆W 1 • ) .
Recall that E[ W 2 (0,x 0 ) (0, x 0 , T, ∆W • )] = 0, we then obtain

D 2 u(0, x 0 ) = E ψ 0 -e βT g(x 0 )1 {N T =0} W 2 (0,x 0 ) (0, x 0 , T 1 , ∆W 1 • ) . (5.16)
(iii) Changing the initial condition (0, x 0 ) in (5.15) and (5.16) by (T 1 , X 1 ) (remember ψ 0 dependent also on the initial condition (0, x 0 )), then plugging the expression of D 1 u(T 1 , X 1 ) and D 2 u(T 1 , X 1 ) into the definition of ψ 0 in (5.12), it follows by identifying the term that

u(0, x 0 ) = E[ψ 1 ].
(iv) Repeating the arguments by replacing the initial condition (0, x 0 ) by (T n+1 , X n+1 ) in (5.15) and (5.16) and then plugging the corresponding expression into the definition of ψ n , etc., we obtain that u(0, x 0 ) = E[ ψ n ] for all n ≥ 0. Then letting n -→ ∞, we obtain

u(0, x 0 ) = lim n→∞ E ψ n = E lim n→∞ ψ n = E ψ ,
which concludes the proof.

5.3 Proof of the representation results in Theorems 2.2, 2.5 and 3.3.

Using the results in Theorem 5.7, we can easily complete the proof of the representation results in Theorems 2.2, 2.5 and 3.3.

Proof of Theorems 2.2 (ii) and 3.3 (ii). (i) In the context of Theorems 2.2 and 3.3, the increment X T k+1 -X T k , conditional on (T k , X T k ), is Gaussian. And the estimator ψ corresponds to the estimator in Theorem 5.7 with automatic differentiation weights function

W 1 θ •, δt, δw := (σ ⊤ 0 ) -1 δw δt and W 2 θ •, δt, δw := (σ ⊤ 0 ) -1 δwδw ⊤ -δtI d δt 2 σ -1 0 .(5.17)
In particular, it is clear that Assumption 5.5 holds true with the above choice of automatic differentiation weight functions in (5.17).

(ii) Besides, the uniform integrability conditions and integrability conditions in Assumption 5.6 can be easily obtained following the lines in the first part of the proof of Theorems 2.2 and 2.2, using Lemma A.2.

(iii) Now, suppose in addition that µ, σ and g are bounded smooth functions with bounded continuous derivatives, so that u ∈ C (iv) Finally, when µ(•) and σ(•) satisfy the Lipschitz condition (3.1) and g is Lipschitz, we can find a sequence of bounded smooth functions (µ ε (•), σ ε (•), g ε (•)) which converges locally uniformly to (µ(•), σ(•), g(•)) as ε → 0. Let X ε be the solution of

dX ε t = µ ε (t, X ε t )dt + σ ε (t, X ε t )dW t .
Then by the stability of SDEs together with dominated convergence theorem, it follows that

V ε 0 := E g ε (X ε T ) -→ V 0 := E g(X T ) , as ε → 0.
Moreover, by Lemma A.2 together with dominated convergence theorem, it is easy to prove that E[ ψ ε ] → E[ ψ] as ε → 0, where ψ ε denotes the estimator of the algorithm (3.3) associated to the coefficient (µ ε , σ ε , g ε ). We then conclude the proof.

Proof of Theorem 2.5 (ii). For the path-dependent case, it is enough to use the same arguments as in Theorem 2.2, together with the PDE system (2.12) in Remark 2.3.

Proof of Theorem 3.5

To introduce the algorithm in the context of Theorem 5.7, we propose to choose μ(•) ≡ 0 and σ(s, y, t, x) = σ(s, y) + ∂ x σ(s, y)(xy).

Before providing the proof of Theorem 3.5, we first give a lemma which justifies our choice of the automatic differentiation weight function W 2 k in (3.11), as well as some related estimations. Let c 1 , c 2 , x ∈ R be constants such that c 1 + c 2 x = 0, we denote by X 0,x solution of the SDE

X 0 = x, dX t = c 1 + c 2 X t dW t , (5.18)
whose solution is given explicitly by

X 0,x t =    -c 1 c 2 + c 1 c 2 + x exp - c 2 2 2 t + c 2 W t , if c 2 = 0, x + c 1 W t , if c 2 = 0. (5.19)
Consider also its antithetic variable X x t defined by

X 0,x t =    -c 1 c 2 + c 1 c 2 + x exp - c 2 2 2 t -c 2 W t , if c 2 = 0, x -c 1 W t , if c 2 = 0. Lemma 5.9. Let x ∈ R, (c 1 , c 2 ) ∈ R 2 be two constants such that c 1 + c 2 x = 0, φ : R → R a bounded continuous function. (i) Then for all t ∈ (0, T ], ∂ 2 xx E φ X 0,x t = E φ X 0,x t 1 (c 1 + c 2 x) 2 -c 2 W t t + W 2 t -t t 2 . (5.20) (ii) Suppose in addition that φ(•) ∈ C 2 b (R).
Then there is some constant C independent of (t, x) such that, for all (t, x) ∈

[0, T ] × R d , E φ X 0,x t -φ(x) 2 W t t 2 + E φ X 0,x t -2φ(x) + φ( X 0,x t ) 2 W 2 t -t t 2 2 ≤ C(c 1 + c 2 x) 2 .
Proof. (i) First, when c 2 = 0, it is clear that result is correct (see e.g. Lemma 2.1 of Fahim, Touzi and Warin [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF]). Next, when c 2 = 0, denote v(x) := E φ X 0,x t , then with the expression of X 0,x t in (5. [START_REF] Henry-Labordère | Counterparty Risk Valuation: A Marked Branching Diffusion Approach[END_REF], it follows that

v(x) = R φ - c 1 c 2 + c 1 c 2 + x e -c 2 2 t/2+c 2 √ ty 1 √ 2π e -y 2 /2 dy.
Suppose that φ(•) ∈ C2 b (R), then using integration by parts, it follows that

v ′ (x) = R φ ′ - c 1 c 2 + c 1 c 2 + x e -c 2 2 t/2+c 2 √ ty e -c 2 2 t/2+c 2 √ ty 1 √ 2π e -y 2 /2 dy = R φ - c 1 c 2 + c 1 c 2 + x e -c 2 1 c 1 + c 2 x y √ t 1 √ 2π e -y 2 /2 dy = E φ(X 0,x t ) 1 c 1 + c 2 x W t t .
Similarly, still using integration by parts, and by direct computation, we obtain

v ′′ (x) = E φ X 0,x t 1 (c 1 + c 2 x) 2 -c 2 W t t + W 2 t -t t 2 .
When φ(•) is only a bounded continuous function, one can approximate φ(•) by a sequence of smooth function φ ε (•) which converges to φ(•) uniformly, and φ ′ ε and φ ′′ ε are bounded continuous. We then obtain

v ε (x) := E φ ε X 0,x t → v(x).
Moreover, the limit lim ε→0 v ′ ε (x), lim ε→0 v ′′ ε (x) exist, thus v ′′ (x) also exists and

v ′′ (x) = lim ε→0 v ′′ ε (x) = E φ X 0,x t 1 (c 1 + c 2 x) 2 -c 2 W t t + W 2 t -t t 2 .
(ii) When c 2 = 0, the estimation in (ii) of the statement is clear true since φ ′ and φ ′′ are uniformly bounded. When c 2 = 0, denote |φ ′ | 0 := sup x |φ ′ (x)|, we obtain, by direct computation, that

E φ X 0,x t -φ(x) 2 W t t 2 ≤ |φ ′ | 0 E X 0,x -x 2 W 2 t t 2 = |φ ′ | 0 E c 1 + c 2 x 2 e -c 2 2 t/2+c 2 Wt -1 c 2 W t -c 2 2 t/2 2 W 2 t (c 2 W t -c 2 2 t/2) 2 t 2 ,
which is clearly uniformly bounded by

C(c 1 + c 2 x) 2 for some constant C independent of (t, x) ∈ [0, T ] × R d .
Next, denote ℓ(y) := x + c 1 c 2 e -c 2 2 t/2+c 2 y -1 , and define ϕ(y) := φ(x + ℓ(y)). Then

ϕ ′′ (y) = φ ′′ (x + ℓ(y))(c 2 + c 1 x) 2 e -c 2 2 t+2c 2 y + φ ′ (x + ℓ(y))(c 2 + c 1 x)c 2 e -c 2 2 t/2+c 2 y .(5.21)
It follows by the definition of ϕ as well as its derivative, together with direct computation, that

E φ X 0,x t -2φ(x) + φ( X 0,x t ) 2 W 2 t -t t 2 2 = E ϕ(W t ) + ϕ(-W t ) -2ϕ(0) 2 W 2 t -t t 2 2 + E 2 ϕ(0) -φ(x) 2 W 2 t -t t 2 2 ≤ E W 2 t (W 2 t -t) t 2 2 sup |z|≤|Wt| ϕ ′′ (z) + E 2 φ x + c 1 + c 2 x c 2 e -c 2 2 t/2 -1 -φ(x) 2 W 2 t -t t 2 2 ,
which is also uniformly bounded by C(c 1 + c 2 x) 2 for some constant C > 0,

Proof of Theorem 3.5. (i) Let us first prove that E ψ 2 < ∞ for ψ defined by (3.12). First, we notice that

W - k = W 2 k for all k = 1, • • • , N T -1, g ∈ C 2 b (R), and with the choice of c k 1 and c k 2 in (3.7), one has c k 1 + c k 2 X T k = σ(T k , X T k )
, which is uniformly bounded. By considering the conditional expectation over ( X T N T , ∆T N T +1 ) using items (ii) of Lemma 5.9, we have E ψ 2 is bounded by

CE β -2N T N T k=2 a(T k , X T k ) -ãk 2a(T k , X T k ) -∂ x σ(T k , X T k ) ∆W T k ∆T k + ∆W 2 T k -∆T k ∆T 2 k 2 ,
for some constant C. Further, by denoting ∆ X

T k := X T k -X T k-1 , one has a(T k , X T k ) -ãk ≤ |σ| 0 + ∂ x σ(T k-1 , X T k-1 )∆ X T k /2 |∂ t σ| 0 ∆T k + ∂ 2 xx σ 0 ∆ X T k 2 ,
where |σ| 0 := sup t,x |σ(t, x)|. Notice that σ ≥ ε > 0, σ and ∂ x σ are uniformly bounded, then to prove that ψ is of finite variance, it is enough to prove that, for some C > 0 large enough, the expectation of

N T k=2 C C + ∂ x σ(T k-1 , X T k-1 )∆ X T k 2 C + ∆ X 2 T k ∆T k 2 C ∆W T k + ∆W 2 T k ∆T k + 1 2 (5.22)
is finite. Similarly to the computation in item (ii) of Lemma 5.9, we have

∆ X T k = X T k -X T k-1 = σ(T k-1 , X T k-1 ) exp -∂ x σ(T k-1 , X T k-1 ) 2 ∆T k /2 + ∂ x σ(T k-1 , X T k-1 )∆W T k -1 ∂ x σ(T k-1 , X k-1
) .

Notice again that σ(•) and ∂ x σ(•) are uniformly bounded, it follows that

E C + ∂ x σ(T k-1 , X T k-1 )∆ X T k C + ∆ X 2 T k ∆T k C ∆W T k + ∆W 2 T k ∆T k + 1 2 X T k-1 , T k-1 , ∆T k ≤ C ′ , for some constant C ′ > 0 independent of X T k-1 , T k-1 , ∆T k .
Then the variance of (5.22) is bounded by CE (C ′ ) N T < ∞ and hence ψ in (3.12) is of finite variance.

(ii) Let us now consider the estimator ψ. By the same computation, we obtain that

E ψ N T , ∆T 1 , • • • , ∆T N T +1 ≤ C N T 1 ∆T N T +1
, for some C > 0, where the r.h.s. is integrable but of infinite variance (see Lemma A.2). Similarly, it is easy to check the uniform integrability condition in Assumption 5.6 for ψ in (3.10).

(iii) Finally, using Lemma 5.9 (i) , it follows that Assumption 5.5 holds true. Moreover, with the regularity condition on σ(t, x) and g in Assumption 3.4, we know u ∈ C 1,3 b (R). We then conclude the proof of u(0, x 0 ) = E[ ψ] = E[ψ] by Theorem 5.7.

Let N = (N s ) s≥0 be a Poisson process with arrival times (T k ) k>0 , denote ∆T k+1 := T k+1 -T k . Let 0 = t 0 < t 1 < • • • < t n = T < ∞ be a discrete time grid, we define further Tk := min(T k , t i ), whenever T k-1 ∈ [t i-1 , t i ) for some i = 1, • • • , n, and ∆ Tk := Tk -Tk-1 for every k = 2, • • • , N T + 1.

Lemma A.2. Let p ∈ (0, 1), then for every constant C > 0, one has

E N T k=1 C (∆ Tk+1 ) p < ∞.
Proof. (i) Notice that we can always add points into the time grid 0 = t 0 < t 1 < • • • < t n = T , which makes ∆ Ti+1 smaller. Therefore, one can suppose without loss of generality that t kt k-1 < (2βC) p-1 for every k = 1, • • • , n.

E   N k i=1 C (∆ T k i+1 ) p   ≤ e -β(t k -t k-1 ) ∞ m=0 (β(t k -t k-1 )) m m! m!2 m C (t k -t k-1 ) p m = e -β(t k -t k-1 ) ∞ m=0 2βC(t k -t k-1 ) 1-p m < ∞,
where the last inequality follows by the fact t kt k-1 < (2βC) p-1 . We then conclude the proof by (A.2). Let X x be the solution of SDE X x 0 = x, dX x t = µ(t, X where, in particular, the partial derivative at the l.h.s. and the integration at the r.h.s. are well defined.

Proof. (i) First, let us notice that (t, x) → (µ, σ)(t, x) is Lipschitz in x, then by standard analysis (see e.g. Chapter 7.8 of [START_REF] Graham | Stochastic Simulation and Monte Carlo Methods[END_REF]), there is some constant C independent of ε > 0 and i = 1 where the first equality follows by the Lipschitz property of x → φ(t, x) and the estimation (A.6) together with (A.4), and in particular, the integrable in the last term of (A.7) is well defined, and hence the limit of the first term of (A.7) exists.

Figure 1 :

 1 Figure 1: Comparison of the computation time of MLMC method and unbiased simulation method for problem (4.3), with the same given error.

. 1 ) 5 . 1 .

 151 Proposition Let µ(•, •) and g(•) be both bounded smooth functions in C 2 b

. 7 )

 7 with terminal condition u(T, x) = g(x), where a(•) := 1 2 σσ ⊤ (•), A : B := Tr(AB ⊤ ) for any two d × d dimensional matrices A, B ∈ M d , and D, D 2 denote the gradient and Hessian operators with respect to the space variable x. Next, let us consider the diffusion process (X 0,x 0 s

Remark 5 . 8 .

 58 (i) The condition that u ∈ C 1,3 b [0, T ] × R dmay be relaxed in the concrete applications of Theorem 5.7. This will be indeed performed in Section 3.3 by exploiting the integrability of the automatic differentiation weights W 1 θ , W 2 θ of Assumption 5.5.

  1,3 b ([0, T ] × R d ). It follows by Theorem 5.7 that V 0 = E[ ψ].

(

  ii) For every k = 1, • • • , n, we denote N k := #{i : T i ∈ [t k-1 , t k )}, and T k i := Tk i with k i := j<k N j + i for i = 1, • • • , N k + 1, and ∆ T k i := T k i -T k i-1 . By the memoryless property of the exponential distribution, it is clear that ∆ T 1 i , i = 2, • • • N 1 + 1 , • • • , ∆ T n i , i = 2, • • • N n + 1 are mutually independent. Moreover, we have law of T k i , i = 1, • • • , N k conditioning on N k = mis the law of order statistics of uniform distribution on [t k-1 , t k ]. Then it follows by Lemma A.1 that for every k = 1, • • • , n,

  x) is bounded and Lipschitz in (t, x), ∂ x σ(t, x) is bounded continuous in (t, x) and uniformly Lipschitz in x. Further, the terminal condition function g

  c 4 is linear operator from R d to M d . However, to the best of our knowledge, the exact simulation of linear SDE (3.14) in high dimensional case, as well as the associated automatic differentiation (Malliavin) weight as in (3.11) (see also Lemma 5.9 below), is still an open question.

Table 1 :

 1 Numerical results for V 0 in (4.2) (case d = 1), US denotes our unbiased simulation algorithm with β = 0.1, the computation times are expressed in second.

	1 -1/8 dt +	1 2	dW t .	(4.1)

Table 2 :

 2 Numerical results for V 0 in (4.2) (case d = 1), US denotes our unbiased simulation algorithm with β = 0.05, the computation times are expressed in second.

		Mean value Statistical error Computation time
	US (N = 10 5 )	0.127032	0.000762635	0.144998
	MLMC	0.127053	0.000536248	0.323337
	US (N = 10 6 )	0.126363	0.000241231	1.40843
	MLMC	0.126747	0.000169842	1.8194
	US (N = 10 7 )	0.126703	7.6418e-05	13.9005
	MLMC	0.126643	5.37691e-05	16.7499

Table 3 :

 3 Numerical results for V 0 in (4.4) (case d = 4), US denotes our unbiased simulation algorithm with β = 0.5, the computation times are expressed in second.

	1 and

  Lemma A.3. Suppose that for all t ∈ [0, T ] and bounded continuous function φ :[0, T ] × R d → R, the derivatives ∂ x i E φ(t, X x t ) , ∂ 2 x i ,x j E φ(t, X x t ) i,j=1,••• ,d exist;and there is some measurable R d -valued function W 1 x, t, (W s ) s∈[0,t] such that for∂ x i E φ(t, X x t ) = E φ(t, X x t ) W 1 i (x, t, (W s ) s∈[0,t] ) , i = 1, • • • , d.Let F (dt) be some probability measure on [0, T ].(i) Suppose that for eachx ∈ R d , i = 1, • • • , d, (W s ) s∈[0,t] ) F (dt) < ∞,and the continuous function φ : [0, T ] × R d → R (which may be unbounded) satisfies(x, t, (W s ) s∈[0,t] ) F (dt) < ∞.for some ε 0 > 0, where (e j ) j=1,••• ,d denotes the canonical basis of R d . Then

	T i (x, t, T 0 E W 1 E φ(t, X x t ) W 1 i Then 0
		∂ x i	0	T	E φ(t, X x t ) F (dt) =	0	T	E φ(t, X x t ) W 1
	and
	0	T	sup ε∈[0,ε 0 ]	1 ε	E W 1
					∂ 2 x i x j	0	T	E φ(t, X x t ) F (dt) =	0	T	∂ 2 x i x j E φ(t, X x t ) F (dt),	(A.5)

x t )dt + σ(t, X x t )dW t , where (µ, σ) : [0, T ]×R d → R d ×M d is continuous and in addition Lipschitz continuous in x. i (x, t, (W s ) s∈[0,t] ) F (dt). (A.3) (ii) Suppose in addition that φ(t, x) is bounded continuous function and Lipschitz in x, and for each x

∈ R d , i, j = 1, • • • d, T 0 E W 1 i (x, t, (W s ) s∈[0,t] ) 2 F (dt) < ∞, i (x + εe j , t, (W s ) s∈[0,t]

) -W 1 i (x, t, (W s ) s∈[0,t] ) F (dt) < ∞, (A.

[START_REF] Bally | The law of the Euler scheme for stochastic differential equations, I: convergence rate of the distribution function[END_REF]

  , • • • d, such that E X x+εe i -X x εwhere the first equality follows by the Lipschitz property of x → φ(t, x) and (A.6). We hence proved (A.3) when x → φ(t, x) is Lipschitz.(iii) When φ is only continuous, it is enough to approximate it by a sequence (φ n ) n≥1 which are all bounded, and Lipschitz in x. Then by the integrability of W 1 i (x, t, (W s ) s∈[0,t] ) as well as that of φ(t, X x t ) W 1 i (x, t, (W s ) s∈[0,t] ) under P(dω)×F (dt), it follows that (A.3) holds true for continuous function φ. (x + εe i , t, •)φ(t, X x t ) W 1 j (x, t, •) F (dt). (A.7)

	(iv) To prove (A.5), let us use (A.3) and obtain that
		lim ε→0	1 ε	∂ x j	0	T	E φ(t, X x+εe i t	) F (dt) -∂ x j	0	T	E φ(t, X x t ) F (dt)
	= lim ε→0	0	T	1 ε	E φ(t, X x+εe i t	) W 1 j (x + εe i , t, •) -φ(t, X x t ) W 1 j (x, t, •) F (dt)
	=	0	T	lim ε→0	1 ε	E φ(t, X x+εe i t	) W 1 j
											2	≤ C(1 + e Ct ).	(A.6)
	(ii) Suppose that φ(t, x) is bounded continuous and Lipschitz in x. It follows that
							lim ε→0	0	T	1 ε	E φ(t, X x+εe i t	) -φ(t, X x t ) F (dt)
						=	0	T	lim ε→0	1 ε	E φ(t, X x+εe i

t )φ(t, X x t ) F (dt) = E T 0 φ(t, X x t ) W 1 i (x, t, (W s ) s∈[0,t] )F (dt) ,

t/2+c 2 √ ty

referees for valuable remarks and suggestions. X. Tan and N. Touzi gratefully acknowledge the financial support of the ERC 321111 Rofirm, the ANR Isotace, and the Chairs Financial Risks (Risk Foundation, sponsored by Société Générale) and Finance and Sustainable Development (IEF sponsored by EDF and CA).

A Appendix

We first provide an estimation on the order statistics of uniform distribution on [0, 1], which induces an estimation on a functional of the arrival times (T k ) k>0 of the Poisson process. We next provide a technical result on the automatic differentiation function related to a SDE.

Proof. First, we notice that for any x ∈ (0, 1),

Then, since the density of the order statistics (U

it follows by direct computation that

where the last inequality follows from (A.1).