Pierre Henry-Labordère

Nadia Oudjane
email: nadia.oudjane@edf.fr

Xiaolu Tan
email: tan@ceremade.dauphine.fr

Nizar Touzi
email: nizar.touzi@polytechnique.edu

Xavier Warin
email: xavier.warin@edf.fr

Branching diffusion representation of semilinear PDEs and Monte Carlo approximation *

Keywords: Semilinear PDEs, branching processes, Monte-Carlo methods

, based on the Malliavin integration by parts, which allows to account for the nonlinearities in the gradient. As a consequence, the particles of our branching diffusion are marked by the nature of the nonlinearity. This new representation has very important numerical implications as it is suitable for Monte Carlo simulation. Indeed, this provides the first numerical method for high dimensional nonlinear PDEs with error estimate induced by the dimension-free Central limit theorem. The complexity is also easily seen to be of the order of the squared dimension. The final section of this paper illustrates the efficiency of the algorithm by some high dimensional numerical experiments.

Introduction

The objective of the present paper is to provide a probabilistic representation for the solution of a nonlinear parabolic second order partial differential equation (PDE) which is suitable for a high dimensional Monte Carlo approximating scheme. Our main results achieve this goal in the context of semilinear PDEs:

-∂ t u -Lu = f (u, Du), u T = g, t < T, x ∈ R d ,
with polynomial non-linearity f t,x (y, z) in the solution and its gradient, diffusion generator L, and bounded terminal condition g.

Previous representation results were obtained in the literature by means of backward stochastic differential equations, as introduced by Pardoux and Peng [START_REF] Pardoux | Adapted solutions of backward stochastic differential equations[END_REF]. The Monte Carlo numerical implications of this representation were introduced by Bally & Pagès [START_REF] Bally | Error analysis of the quantization algorithm for obstacle problems[END_REF], Bouchard & Touzi [START_REF] Bouchard | Discrete-time approximation and Monte-Carlo simulation of backward stochastic differential equations[END_REF] and Zhang [START_REF] Zhang | A numerical scheme for backward stochastic differential equations[END_REF], and generated a large stream of the literature. However, these methods can be viewed as a Monte Carlo version of the finite elements methods, and as such, are subject to the problem of curse of dimensionality. Our primary goal is to avoid this numerical problem so as to be capable to handle genuinely high-dimensional problems. This however will be achieved at the cost of some limitations...

Our main representation result is obtained by using the branching diffusion trick to absorb the nonlinearity, as illustrated by Skorokhod [START_REF]Skorokhod Branching diffusion processes[END_REF], Watanabe [START_REF] Watanabe | On the branching process for Brownian particles with an absorbing boundary[END_REF] and McKean [START_REF] Mckean | Application of Brownian motion to the equation of Kolmogorov-Petrovskii-Piskunov[END_REF] in the context of the KPP equation, see also the extensions in Rasulov, Raimova & Mascagni [START_REF] Rasulov | Monte Carlo solution of Cauchy problem for a nonlinear parabolic equation[END_REF] and our previous paper [START_REF] Henry-Labordère | A numerical algorithm for a class of BSDEs via the branching process[END_REF] where the representation is also shown to allow for path-dependency.

Since the gradient is also involved in the nonlinearity, our representation result is a significant improvement of the classically well-know representation of KPP equations. We observe that the polynomial nonlinearity naturally induces some restrictions needed to ensure the non-explosion of the corresponding solution. As a consequence, our representation holds under technical conditions of small maturity or small nonlinearity of the PDE.

The main idea for our representation is to use the Monte Carlo automatic differentiation technique in addition to the branching diffusion representation. The automatic differentiation in Monte Carlo approximation of diffusions was successfully used in the previous literature by Fournié et al. [START_REF] Fournié | Applications of Malliavin Calculus to Monte Carlo Methods in Finance[END_REF], Bouchard, Ekeland & Touzi [START_REF] Bouchard | On the Malliavin approach to Monte Carlo approximation of conditional expectations[END_REF], Henry-Labordère, Tan & Touzi [START_REF] Henry-Labordère | Unbiased simulation of stochastic differential equations[END_REF], and Doumbia, Oudjane & Warin [START_REF] Doumbia | Computing expectations for general SDE with pure Monte Carlo methods[END_REF]. The resulting branching diffusion in the representation differs from that of the original founding papers [START_REF]Skorokhod Branching diffusion processes[END_REF][START_REF] Watanabe | On the branching process for Brownian particles with an absorbing boundary[END_REF][START_REF] Mckean | Application of Brownian motion to the equation of Kolmogorov-Petrovskii-Piskunov[END_REF] by introducing marks for the particles born at each branching. The mark of the particle determines the nature of the differentiation, and thus induces the corresponding automatic differentiation weight.

We next illustrate the main idea behind our representation in the context of the following extension of the one-dimensional Burgers equation:

Lu := 1 2
∆u, and f (y, z) = 1 2 (y 2 + yz).

Let W 1 be a Brownian motion, and τ 1 an independent random variable with density ρ > 0 on R + , and denote F (t) := ∞ t ρ(s)ds. We also introduce another independent random variable I 1 which takes the values 0 and 1 with equal probability. Then, denoting by E t,x the expectation operator conditional on the starting data W t = x at time t, we obtain from the Feynman-Kac formula the representation of the solution u as:

u(0, x) = E 0,x F (T ) g(W T ) F (T ) + T 0 f (u, Du)(t, W t ) ρ(t) ρ(t)dt = E 0,x φ T (1) , W 1 T (1) ,
where T (1) := τ 1 ∧ T , and φ(t, y) := 1 {t≥T } F (T ) g(y)+ 1 {t<T } ρ(t) (uD I 1 u)(t, y).

(1.1)

We next consider the two alternative cases for the value of I 1 .

• On the event set {I 1 = 0}, it follows from the Markov property that:

(uD I 1 u)(t, y) = u(t, y) 2 = E t,y φ(t + τ 1 , W 1 t+τ 1 ) 2 .
The tricky branching diffusion representation now pops up naturally by rewriting the last expression in terms of independent copies (W 1,1 , τ 1,1 ) and (W 1,2 , τ 1,2 ) as:

(uD I 1 u)(t, y) = E t,y φ t + τ 1,1 , W 1,1 t+τ 1,1 E t,y φ t + τ 1,2 , W 1,1 t+τ 1,2 = E t,y φ t + τ 1,1 , W 1,1 t+τ 1,1 φ t + τ 1,2 , W 1,2 t+τ 1,2 .

where E t,y denotes the expectation operator conditional on W 1,1 t = W 1,2 t = y. Substituting this expression in (1.1) and using the tower property, we see that the branching mechanism allows to absorb the nonlinearity.

• On the event set {I 1 = 1}, we arrive similarly to the expression (uD I 1 u)(t, y) = E t,y φ t + τ 1,1 , W 1,1 t+τ 1,1 ∂ y E t,y φ t + τ 1,2 , W 1,2 t+τ 1,2 .

Our main representation is based on the following automatic differentiation:

∂ y E t,y φ t + τ 1,2 , W 1,2 t+τ 1,2 = E t,y W 1,2 t+τ 1,2 -W 1,2 t τ 1,2 φ t + τ 1,2 , W 1,1 t+τ 1,2 ,
which is an immediate consequence of the differentiation with respect to the heat kernel, i.e. the marginal density of the Brownian motion. By the independence of W 1,1 and W 1,2 , this provides:

(uD I 1 u)(t, y) = E t,y W 1,2 t+τ 1,2 -W 1,2
The two previous cases are covered by denoting T (1,i) := T ∧ (τ 1 + τ 1,i ) for i = 0, 1, and introducing the random variable:

W 1 := 1 {I 1 =0} +1 {I 1 =1} ∆W 1,2 T (1,2) ∆T (1,2) , with ∆W 1,2 T (1,2) := W 1,2 T (1,2) -W 1,2 T (1) , ∆T (1,2) := T (1,2) -T (1) , so that u(0, x) = E 0,x 1 {T (1) =T } g(W T ) F (T ) + 1 {T (1) <T } W 1 ρ(T (1) ) 2 i=1 1 {T (1,i) =T } g(W 1,i T ) F (∆T (1,i) ) + 1 {T (1,i) <T } (uD I 1,i u) T (1,i) ,W 1,i T (1,i) ρ(∆T (1,i) )
.

Our main representation result is obtained by iterating the last procedure, and solving the integrability problems which arise because of the singularity introduced by the random variable W 1 . The automatic differentiation, which is the main additional ingredient to the branching diffusion representation, is illustrated in the previous example when the operator L corresponds to the Brownian motion. This extends to the case of a more general diffusion operator by the so-called Bismuth-Elworthy-Li formula based on the Malliavin integration by parts formula, see Fournié et al. [START_REF] Fournié | Applications of Malliavin Calculus to Monte Carlo Methods in Finance[END_REF] for its use in the context of Monte Carlo approximation and the extension to other sensitivities.

Our main result provides a probabilistic representation of the solution of the semilinear PDE, with polynomial nonlinearity, in terms of a branching diffusion. This requires naturally a technical condition ensuring the existence of a non-exploding solution for the PDE which can be either interpreted as a small maturity or a small nonlinearity condition. This new representation provides a new ingredient for the analysis of the corresponding PDE as it can be used to argue about existence, uniqueness, and regularity. We shall indeed prove a C 1 -regularity result in order to prove the main Theorem 3.5.

Moreover, our new representation has an important numerical implication as it is suitable for high dimensional Monte Carlo approximation. This is in fact the first high dimensional general method for nonlinear PDEs ! The practical performance of the method is illustrated on a numerical example in dimension d = 20. The convergence of the numerical method is a direct consequence of the law of large numbers. The rate of convergence is also a direct consequence of the central limit theorem, and is therefore dimension-free. The complexity of the method is easily shown to be of the order of d 2 , which cannot be avoided by the very nature of the equation whose second order term involves d × d matrices calculations.

The paper is organized as follows. Section 2 introduces the marked branching diffusion. The main representation result is stated in Section 3. We next provide further discussions in Section 4 on the validity of our representation for systems of semilinear PDEs, and the possible combination with the unbiased simulation technique of [START_REF] Henry-Labordère | Unbiased simulation of stochastic differential equations[END_REF][START_REF] Doumbia | Computing expectations for general SDE with pure Monte Carlo methods[END_REF]. The Monte Carlo numerical implications of our representation in high dimension are reported in Section 5 with an illustration by a numerical experiment in dimension 20. Finally, we provide more numerical examples in Section 6.

b i ) i=1,••• ,m , where c : [0, T ] × R d → R and b i : [0, T ] × R d → R d . For every = ( 0 , 1 , • • • , m ) ∈ L, denote | | := m i=0 i . A generator function f : [0, T ]×R d ×R×R d is then defined by f (t, x, y, z) := =( 0 , 1 ,••• , m)∈L c (t, x) y 0 m i=1 b i (t, x) • z i .
(2.1)

Given two matrix A, B ∈ M d , denote A : B := Trace(AB ). We will consider the following semilinear PDE:

∂ t u + µ • Du + 1 2 σσ : D 2 u + f (•, u, Du) = 0, on [0, T ) × R d , and u(T, .) = g, (2.2)
for some bounded Lipschitz function g : R d -→ R.

Remark 2.1. The nonlinearity (2.1) includes the simplest case of a source term. Indeed, for = (0, 0,

• • • , 0), we have c (t, x)y 0 m i=1 b i (t, x) • z i = c (t, x).

Age-dependent marked branching process

In preparation of the representation result, let us first introduce a branching process, characterized by a distribution density function ρ : R + → R + , a probability mass function (p ) ∈L (i.e. p ≥ 0 and ∈L p = 1). Instead of the usual exponential arrival time, we shall consider a branching particle process with arrival time of distribution density function ρ. At the arrival time, the particle branches into | | offsprings with probability p , among which, i particles carry the mark i, i = 0, . . . , m. Then regardless of its mark, each descendant particle performs the same but independent branching process as the initial particle.

To construct the above process, we will consider a probability space (Ω, F, P) equipped with

• a sequence of i.i.d. positive random variables (τ m,n ) m,n≥1 of density function ρ,

• a sequence of i.i.d. random elements (I m,n ) m,n≥1 with P(I m,n = ) = p , ∈ L.

In addition, the sequences (τ m,n ) m,n≥1 and (I m,n ) m,n≥1 are independent.

We now construct an age-dependent branching process, with (τ m,n ) m,n≥1 and (I m,n ) m,n≥1 , using the following procedure.

1. We start from a particle marked by 0, indexed by (1), of generation 1, whose arrival time is given by T (1) := τ 1,1 ∧ T .

Let

k = (k 1 , • • • , k n-1 , k n ) ∈ N n be
, • • • , k n , i) for i = 1, • • • , |I k |. 3. When I k = ( ˆ 0 , ˆ 1 , • • • , ˆ m ),
we have | ˆ | offspring particles, among which we mark the first ˆ 0 particles by 0, the next ˆ 1 particles by 1, and so on, so that each particle has a mark i

for i = 0, • • • , m. 4. For a particle k = (k 1 , • • • , k n , k n+1 ) of generation n + 1, we denote by k-:= (k 1 , • • • , k n )
the "parent" particle of k, and the arrival time of k is given by

T k := T k-+ τ n+1,π n+1 (k) ∧ T .
5. In particular, for a particle k = (k 1 , • • • , k n ) of generation n, and T k-is its birth time and also the arrival time of k-. Moreover, for the initial particle k = (1), one has k-= ∅, and T ∅ = 0.

The above procedure defines a marked age-dependent branching process. We denote further

θ k := mark of k, K n t := k of generation n s.t. T k-≤ t < T k , when t ∈ [0, T ), {k of generation n s.t. T k = T }, when t = T,
and also

K n t := ∪ s≤t K n s , K t := ∪ n≥1 K n t and K t := ∪ n≥1 K n t .
Clearly, K t (resp. K n t ) denotes the set of all living particles (resp. of generation n) in the system at time t, and K t (resp. K n t ) denotes the set of all particles (resp. of generation n) which have been alive before time t.

Example 2.2. Let us consider the case d = 1, with f (t, x, y, z) := c 0,0 (t, x) + c 1,0 (t, x)y + c 1,1 (t, x)yz.

In this case, m = 1, L = { ¯ 1 = (1, 0), ¯ 2 = (1, 1)}.
For the sake of clarity, we present an typical path of the associated age-dependent process, with graphical illustration below. The process starts from time 0 with one particle indexed by [START_REF] Athreya | Branching processes[END_REF]. At terminal time T , the number of particles alive is 3, with

K T = (1, 2, 1), (1, 1, 1, 1), (1, 1, 1, 2) ,
• At time T (1) , particle (1) branches into two particles (1, 1) and (1, 2).

• At time T (1,1) , particle (1, 1) branches into (1, 1, 1) and (1, 1, 2).

• At time T (1,2) , particle (1, 2) branches into (1, 2, 1).

• At time T (1,1,2) , particle (1, 1, 2) dies out without any offspring particle.

• At time T (1,1,1) , particle (1, 1, 1) branches into (1, 1, 1, 1) and (1, 1, 1, 2).

• The particles in blue are marked by 0, and the particles in red are marked by 1.

E 0 T (1) T (1,1) T (1,2) T (1,1,2) T (1,1,1) T ¨¨¨¨¨¨r r r r r r r r r r r 3 3 3 3 3 3 3 3 3 3 3 3 (1) (1, 2)
(1, 1)

(1, 1 , 1) (1, 1, 2) (1, 2, 1) (1, 1, 1, 1) (1

The marked branching diffusion

We next equip each particle with a Brownian motion in order to define a branching Brownian motion.

We consider a sequence of independent d-dimensional Brownian motion (W m,n ) m,n≥1 , which are also independent of (τ m,n , I m,n ) m,n≥1 . Define W (1) t

= ∆W

(1) t := W 1,1 t for all t ∈ 0, T (1) and then for each k

= (k 1 , • • • , k n ) ∈ K T \ {(1)}, define W k t := W k- T k-+ ∆W k t-T k-, with ∆W k t-T k-:= W n,πn(k) t-T k-, for all t ∈ [T k-, T k ]. (2.3) Then (W k • ) k∈K T is a branching Brownian motion. For each k ∈ K T , we define an associated diffusion process (X k t ) t∈[T k-,T k ] by means of the following SDE X k t = X k- T k-+ t T k- µ s, X k s ds + t T k- σ s, X k s dW k s , t ∈ [T k-, T k ], P-a.s., (2.4)
where for particle (1), we fix the initial condition X

(1) 0 = x 0 for some constant x 0 ∈ R d . The well-posedness of the last SDE is guaranteed by standard conditions on the coefficients µ, σ contained in Assumption 3.1.

The process (X k • ) k∈K T is our main marked branching diffusion process. We finally introduce the sub-σ-fields

F 0 := σ τ m,n , I m,n : m, n ≥ 1 , F m := σ W i,n , τ i,n , I i,n : n ≥ 1, i ≤ m , m ≥ 1. (2.5)

The main representation

We shall provide a representation result for a class the semilinear PDEs (2.2) under general abstract conditions. More explicit sufficient conditions are provided later.

Branching diffusion representation of semilinear PDEs

We first collect the conditions on the marked branching diffusion which are needed for our main results.

F (T ) := ∞ T ρ(t)dt > 0. (ii) (µ, σ) : [0, T ] × R d → R d × M d are bounded continuous, and Lipschitz in x. (iii) c : [0, T ] × R d → R and b i : [0, T ] × R d → R d are bounded continuous.
Our next assumption is the key automatic differentiation condition on the underlying diffusion X t,x s defined by

X t,x s = x + s t µ r, X t,x r dr + s t σ r, X t,x r dW r , s ∈ [t, T ], (3.1)
where W is a d-dimensional Brownian motion independent of the branching diffusion.

Assumption 3.2.

There is a measurable functional W(t, s, x,

(W r -W t ) r∈[t,s] ) sat- isfying (t, x) → W(t, s, x, (W r -W t ) r∈[t,s]
) is continuous, and for any s ∈ [t, T ] and bounded measurable function φ : R d → R, one has

∂ x E φ X t,x s = E φ X t,x s W(t, s, x, (W r -W t ) r∈[t,s] ) . Remark 3.3. In case (µ, σ) ≡ (µ 0 , σ 0 ) for some constant (µ 0 , σ 0 ) ∈ R d × M d
, where σ 0 is not generate, then an example of such automatic differentiation function can be given by

W t, s, x, (W r -W t ) r∈[t,s] := (σ 0 ) -1 W s -W t s -t .
For general coefficient functions (µ, σ) satisfying some regularity and non-degeneracy conditions, one can find such functional W using Malliavin calculus (see more discussions in Section 3.2). Now, for each particle k ∈ K T , we recall that it is born at time T k-and dies out at time T k , its mark is given by θ k and its branching type is given by I k . Let us denote

W k := 1 {θ k =0} + 1 {θ k =0} b θ k (T k-, X k T k-) • W T k-, T k , X k T k-, ∆W k • . (3.2)
We next introduce for a smooth function u

∈ C 1,2 ([0, T ] × R d ) ψ n := k∈∪ n j=1 K j T g(X k T ) -g(X k T k-)1 {θ k =0} F (∆T k ) W k k∈∪ n j=1 (K j T \K j T ) c I k (T k , X k T k ) p I k W k ρ(∆T k ) k∈K n+1 T 1 {θ k =0} u + m i=1 1 {θ k =i} b i • Du (T k-, X k T k-) , (3.3)
for all n ≥ 1, and the corresponding limit

ψ := k∈K T g(X k T ) -g(X k T k-)1 {θ k =0} F (∆T k ) W k k∈K T \K T c I k (T k , X k T k ) p I k W k ρ(∆T k ) . (3.4)
Further, notice that the above branching diffusion process (T k , X k • ) k∈K T and random variables ψ, ψ n are defined with initial condition (0, x 0 ) on interval [0, T ]. By exactly the same way, we can define the system with initial condition (t, x) on interval [t, T ], let us denote them respectively by

(T t k , W t,k • , X t,x,k • ) k∈K t T , ψ t,x
, and ψ t,x n . We now provide a first result, under strong regularity conditions, which provides a better understanding of our representation. We emphasize that our main representation result in Theorem 3.5 below will be established under more transparent conditions. Proposition 3.4. Let Assumptions 3.1 and 3.2 hold true. Suppose that the PDE (2.2)

has a solution u ∈ C 1,2 ([0, T ] × R d ) with E T t f (•, u, Du)(s, X t,x s ) ds < ∞, for some (t, x) ∈ [0, T ] × R d . Assume further that (ψ t,x n ) n≥1 is uniformly integrable. Then ψ t,x ∈ L 1 and u(t, x) = E ψ t,x .
Proof. (i) It suffices to consider (t, x) = (0, x 0 ). Since g is bounded, it follows from the integrability condition on the process f (., u, Du) and the Feynma-Kac formula that

u(0, x 0 ) = E 1 F (T ) g X 0,x 0 T F (T ) + T 0 1 ρ(s) f •, u, Du s, X 0,x 0 s ρ(s)ds = E 1 F (T (1) ) g X (1) T 1 {T (1) =T } + 1 ρ(T (1) ) c I (1) p I (1) u I (1),0 m i=1 b i • Du I (1),i T (1) , X (1)
T (1) 1 {T (1) <T } (3.5) = E[ψ 1 ].
(ii) Next, let b 0 ∈ R d be a constant vector, and assume in addition that the random variable

ψ 1 b 0 • W(0, T (1) , x 0 , ∆W (1)
• ) is integrable. Then under Assumptions 3.2,

b 0 • D x u(0, x 0 ) = E ψ 1 b 0 • W 0, T (1) , x 0 , ∆W (1)
• = E ψ 1 - 1 F (T (1) ) g(x 0 )1 {T (1) =T } b 0 • W 0, T (1) , x 0 , ∆W (1)
• , (3.6)
where the first equality follows by Lemma A.3 of [START_REF] Henry-Labordère | Unbiased simulation of stochastic differential equations[END_REF] and the second equality follows from the fact that E W 0, T, x,

(W s ) s∈[0,T ] = 0. (iii) For k ∈ K 2 T , change the initial condition from (0, x 0 ) to (T k-, X k T k-) = (T (1) , X (1)
T (1) ) in formula (3.5) and (3.6). Then, with F 1 defined in (2.5),

u(T k-, X k T k-) = E 1 {k∈K 2 T } g(X k T ) F (∆T k ) + 1 {k∈K 2 T \K 2 T } Ψ k F 1 ,
by the Markov property, and by Assumption 3.2,

Du(T k-, X k T k-) = E 1 {k∈K 2 T } g(X k T ) -g(X k T k-) F (∆T k ) + 1 {k∈K 2 T \K 2 T } Ψ k W(T k-, T k , X k T k-, ∆W k • ) F 1 ,
where

Ψ k := 1 ρ(∆T k ) c I k (T k ,X k T k ) p I k k -=k 1 {θ k =0} u+ m i=1 1 {θ k =i} b i •Du (T k -, X k T k - ) .
Plugging these expressions in the definition of ψ 1 in (3.3), it follows from the integrability of ψ 2 and the tower property of conditional expectations that u(0, x 0 ) = E ψ 2 .

(iv) Iterating this procedure, we see that

u(0, x 0 ) = E ψ n , for all n ≥ 1, and therefore u(0, x 0 ) = lim n→∞ E ψ n = E ψ ,
where the last equality follows by the uniform integrability condition of (ψ n ) n≥1 .

We now state our main representation result under abstract conditions on the automatic differentiation weight function W(•) involving the slight modification of ψ:

ψ := k∈K T g(X k T ) -g(X k T k-)1 {θ k =0 or k=(1)} F (∆T k ) W k k∈K T \K T c I k (T k , X k T k ) p I k W k ρ(∆T k ) , (3.7)
with ψt,x defined by an obvious change of origin. Explicit sufficient conditions for the validity of the next result will be reported in Section 3.2 below.

Theorem 3.5. Let Assumptions 3.1 and 3.2 hold true, and suppose in addition that for all (t, x) ∈ [0, T ] × R d , there is some ε > 0 such that (ψ s,y ) (s,y)∈Bε(t,x) and ψs,y W(s, T s (1) , y, ∆W

u : [0, T ] × R d → R is continuous. Similarly, let us define v i (t, x) := E ψt,x b i (t, x)•W(t, T t (1) , x, ∆W t,(1) • ) = E ψ t,x b i (t, x)•W(t, T t (1) , x, ∆W t,(1) • ) ,
which is also continuous by the uniformly integrability condition.

(ii) Let us define φ : [0, T ] × R d → R by φ I (1) , T (1) , X (1)
T (1) := 1 F (T (1) ) g X (1) T 1 {T (1) =T } (3.8) + 1 ρ(T (1) ) c I (1) p I (1) ψ T (1) ,X (1)
T (1)

I (1),0 m i=1 v i T (1) , X (1)
T (1)

I (1),i 1 {T (1) <T } = E ψ F 1 ,
where F 1 is defined in (2.5). Notice that φ(i, t, x) is continuous in x, then it follows by Assumption 3.2 and Lemma A.3 of [START_REF] Henry-Labordère | Unbiased simulation of stochastic differential equations[END_REF] that

Du(0, x 0 ) = E φ I (1) , T (1) , X (1)
T (1) W(0, T (1) , x 0 , ∆W

• ) = E ψW(0, T (1) , x 0 , ∆W

• ) .

By changing the initial condition from (0, x 0 ) to (t, x) and notice that ) . Moreover, one has

E[W(t, T t (1) , x, ∆W t,(1) • )1 {T t (1) =T } ] = 0, it follows that Du(t, x) = E ψt,x W(t, T t (
v i (t, x) = b i (t, x) • Du(t, x).
(iii) Using the expression in (3.8) and the law of I (1) and T (1) , and with similar arguments as in (3.5), it follows that

u(t, x) = E[ψ t,x ] = E g X t,x T + T t f •, u, Du (s, X t,x s )ds . Let h > 0, denote h h := (t + h) ∧ inf{s > t : |X t,x
s -x| ≥ 1}, then by the flow property of X t,x , one has

u(t, x) = E u(h h , X t,x h h ) + h h t f •, u, Du (s, X t,x s )ds ,
and we may verify by standard arguments that u is a viscosity solution of PDE (2.2).

More explicit sufficient conditions

We now provide some explicit sufficient conditions which guarantee the validity of the conditions of Theorem 3. We first recall the Bismut-Elworthy-Li formula from Malliavin calculus, which was used by Fournié, Lasry, Lebuchoux, Lions and Touzi [START_REF] Fournié | Applications of Malliavin Calculus to Monte Carlo Methods in Finance[END_REF] as an automatic differentiation tool, see also [START_REF] Bouchard | On the Malliavin approach to Monte Carlo approximation of conditional expectations[END_REF], [START_REF] Bouchard | Discrete-time approximation and Monte-Carlo simulation of backward stochastic differential equations[END_REF] and [START_REF] Fahim | A probabilistic numerical method for fully nonlinear parabolic PDEs[END_REF] for subsequent usefulness of the automatic differentiation in the context of the Monte Carlo approximation of nonlinear PDEs. We emphasize that such automatic differentiation function is not unique.

Assumption 3.6. The coefficients µ, σ are bounded continuous, with bounded continuous partial gradients Dµ, Dσ, and σ is uniformly elliptic.

Notice that (X t,x s ) s∈[t,T ] , as defined by (3.1), is completely determined by (t, x, (W s -W t ) s∈[t,T ] ). We then introduce the corresponding first variation process Y :

Y t := I d , dY s = Dµ(s, X t,x s )Y s ds + d i=1 Dσ i (s, X t,x s )Y s dW i s , for s ∈ [t,
W t, s, x, (W r -W t ) r∈[t,s] = σ 0 -1 W s -W t s -t .
With the above choice of automatic differentiation weight function (3.10), we can now derive some upper bounds for random variables (ψ n , n ≥ 1). Recall that L g is the Lipschitz constant of g, denote by

B ∞ 0 (L g ) := {(x 1 , • • • , x d ) ∈ R d : |x i | ≤ L g } and W t,x,s := W t, s, x, (W r -W t ) r∈[t,s] . Then for n ≥ 1, q > 1, we introduce two constants C 1,q and C 2,q by C 1,q := |g| q ∞ ∨ sup 0≤t<s≤T, x∈R d , i=1,••• ,m, b 0 ∈B ∞ 0 (Lg) E b 0 • (X t,x s -x) b i (t, x) • W t,x,s q and C 2,q := sup 0≤t<s≤T, x∈R d , i=1,••• ,m E √ s -t b i (t, x) • W t,x,s q ,
and then

C 1,q := C 1,q F (T ) q-1 , C 2,q := C 2,q sup ∈L, t∈(0,T ] |c | ∞ p t -q 2(q-1)

ρ(t)

q-1 .

Remark 3.9. (i) Under Assumption 3.6, the tangent process Y is defined by a linear SDE, which has finite moment of any order q ≥ 1. Then the two constant C 1,q and C 2,q are both finite. And for all k ∈ K T , one has

max |g| q ∞ , E (D g • ∆X k )W k q F 0 ≤ C 1,q , E ∆T k |W k | q F 0 ≤ C 2,q , (3.11)
where the sub-σ-field F 0 is defined in (2.5).

(ii) Notice that for a random variable N ∼ N (0, 1) and non-negative integer q ≥ 0, one has

E[|N | q ] = 2 q 2 Γ q+1 2 / √ π. Then if (µ, σ) ≡ (0, σ 0 ), for some constant (µ 0 , σ 0 ) ∈ R d × M d
, and W as in Remark 3.8, it follows by direct computation that

C 1,q ≤ |g| q ∞ ∨ sup b 0 ∈B ∞ 0 (Lg) b 0 σ 0 σ 0 b 0 + max i=1,••• ,m b i (σ 0 σ 0 ) -1 b i ∞ q 2 q-1 Γ 2q + 1 2 / √ π, and
C 2,q = max i=1,••• ,m b i σ 0 σ 0 -1 b i q 2 ∞ 2 q 2 Γ q + 1 2 / √ π.
We are now ready for the main explicit sufficient conditions for the validity of the representation Theorem 3.5. Notice that the following conditions can be interpreted either as a small maturity or small nonlinearity restriction. Assumption 3.10. For some q > 1, one of the following two items holds true.

(i) Both C 1,q 1 F (T ) q and sup ∈L,t∈(0,T ] C 2,q |c |∞ p 1 √ tρ(t) q are bounded by 1. (ii) T < ∞ C 1,q C 2,q ∈L |c | ∞ x | | -1 dx. Remark 3.11. (i) To ensure that sup ∈L,t∈(0,T ] C 2,q |c |∞ p 1 √ tρ(t)
q is bounded by 1, it is necessary to choose (p ) ∈L such that |c |∞ p is uniformly bounded, and to choose a density function such that ρ(t) ≥ Ct -1/2 . (ii) To ensure that C 2,q is finite, one needs to choose the density function ρ such that ρ(t) ≥ Ct -q 2(q-1) , and hence it is necessary that q ∈ (2, ∞) so that q 2(q-1) ∈ (

) (t,x)∈[0,T ]×R d is uni- formly integrable. Consequently, u(t, x) := E[ψ t,x ] is a viscosity solution of PDE (2.2). (ii) If Assumption 3.10 holds with some q ≥ 2, then E |ψ t,x | 2 < ∞.
Proof. (i) First, using Proposition 3.7, it is clear that Assumption 3.2 holds true with the choice of automatic differentiation function in (3.10).

(ii) Next, for q ≥ 1, let us introduce

χ q ∞ := k∈K T C 1,q 1 F (∆T k ) q k∈K T \K T C 2,q |c I k | ∞ p I k 1 √ ∆T k ρ(∆T k ) q .
By conditioning on F 0 , it follows from (3.11), together with direct computation, that

E[|ψ| q ] ≤ E[χ q ∞ ] and E ψW(t, T t (1) , x, ∆W t (1) ) q ≤ CE[χ q ∞ ], (3.12)
for some constant depending only on the Lipschitz constant L g . (iii) When Assumption 3.10 (i) holds true for some q > 1, then it is clear that E[|ψ| q ] ≤ 1. Notice that the above argument is independent of the initial condition (0, x 0 ), it follows that ψ t,x , ψt,x W(t, T t (1) , x, ∆W t (1) ) (t,x)∈[0,T ]×R d is uniformly integrable. (iv) When Assumption 3.10 (ii) holds true for some q > 1. Consider the ODE on [0, T ]:

η(T ) = C 1,q , η (t) + ∈L C 2,q c ∞ η(t) | | = 0.
Under Assumption 3.10 (ii), it is clear that the above ODE admits a unique finite solution on [0, T ]. We next introduce a sequence of random variables

χ q n := k∈∪ n j=1 K j T C 1,q F (∆T k ) k∈∪ n j=1 (K j T \K j T ) C 2,q |c I k | ∞ p I k 1 ρ(∆T k ) k∈K n+1 T η(T k-) , and
χ q ∞ := lim n→∞ χ q n = k∈K T C 1,q F (∆T k ) k∈K T \K T C 2,q |c I k | ∞ p I k 1 ρ(∆T k ) .
Then by the same arguments as in the proof of Proposition 3.4, it is easy to check that

η(0) = η(T ) + T 0 ∈L C 2,q |c | ∞ η(t) | | dt = E χ q 1 = E χ q n , ∀n ≥ 1;
and hence by direct computation, it follows that

E ψ q ≤ E χ q ∞ ≤ E χ q ∞ ≤ lim inf n→∞ E χ q n = η(0) < ∞.
Changing the origin from (0, x 0 ) to (t, x), we see that sup

(t,x)∈[0,T ]×R d E |ψ t,x | q ≤ sup t∈[0,T ] η(t) < ∞,
and hence ψ t,x (t,x)∈[0,T ]×R d is uniformly integrable. The same arguments using (3.12) show that ψt,x W(t, T t (1) , x, ∆W t (1) ) (t,x)∈[0,T ]×R d is uniformly integrable.

Further discussions

Representation of the PDE system Let us consider a PDE system (v j ) j=1,••• ,n , where for each j, v j :

[0, T ] × R d → R satisfies ∂ t v j + µ j • Dv j + 1 2 a j : D 2 v j + f (•, v 1 , • • • , v n , Dv 1 , • • • , Dv n ) = 0,
for some diffusion coefficient function (µ j , a j ) : [0, T ] × R d -→ R d × S d , and some polynomial function

f : [0, T ] × R d × R n × (R d ) n -→ R.
Our methodology immediately applies to this context, and provides a stochastic representation for the solution of the above PDE system, by means of a regime-changed branching diffusions: at every branching time, the independent offspring particles perform subsequently different branching diffusion regime.

Representation in view of unbiased simulation With the same idea of proof, we can also obtain an alternative representation result, with a frozen coefficient SDE in place of SDE (2.4). When the coefficient function a ≡ a 0 for some constant a 0 ∈ S d , this has significant application in terms of Monte Carlo approximation, as it leads to a representation random variable which can be simulated exactly, while the branching diffusion process X k • in (2.4) needs a time discretization technique and hence creates some discretization error in the simulation. Let us present this alternative representation formula in the case of constant diffusion coefficient case, i.e. a ≡ a 0 = σ 0 σ 0 for some non-degenerate constant matrix σ 0 ∈ M d .

Let L := L∪{∂}, where ∂ represents an artificial index; p = (p ) ∈L be a probability mass function and ( I m,n ) m,n≥1 be a sequence of i.i.d. random variables of distribution p, and independent of the sequences of i.i.d Brownian motion ( W m,n ) m,n≥1 and i.i.d positive random variable ( T m,n ) m,n≥1 of density function ρ. Then following exactly the same procedure in Section 2.2, we can construct another age-dependent branching process, denoted by ( T k ) k∈ K T with branching type I k := I n,πn(k) . Here, when

I k = ( ˆ 0 , • • • , ˆ m ) ∈ L, it produces | ˆ |
offspring particles, marked by i = 0, • • • , m exactly as in Step 3 in the construction of age-dependent process K T in Section 2.2; when I k = ∂, it produces only one offspring particle, marked by m + 1. Then for every k ∈ K T , we equipped it with an independent Brownian motion W k

• as in (3.2). Next, let us define X

(1) 0 = x 0 , and subsequently for every k ∈ K T ,

X k T k := X k T k- + µ( T k-, X k T k- )∆ T k + σ 0 ∆ W k ∆ T k , with X k T k- := X k- T k- . (4.1)
For this case, the automatic differentiation functions take a particularly simple formula, which is compatible with the purpose of the unbiased simulation algorithm. Let us introduce

W k := 1 {θ k =0} + b θ k ( T k-, X k T k- ) • (σ 0 ) -1 ∆ W k ∆ T k ∆ T k 1 {θ k ∈{1,••• ,m}} (4.2)
+ µ( T k-, X k T k- ) -µ( T (k-)-, X k- T (k-)- ) • (σ 0 ) -1 ∆ W k ∆ T k ∆ T k 1 {θ k =∂} .
Finally, setting c ∂ ≡ 1, and replacing (X k • , W k , p I k ) in the definition of ψ and ψ (in and below (3.4)) by ( X k

• , W k , p I k ), we obtain

ψ := k∈ K T g( X k T ) -g( X k T k-)1 θ k =0 F (∆T k ) W k k∈ K T \ K T c I k ( T k , X k T k ) p I k W k ρ(∆ T k ) , (4.3)
and similarly ˜ ψ. Next, given a constant vector µ 0 ∈ R d , we keep the same branching Brownian motion (

W k • ) k∈ K T
, and then introduce another diffusion process X µ 0 ,k

• by

X µ 0 ,k T (1)
:= x 0 + µ 0 ∆ T (1) + σ 0 ∆ W (1)
∆T (1) , and the subsequent process

X µ 0 ,k T k for k ∈ K T \ {(1)
} by the same induction relation as in (4.1). We then introduce W µ 0 k as in (4.2) by replacing X k by X µ 0 ,k , and replacing

µ( T (k-)-, X k- T (k-)- ) by µ 0 when k = (1). Replacing ( X k , W k ) by ( X µ 0 ,k , W µ 0 k ) in (4.
3), it defines a new random variable ψ µ 0 . Finally, by changing the initial condition (0, x 0 ) and time interval [0, T ] to (t, x) and [t, T ], one obtains

W t,k , T t k , ψ t,x , ˜ ψ t,x
, ψ µ 0 ,t,x etc.

Proposition 4.1. Suppose that Assumptions 3.1 holds true, and the semilinear PDE (2.2) has uniqueness for bounded viscosity solution. Suppose in addition that for every (t, x) ∈ [0, T ] × R d , and µ 0 lies in a neighborhood of µ(t, x), one has ψ µ 0 ,t,x and ψ µ 0 ,t,x ∆ W t,(1)

∆ T t (1)
/∆ T t,( 1) is integrable, and the family of random variables

( ψ t,x ) (t,x)∈[0,T ]×R d and ˜ ψ t,x ∆ W t,(1) ∆ T t (
:= E[ ψ µ 0 ,t,x ] and v(t, x) := E ˜ ψ t,x ∆ W t,(1) ∆ T t (1)
/∆ T t, [START_REF] Athreya | Branching processes[END_REF] .

Notice that v is uniformly bounded and continuous. Recall that W is a standard ddimensional Brownian motion independent of the branching diffusion process, we also introduce

X t,x s := x + µ 0 (s -t) + σ 0 (W s -W t ), s ∈ [t, T ],
where µ 0 ∈ R d is a constant vector in a neighborhood of µ(t, x). Then one obtains as in (3.5) that ũ(µ 0 , t, x) = E 1

F ( T t (1) ) g X t,x T t (1)
1 { T t (1) =T } + 1 { I (1) =∂} ρ( T t (1) )p ∂ (µ -µ 0 ) • v T t (1) , X t,x T t (1) + 1 { I (1) =∂} ρ( T t (1) ) c I (1) p I (1) û I (1),0 m i=1 (b • v) I (1),i T t (1) , X t,x T t (1) 1 { T t (1) <T } = E g X t,x T + T t (µ -µ 0 ) • v + f (•, û, v) s, X t,x s ds .
By standard argument, (t, x) → ũ(µ 0 , t, x) is a viscosity solution of

-∂ t u + µ 0 • Du + 1 2 a 0 : D 2 u + (µ -µ 0 ) • v + ∈L c û 0 m i=1 (b i • v) i = 0,
with terminal condition g. Since û and v are bounded continuous, the above PDE has uniqueness for bounded viscosity solution, which induces that ũ(µ 0 , t, x) is independent of µ 0 and û(t, x) = ũ(µ 0 , t, x) for µ 0 in a neighborhood of µ(t, x).

(ii) We can then compute the derivative D x ũ(µ 0 , t, x) and then set µ 0 := µ(t, x), it follows that

Dû(t, x) = D x ũ(µ(t, x), t, x) = v(t, x),
which is also bounded continuous. This implies that û(t, x) is a viscosity solution of (2.2), and we hence conclude the proof by uniqueness of the viscosity solution of (2.2).

The integrability and square integrability of ψ can be analyzed in exactly the same way as in Theorem 3.12. We just notice that the above defined random variable ψ can be simulated exactly from a sequence of Gaussian random variable, discrete distributed random variables I m,n and r.v. T m,n of distribution density function ρ. It is then in particular interesting to serve as a Monte-Carlo estimator for u(0, x 0 ).

On the representation of fully nonlinear PDEs Formally, one can also obtain a representation result for fully nonlinear PDE, using the same automatic differentiation functions of order 2. However, this raises a serious integrability problem which can not be solved by conditions as in Assumption 3.10. To illustrate the main difficulty, let us consider the following PDE in the one-dimensional case d = 1:

u(T, x) = g(x), ∂ t u + 1 2 D 2 u + f 0 (D 2 u) = 0, on [0, T ] × R d , (4.4)
where f 0 (γ) = c 0 γ for some constant c 0 > 1 2 . Notice that there is only one term in function f 0 , then a natural guess for the representation is to consider a branching Brownian motion with exactly one offspring particle at every arrival time. This can be seen as a Brownian motion W equipped with a sequence of random time mark (T i ) i=1,••• ,N T , where

T i := T ∧ i j=1 τ j,1 , N T := inf i : T i ≥ T .
Notice that for any t > 0 and bounded measurable function φ : R → R, one has

∂ 2 xx E φ(x + W t ) = E φ(x + W t ) W 2 t -t t 2 .
Then arguing as in Theorem 3.5, we may expect that u(0, x 0 ) = E[ ψ , with

ψ := g(x + W T ) 1 F (T -T N T -1 ) N T -1 i=1 c 0 (W T i+1 -W T i ) 2 -(T i+1 -T i ) (T i+1 -T i ) 2 ρ(T i -T i-1 ) ,
provided that ψ is integrable. However, the integrability of ψ could fail in general. For simplicity, let g ≡ 1, and notice that F ≤ 1. Then by taking conditional expectation, one has, for some constant C > 0 and

c 1 := E c 0 (W 2 1 -1) , that E | ψ| ≥ E N T -1 i=1 c 1 (T i+1 -T i )ρ(T i -T i-1 ) ≥ E c 1 ρ(T 1 ) c 1 (T 2 -T 1 )ρ(T 2 -T 1 ) c 1 T 3 -T 2 1 {T 1 ≤T /2,T 2 -T 1 <T /2,T 3 -T 2 ≥T } ≥ CE c 1 (T 2 -T 1 )ρ(T 2 -T 1 ) 1 {T 2 -T 1 <T /2} = C T /2 0 1 t dt = ∞.
Of course, for linear PDEs as in (4.4), one can simulate a Brownian motion with volatility coefficient 1 + 2c 0 whenever 1 + 2c 0 > 0 to obtain the solution. But it is not the case for general fully nonlinear PDEs.

On the representation results by BSDE Another probabilistic representation

of semilinear parabolic PDE is the Backward Stochastic Differential Equation (BSDE) proposed by Pardoux and Peng [START_REF] Pardoux | Adapted solutions of backward stochastic differential equations[END_REF]. Namely, given a classical solution u of semilinear PDE (2.2), we define

(Y t , Z t ) := u(t, X 0,x 0 t ), σDu(t, X 0,x 0 ) . Then (Y, Z) provides a solution to BSDE Y t = g(X 0,x 0 T ) + T t f s, X 0,x 0 s , Y s , σ -1 (s, X 0,x 0 s )Z s ds -Z s dW s , t ∈ [0, T ], P-a.s.
Based on the discretization technique on the BSDE, one can then obtain a probabilistic numerical solution for semilinear parabolic PDEs, see e.g. Bouchard and Touzi [START_REF] Bouchard | Discrete-time approximation and Monte-Carlo simulation of backward stochastic differential equations[END_REF], and Zhang [START_REF] Zhang | A numerical scheme for backward stochastic differential equations[END_REF], etc. Generally speaking, these numerical schemes for BSDE need a (time-consuming) simulation-regression technique to compute the conditional expectation appearing in the schemes.

Our representation result induces a pure Monte Carlo simulation algorithm, which avoids the regression procedure in the numerical schemes of BSDEs. Nevertheless, our numerical method provides only the solution of PDE at time 0, and it needs some restrictive conditions on the coefficient functions f such as Assumption 3.10 to obtain a finite variance estimator. We will provide more numerical examples as well as some variance reduction techniques in Section 5 below.

A Monte Carlo algorithm 5.1 The implementation of the numerical algorithm

The above representation result in Theorem 3.5 induces a Monte Carlo algorithm to compute the solution of PDE (2.2), by simulating the random variable ψ or ψ. We provide here some discussion on the implementation of the numerical algorithm.

The choice of density function ρ As discussed in Remark 3.11, to ensure Assumption 3.10, a necessary condition is to choose ρ(t) ≥ Ct -1/2 . A natural candidate as distribution, which is also easy to be simulated, is the gamma distribution Γ(κ, θ), with κ ≤ 1 2 , whose density function is given by

ρ 0 (t) = 1 Γ(κ)θ κ t κ-1 exp(-t/θ)1 {t>0} , (5.1)
where Γ(κ) := ∞ 0 s κ-1 e -s ds. In particular, one has

F k := ∞ ∆T k ρ 0 (t)dt = 1 - γ(κ, ∆T k /θ) Γ(κ)
, where γ(κ, t) := t 0 s κ-1 e -s ds.

Complexity The dimension d of the problem, the choice of (p ) ∈L and ρ will of course influence the complexity of algorithm. First, the complexity is proportional to the number of particles in the branching process, i.e. #K T , and for each particle, the complexity of simulation and calculation is of order Cd

m(t) = ∞ k=0 n k 0 F * ,k (t), where F * ,k (t) := P τ 1,1 + • • • + τ 1,k < t . (5.2)
(ii) Let ρ be given by (5.1), then F * ,k (t) = 1 Γ(kκ) γ(kκ, x/θ) and hence

m(t) = ∞ k=0 γ(kκ, t/θ)n k 0 Γ(kκ) .
Proof. (i) Using Lemma 4.4.3 of Athreya and Ney [START_REF] Athreya | Branching processes[END_REF], one has that satisfies the equation m(t) = 1 + n 0 t 0 m(t -s)ρ(s)ds, whose solution is given explicitly by (5.2). Further, when ρ is the density function of Gamma distribution, the function F * ,k (t) can be computed explicitly.

A high dimensional numerical example

We first focus on a simple numerical example in high dimension. Let (µ, σ) ≡ (0, σ 0 ) for some constant matrix σ 0 = 1

√ d I d , and f (t, x, y, z) = k(t, x) + cy(b • z), where b := 1 d (1 + 1 d , 1 + 2 d , • • • , 2) and k(t, x) := cos(x 1 + • • • + x d ) α + σ 2 2 + c sin(x 1 + • • • + x d ) 3d + 1 2d e α(T -t) e α(T -t) .
With terminal condition g(x) = cos(x

1 + • • • + x d ), the explicit solution of semilinear PDE (2.2) is given by u(t, x) = cos(x 1 + • • • + x d )e α(T -t) .
In our numerical experiment, we set α = 0.2, c = 0.15, T = 1, and x 0 = 0.51I d , where 1I d stands for the unit vector in R d for d = 5, 10 and 20. We would like to emphasize that, to the best of our knowledge, no alternative methods are available for solving such a high-dimensional semilinear PDE. In Table For numerical implementations, we use gamma distribution (5.1), with κ = 0.5 and θ = 2.5. On each test performed, a computation is achieved with n particles. An estimation E with n particles is then calculated. The standard deviation of E is estimated with 1000 runs of n particles and its log-plot is reported below on the different figures for different values of n. We also show on some figures the convergence of the solution obtained on the average of the 1000 runs.

On Figures 1, 2, 3, we illustrate that the Monte Carlo method converges easily to our analytic solution. Computational costs are estimated on one core of a Laptop core I7 processor 2.2 GHz and are reported in Table 2 for a number of simulations equal to 96000 permitting to get a solution with an error less than 0.1%.

X 0,x t = (1 -e -t )1I + e -t x + σ 1 -e -2t 2 Z, Z ∼ N (0, I d ). (6.1)
The Malliavin weight used in the algorithm can be computed explicitly and is given by Z/ σ e 2t -1 2 . We will compare numerical results from four different schemes.

• (scheme a) using the representation (3.4) with the explicit solution (6.1) of the SDE (3.1).

• (scheme b) using the representation (4.3) with freezing coefficient techniques.

• (scheme c) using the representation (4.3), enhanced by the resampling scheme (see Appendix for more details).

• (scheme d) using the representation (3.4), enhanced by the resampling scheme.

The density function ρ is that of the gamma law with parameters κ and θ. If not indicated, the parameters of the law are set to κ = 0.5 and θ = 2.5 and the probability p l are chosen equal. On each test, a calculation is achieved with n particles (starting with n = 1562 for scheme a and with n = 100000 for schemes b and c). The n particles are shared on 96 processors and each processor i calculates an estimation E i of the solution with n 96 particles. Then an estimation E with n particles is achieved with E = 1 96 i E i . When importance sampling is used, in order to avoid communications that breaks parallelism, it is used on each processor so with n 96 particles on each processor. The standard deviation of E is estimated with 1000 runs of n particles and its log is reported on the different figures below for different values of n. We expect that by quadrupling the values of n, the standard deviation std divides by a factor 2 and the plot (log(n), log(std)) should be linear with a slope equal to -1 2 . The theoretical rate of convergence is also plotted on each figure (as in our previous example, the solutions are obtained on the average of the 1000 runs).

Some examples in one space dimension

• For d = 1, we take f (t, x, y, z) := 0.2y 2 + 0.3y 3 . Results on Figure 4 show that the method converges. Scheme a is far more effective than scheme b and that the importance sampling of scheme c is effective. The log of the standard deviation decreases for all schemes linearly with the log of the particle number as predicted by the theory. Note that the computational cost for 1000 runs with 25000 particles on one core is equal to 490 seconds for scheme a, 200 seconds with scheme b and 260 seconds with scheme c. • As a second example in d = 1, we take a Burgers type nonlinearity f (t, x, y, z) = 0.15yz. Results on Figure 5 show that all the schemes converge to our numerical finite difference solution. Note that the computational cost for 1000 runs with 25000 particles on one core is roughly equal to 200 seconds for scheme a, 100 seconds for scheme b, 300 seconds for scheme c. • As a third example in d = 1, we keep the same nonlinearity with f (t, x, y, z) = 0.3yz. We expect that the variance of the results will be higher than in the previous case. This is observed in Figure 6. Scheme a still converges. Scheme b converges slowly and Importance Sampling of scheme c permits to get faster convergence and to recover the good rate in the log of the standard deviation decay. The computational times are the same as in our previous test.

• As a fourth example in d = 1, we take a nonlinearity with f (t, x, y, z) = 0.08z 2 . Results are shown in Figure 7. The importance sampling of scheme c is required • As a last example in d = 1, we keep the same type of nonlinearity f (t, x, y, z) = 0.2z 2 . Schemes b and c don't converge anymore. We only test scheme a using different values for the parameters κ and θ (see Figure 8). The change in θ does not seem to change convergence properties. The change in κ (from 0.5 to 0.4) does not seem to modify our results. However, some tests, not reported here, show that the variance can increase a lot using κ around 0.25. Then, as the average jump size is proportional to θ, it is more efficient to take some quite high values for θ in order to reduce the computational time. For the same reason, it is optimal to choose a κ equal to 0.5. In Table 3, we report the computational time, associated to different choices of (κ, θ), as a multiplicative factor of the computational effort with benchmark parameters κ = 0.5, θ = 2.5.

We notice that for all the parameters, the decay in the variance is far from the expected theoretical one (see Figure 8). We then use our benchmark parameters and compare the results obtained using scheme a and scheme d (importance sampling is used here).

Results are reported on Figure 9. They illustrate that the importance sampling method allows to improve the convergence rate.

Some examples in two space dimensions

Although the efficiency of our algorithm was illustrated on our previous experiments, this Monte-Carlo method cannot compete a PDE deterministic methods in d = 1. In this section, we focus on d = 2, where advantages of PDE implementation remain but are not so obvious.

• For the first example in d = 2, we take f (t, x, y, z) := 0.15y1I.z. On Figure 10, we give the results obtained using our three schemes showing that Importance Sampling is needed. Note that the computation cost for 1000 runs with 25000 particles on one core is roughly equal to 230 seconds for scheme a, 90 seconds for scheme b, 580 seconds for scheme c. • For the third example, we test the influence of the coefficients on Scheme a for a non linearity f (t, x, y, z) := K(z.1I) 2 with K = 0.05, K = 0.1, K = 0.2. Using Scheme b and c, we cannot get proper convergence due to high variances observed. On Figure 12, we give the convergence obtained with the different K values and on Figure 13 the standard deviation associated. As the coefficients grow, the variance of the results gets higher preventing the method from converge when K = 0.2.

• At last we test the influence of the function g. The representation of the solution involves the product of g functions so we expect that the variance of the result is highly sensitive to the scaling of this function. Here we choose to keep f (t, x, r, p) := 0.05(Du.1I) 2 and take different values for the g function. On figure 14 we take g(x) = 2( 1

An example in three space dimensions

We take f (t, x, y, z) := 0.15(z.1I) 2 . Results are given on Figure 16, still showing that importance sampling is necessary while using discretization of the scheme and that the exact scheme has a lower variance.

A Resampling scheme for branching processes

Notice that our estimator (3.4) and (4.3) are provided as a product of some random variables. Then similar to Doumbia, Oudjane and Warin [START_REF] Doumbia | Computing expectations for general SDE with pure Monte Carlo methods[END_REF], one can use the resampling scheme (or interacting particle systems), see Del Moral [START_REF] Moral | Feynman-Kac formulae, Genealogical and interacting particle systems with applications, Probability and its Applications[END_REF]. Intuitively, this scheme replaces the expectation of a product by a product of expectations, which potentially stabilizes the Monte-Carlo estimator.

Let us first introduce the Markov chain (X n ) n≥1 , taking values in ∪ p≥1 ([0, T ] 2 × R 2d × {0, • • • , m} × L) p such that X 1 = (0, 0, x 0 , x 0 , 0, I 0 ) with I 0 = (1, 0, • • • 0) ∈ N m+1 and for any n ≥ 1, one defines

X n+1 := X n , (T k -, T k , X k T k -, X k T k , θ k , I k ) k∈∪ n p=1 K p T
.

Notice that this Markov chain has an absorbing state since for any ω ∈ Ω there is a generation n(ω) for which all branches have died (either having no offspring before reaching T or having reached T ) implying K n+1 T (ω) = ∅ and consequently X n+1 (ω) = X n (ω). Then (X n ) n≥0 is a Markov chain. We next introduce

G n (X n ) := k∈K n T g(X k T ) -g(X k T k-)1 {θ k =0} F (∆T k ) W k k∈(K n T \K n T ) c I k (T k , X k T k ) p I k W k ρ(∆T k ) , (A.1) so that ψ = ∞ n=1 G n (X n ).
Notice that the above representation consists of a product from contributions from each generation n ≥ 1. Since the number of generation prior to the maturity T is finite a.s., the last product only involves finite number of terms, a.s. We also observe that except for the trivial case of constant function g, E[|G n (X n )|] = 0. By iteration, it is easy to see that

E P 0 [ψ] = ∞ n=1 E P n-1 G n (X n ) E P∞ ∞ n=1 sgn(G n (X n )) ,
where given P 0 , one defines P n by dPn dP n-1 := |Gn(Xn)|

E P n-1 [|Gn(Xn)|]
, for n ≥ 1. The particle algorithm consists in simulating the dynamics of an interacting particle system of size N , (ξ Finally u(0, x 0 ) is approximated as a product of empirical averages:

∞ n=1 1 N N i=1 |G n (ξ i,N n )| 1 N N i=1 ∞ n=1 sgn G n (ξ i,N n ) . (A.2)
Notice again that, for every simulation (ξ i,N ), the maturity T is attained for some finite generation, then the above product ∞ n=1 can be restricted to the a finite product

2

  The marked branching diffusion 2.1 Semilinear PDE with polynomial nonlinearity Let d ≥ 1, M d denotes the set of all d × d matrices, and (µ, σ) : [0, T ] × R d → R d × M d the coefficient functions. For a function u : [0, T ] × R d → R, we denote by Du and D 2 u the gradient and the Hessian of the function u(t, x) w.r.t. variable x. Let m ≥ 0 be a positive integer, we consider a subset L ⊂ N m+1 , and a sequence of functions (c ) ∈L and (

Assumption 3 . 1 .

 31 (i) The probability mass function (p ) ∈L satisfies p > 0 for all ∈ L, and ∈L | | p < ∞. The density function ρ : R + → R + is continuous and strictly positive on [0, T ], and such that

5 .

 5 Define |ϕ| ∞ := sup x∈R d |ϕ(x)| for any bounded function ϕ : R d → R, and |φ| ∞ := sup d i=1 |φ i | ∞ for any bounded vector function φ : R d → R d .

2 .

 2 Let us denote n 0 := ∈L p | | and m(t) := E #K t . Proposition 5.1. (i) The function m(t) is given by

Figure 1 :

 1 Figure 1: Estimation and standard deviation observed in d = 5 depending on the log of the number of particles used.

Figure 2 :

 2 Figure 2: Estimation and standard deviation observed in d = 10 depending on the log of the number of particles used.

Figure 3 :Table 2 :

 32 Figure 3: Estimation and standard deviation in d = 20 depending on the log of the number of particles used.

Figure 4 :

 4 Figure 4: Estimation and standard deviation obtained in d = 1 for f (t, x, y, z) := 0.2y 2 + 0.3y 3

Figure 5 :

 5 Figure 5: Estimation and standard deviation obtained in d = 1 for f (t, x, y, z) = 0.15yz.

Figure 6 :

 6 Figure 6: Estimation and standard deviation obtained in d = 1 for f (t, x, y, z) = 0.3yz.

Figure 7 :

 7 Figure 7: Estimation and standard deviation obtained in d = 1 for f (t, x, y, z) = 0.08z 2 .

Figure 8 :Figure 9 :

 89 Figure 8: Estimation and standard deviation obtained in d = 1 for f (t, x, y, z) = 0.2z 2 . Different values for κ and θ are used.

Figure 10 :

 10 Figure 10: Estimation and standard deviation obtained in d = 2 for f (t, x, y, z) := 0.15y1I.z.

Figure 11 :

 11 Figure 11: Estimation and standard deviation obtained in d = 2 for f (t, x, y, z) := 0.04(z.1I) 2 .

  i) -1) + and give the convergence of schemes a and

Figure 12 :

 12 Figure 12: Convergence of scheme a for different K values.

Figure 13 :

 13 Figure 13: Standard deviation of the scheme a for different K values.

Figure 14 :

 14 Figure 14: Estimation and standard deviation observed in dimension 2 for case 4, g(x) = 2( 1 d

Figure 15 :

 15 Figure 15: Estimation and standard deviation observed in dimension 2 for case 4, g(x) = 3( 1 d

Figure 16 :

 16 Figure 16: Estimation and standard deviation observed in dimension 3 for case 1.

;

  1,N p , • • • ξ N,N p ), on ∪ p≥1 ([0, T ] 2 × R 2d × {0, • • • , m} × L) p , from step n = 1 to n = ∞and then to approximate each expectation E P n-1 G n (X n ) by the empirical mean value of the simulation. The algorithm can be given as an iteration of the following two steps, initiated by n = 1,Selection step Given N copies of simulation (ξ i,N n ) i=1,••• ,N of X n, one draws randomly and independently N particles among the current particle system with a probability|Gp(ξ i,N p )| N j=1 |Gp(ξ i,N p |)Evolution step Each new selected particle evolves randomly and independently according to the transition of the Markov chain (X n ) between n and n + 1.

  n N := inf{n ≥ 1 | ξ i,N n has reached T for all i = 1, • • • N }.

  a particle of generation n, with arrival time T k . When T k < T , we let I k = I n,πn(k) , where π n is an injection from N n to N, and at time T k , it branches into |I k | offspring particles, which constitute n+1-the generation, and are indexed by (k 1

  Then the age-dependent branching process is well defined on [0, T ], i.e. the number of particles in K t is finite a.s. for all t ∈ [0, T ].

	, 1, 1, 2)
	Proposition 2.3. Assume that
	Proof. See e.g. Theorem 1 of Athreya and Ney [1, Chapter IV.1], or Harris [12, pp.
	138-139].

∈L | |p < ∞.

  Proposition 3.7. Let Assumption 3.6 hold true, then Assumption 3.2 holds true with the choice of automatic differentiation function W defined by W t, s, x, (W r -W t ) r∈[t,s] := 1 s -t Remark 3.8. When µ ≡ 0 and σ(t, x) ≡ σ 0 for some non-degenerate constant matrix σ 0 ∈ M d , one then has Y t ≡ I d and so that

	s	σ -1 (r, X	t,x r )Y r dW r .	(3.10)
	t

T ], P-a.s.,

(3.9)

where I d denotes the d × d identity matrix, and σ i (t, x) ∈ R d denotes the i-th column of matrix σ(t, x). Then one has the following result (see e.g. Exercise 2.3.5 of Nualart

[19, p.p. 125]

, or Proposition 3.2. of

[START_REF] Fournié | Applications of Malliavin Calculus to Monte Carlo Methods in Finance[END_REF]

).

  1 2 , 1). Theorem 3.12. Consider the automatic differentiation function (3.10), and suppose that Assumptions 3.1, 3.6 and 3.10 hold true. (i) Then Assumptions 3.2 holds, and ψ t,x , ψt,x W(t, T t (1) , x, ∆W t (1)

  Then the derivative Dû exists, û and Dû are both continuous; and moreover, u is the unique bounded viscosity solution of semilinear PDE (2.2).

	Sketch of proof. (i) First, by the uniform integrability condition, û is bounded
	continuous. Let us introduce
	ũ(µ 0 , t, x)

1) /∆ T t,(1) (t,x)∈[0,T ]×R d are uniformly integrable with uniformly bounded expectation, define û(t, x) := E[ψ t,x ].

Table 1 :

 1 1, we report the analytic solution of the semilinear PDE and that of the corresponding linear PDE by setting c = 0. The different results indicate that the nonlinearity term has an impact. Analytical solution for the linear PDE (i.e., c = 0) versus analytical solution for the semilinear PDE in d = 5, 10 and 20.

	Dimension	5	10	20
	Linear Solution	-1.0436 0.3106 -0.9661
	Non linear solution -0.97851 0.34646 -1.0248

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).