L Abbas-Turki
email: lokmane.abbasturki@sorbonne-universite.fr.

B Diallo

G Pagès

Conditional Monte Carlo Learning for Diffusions I: main methodology and application to backward stochastic differential equations

We present a new algorithm based on a One-layered Nested Monte Carlo (1NMC) to simulate functionals U of a Markov process X. The main originality of the proposed methodology comes from the fact that it provides a recipe to simulate U t≥s conditionally on X s . Because of the nested structure that allows a Taylor-like expansion, it is possible to use a very reduced basis for the regression. Although this methodology can be adapted for a large number of situations, we only apply it here for the simulation of Backward Stochastic Differential Equations (BSDEs). The generality and the stability of this algorithm, even in high dimension, make its strength. It is heavier than a straight Monte Carlo (MC) but it is far more accurate to simulate quantities that are almost impossible to simulate with MC. The parallel suitability of 1NMC makes it feasible in a reasonable computing time. This paper explains the main version of this algorithm and provides first results of error estimates. We also give various numerical examples with a dimension equal to 100 that are executed from few seconds to few minutes on one Graphics Processing Unit (GPU).

Introduction

Numerous contributions in numerical methods based on Monte Carlo reached recently their limits in dealing with the curse of dimensionality [START_REF] Bellman | Dynamic programming[END_REF]. This paper is a natural progress of an increasing interest in Nested Monte Carlo (NMC) started in [START_REF] Gordy | Nested Simulation in Portfolio Risk Measurement[END_REF][START_REF] Lee | Monte Carlo Computation of Conditional Expectation Quantiles[END_REF][START_REF] Lee | Computing the distribution function of a conditional expectation via Monte Carlo: Discrete conditioning spaces[END_REF] and used with regression in [START_REF] Abbas-Turki | XVA principles, nested Monte Carlo strategies, and GPU optimizations[END_REF][START_REF] Broadie | Risk estimation via regression[END_REF] and with Multilevel method in [START_REF] Giorgi | Weak Error for Nested Multilevel Monte Carlo[END_REF]. In this contribution, Onelayered Nested Monte Carlo (1NMC) is set for only two layers of Monte Carlo that offer the possibility of having two layers of approximation. In contrast to previous works, our methodology is based on a judicious combination between 1NMC and various localized regressions.

Although we are not the first to propose a learning procedure for BSDEs or option pricing/hedging [START_REF] Buehler | Deep hedging[END_REF][START_REF] Han | Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations[END_REF][START_REF] Goudenège | Variance Reduction Applied to Machine Learning for Pricing Bermudan/American Options in High Dimension[END_REF], we are the first to do it using NMC and local regressions instead of a neural network. Thus, we replace standard learning using a global minimization to train a neural network by local regressions that provide very accurate results due to 1NMC. Combining regression with 1NMC makes the local projections accurate, even on a reduced basis, as long as the projected values are sufficiently smooth. Of course, local neural networks can replace local regressions but, as will be shown, the latter are already quite sufficient for the considered examples and have very efficient implementation using batch parallel processing [START_REF] Abbas-Turki | Resolving small random symmetric linear systems on graphics processing units[END_REF].

The presented procedure is called Conditional Monte Carlo Learning for Diffusions (CMCLD) as the regression-based analytical representations are established on the top of a set of outer trajectories of diffusions. If these outer trajectories are erased the whole learning procedure is lost. The conditional learning is the main drawback of this methodology but we are working on making the learning unconditional in future works. The goal of this paper is to lay the foundation of the simplest version of this methodology making possible further extensions of it in Conditional Monte Carlo Learning for Diffusions part II (CMCLDII) [START_REF] Abbas-Turki | Conditional Monte Carlo Learning for Diffusions II: extended methodology and application to risk measures and early stopping problems[END_REF].

On the top of what is presented in this part I, part II [START_REF] Abbas-Turki | Conditional Monte Carlo Learning for Diffusions II: extended methodology and application to risk measures and early stopping problems[END_REF] will be dedicated to further error controls and slight modifications that make possible accurate simulation of tail values as well as optimal stopping problems in high/very high dimensions. Part I and II have a common purpose of explaining in depth CMCLD and its benefit either when the dimension is high/very high or when the problem can be formulated as a composition of various functionals of a Markov process.

On a filtered probability space (Ω, F, (F t ) 0≤t≤T , P), the generic example is driven by an F t -Markov process (X t ) t∈[0,T ] taking its values on R d 1 . Given the fine time discretization S = {t 0 , ..., t 2 L } = 0, T /2 L , ..., T , let U s be a functional of X defined for s ∈ S by

(f ) U s = u s (X s ) = E s   2 L k=s2 L /T f (t k , X t k , X t k+1 )   = E s≤t k ≤T f (t k , X t k , X t k+1 ) F s ,
where E s (•) = E • F s , the expectation is always considered under P, each deterministic function f (t k , •, •) is B(R d 1 ) ⊗ B(R d 1 )-measurable and assumed to satisfy the square integrability condition E(f 2 (t k , X t k , X t k+1 )) < +∞ with convention f (t 2 L , X t 2 L , X t 2 L +1 ) = f (t 2 L , X t 2 L ). From the Markov assumption, for each s ∈ S, the B(R d 1 )-measurable function u s is deterministic. The simulation of U is generic to all BSDE examples presented in this paper. As nested simulations involve heavy notations, it is easier to present the whole algorithm for the simulation of U then apply it on specific examples.

When previous contributions target estimations of U t k for k = 0, ..., 2 L knowing some realization of {X t j } 0≤j≤k (m 0 = 1, ..., M 0 ), our purpose is to simulate approximations {U m 0 ,m 1 t k ,s } s≥t k+1 , with (m 0 = 1, ..., M 0 ) and (m 1 = 1, ..., M 1 ), of {U s } s≥t k+1 conditionally on the realization {X m 0 t j } 0≤j≤k . This task requires the simulation of a first layer (X m 0 ) m 0 =1,...,M 0 of trajectories that are kept in the machine's random-access memory, then a second (X m 0 ,m 1 ) m 1 =1,...,M 1 unstored layer of trajectories, on the top of the first layer indexed by m 0 , only used to learn how should we perform approximations U m 0 ,m 1 called inner layer simulation of the process U .

Although more complex than a regression/MC method (cf. [START_REF] Lemor | Rate of convergence of an empirical regression method for solving generalized backward stochastic differential equations[END_REF][START_REF] Gobet | Adaptive importance sampling in leastsquares Monte Carlo algorithms for backward stochastic differential equations[END_REF]), getting the inner layer simulation provides much more information on the process that has to be simulated.

In particular, it is possible to set the outer layer simulation U m 0 t k of U t k to be equal to 1

M 1 M 1 m 1 =1 f (t k , X m 0
t k , X m 0 ,m 1 t k ,t k+1 ) + U m 0 ,m 1 t k ,t k+1 . Knowing the inner layer simulation U m 0 ,m 1 , we can compute quantiles on U or, even more remarkable, can simulate another process U that satisfies equation ( f ) (Replace f by f in equation (f )) with an f that can be a function of U like for instance f (t k , x, y) = f (t k , U t k (x), U t k+1 (y)). Consequently, when sufficient assumptions are satisfied, we can learn how to compute functionals of functionals of X with the same 1NMC. This latter fact makes possible the simulation of Valuation Adjustments [START_REF] Abbas-Turki | XVA principles, nested Monte Carlo strategies, and GPU optimizations[END_REF] as long as one can write them as a composition of functionals then start simulating by the innermost functional till the most outer composition.

The outer simulation U m 0 is then defined using the inner one U m 0 ,m 1 , the latter obtained from regressing on the inner trajectories X m 0 ,m 1 . When U m 0 t k is defined for all k = 0, ..., 2 L , U m 0 ,m 1 has not to be defined on the whole fine discretization S = {t 0 , ..., t 2 L } = 0, T /2 L , ..., T . We introduce then a coarser discretization set S that controls the depth of our learning procedure by adding/reducing regression steps needed for sufficiently accurate values of U m 0 . It is possible to improve the depth in a parareal fashion [START_REF] Lions | A "parareal" in time discretization of PDE's[END_REF] which increases further the parallel scalability of the algorithm. Controlling the depth of our learning procedure makes possible the balance between complexity and accuracy. This compromise can be even formulated iteratively as explained in [START_REF] Abbas-Turki | Conditional Monte Carlo Learning for Diffusions II: extended methodology and application to risk measures and early stopping problems[END_REF].

Having inner and outer simulations allows also to control the bias throughout the whole algorithm. Indeed, for any s ∈ [0, T [ and any s ∈]s, T ] one can estimate both parts of equality

E(U s ) = E U s + s t l+1 >s f (t l , X t l , X t l+1 ) (1.1) using M 0 m 0 =1 U m 0 s M 0 , 1 M 0 M 0 m 0 =1 U m 0 s + s t l+1 >s f (t l , X m 0 t l , X m 0 t l+1
) and their difference that represents the estimated value of the average bias. Controlling the difference is numerically possible using a particular choice of the time step s, associated to any s ∈ S, at which we define the terminal condition involved in the inner simulation. To simulate U s , we can subsequently set a terminal condition u m 0 ,S s,s (•) to be only equal to f (T, •) when s = T . The choice of s and of the terminal condition u m 0 ,S s,s (•) are discussed in the example of Section 2.2, are presented in definitions 3.1 and 3.2 with an upper bound error established in CMCLDII [START_REF] Abbas-Turki | Conditional Monte Carlo Learning for Diffusions II: extended methodology and application to risk measures and early stopping problems[END_REF]. Combining the tower property in equality (1.1) with the nested structure of our algorithm is very beneficial to stop the bias propagation. Further benefits of this combination on the variance of regressed values is presented in [START_REF] Abbas-Turki | Conditional Monte Carlo Learning for Diffusions II: extended methodology and application to risk measures and early stopping problems[END_REF].

Starting with a simple application of the generic example, Section 2 introduces the method as well as notations. In Section 3, we define both inner U m 0 ,m 1 and outer U m 0 simulations of the process U given in (f ) then we adapt it to BSDEs with a Markov forward process. Using similar arguments to the one presented in [START_REF] Broadie | Risk estimation via regression[END_REF], Section 4 explains the impact of inner regressions on outer simulations. Section 5 shows the robustness of our methodology on highly dimensional problems beyond what is known to be possible in previous contributions.

2 Notations and the methodology applied on a simple example

We introduce in Section 2.1 some needed notations and set the stage for Section 2.2 that explains CMCLD starting from a simple application.

2.1 Notations: 1NMC, stabilized regressions, coarse discretization and δ time operator

1NMC stands for one-layered nested MC and thus involves two layers of Monte Carlo simulation of a Markov process X. Using a sufficiently fine discretization S = {t 0 , ..., t

2 L } = {0, ∆ t , 2∆ t , ..., T } with ∆ t = T /2 L , one simulates M 0 realizations (X m 0 t k ) m 0 =1,...,M 0 k=1,...,2 L
of the Markov process X starting at a deterministic point X 0 = x 0 ∈ R d 1 with the following induction

X m 0 t k = E t k-1 (X m 0 t k-1 , ξ m 0 t k ), when k ≥ 1 and X m 0 t 0 = x 0 , (2.1)
where (ξ

m 0 t k ) m 0 =1,...,M 0 k=1,...,2 L
are independent realizations of an R d 2 random vector ξ and the functions (

E t k ) k=0,...,2 L -1 : R d 1 +d 2 → R d 1 are Borel-measurable. We use X m 0 ,1 t k , ..., X m 0 ,d 1 t k to denote the d 1 components of the vector X m 0 t k . The sample (X m 0 t k ) m 0 =1,...,M 0 k=1,...,2 L
stays on the machine memory and is supposed to approximate accurately (X t ) t∈[0,T ] in a sense explained in Section 4.

For a decreasing sequence (s j ) j=0,...,2 L that takes its values in the time discretization set S, an extra simulation conditional to the starting X m 0 s j is needed for the learning procedure. Introducing independent realizations (ξ m 0 ,m 1 t j ,t k ) (m 0 ,m 1 )∈{1,...,M 0 }×{1,...,M 1 +M 1 } k∈{j,...,2 L },j∈{1,...,2 L } of the random vector ξ that are also independent from (ξ m 0 t k ) m 0 =1,...,M 0 k=1,...,2 L , we set for

t k-1 ≥ s j X m 0 ,m 1 s j ,t k = E t k-1 (X m 0 ,m 1 s j ,t k-1 , ξ m 0 ,m 1 s j ,t k ) and X m 0 ,m 1 s j ,s j m 1 =1,...,M 1 +M 1 = X m 0 s j . (2.2)
We use

X m 0 ,m 1 ,1 s j ,t k , ..., X m 0 ,m 1 ,d 1 s j ,t k
to denote the d 1 components of the vector X m 0 ,m 1 s j ,t k . For s j ≤ s l ≤ s k , we also introduce the notation X m 0 ,m 1 s j ,s l :s k and ξ m 0 ,m 1 s j ,s l :s k for respectively (X m 0 ,m 1 s j ,s l , X m 0 ,m 1 s j ,s l +∆t , ..., X m 0 ,m 1 s j ,s k -∆t , X m 0 ,m 1 s j ,s k ) and (ξ m 0 ,m 1 s j ,s l , ξ m 0 ,m 1 s j ,s l +∆t , ..., ξ m 0 ,m 1 s j ,s k -∆ , ξ m 0 ,m 1 s j ,s k ). In (2.2), we simulate M 1 + M 1 conditional realizations of X in order to keep those indexed from m 1 = M 1 + 1 to m 1 = M 1 + M 1 for the approximation of regression matrices. Consequently, we make explicit the independence between trajectories used for the estimation of regression matrices and those used in the backward induction. This makes also explicit the difference between the number of trajectories needed for regression matrices approximation and the number of those used for the backward induction. To reduce the complexity of the algorithm and memory occupation, trajectories used for regression matrices can be simulated offline then erased from the memory. Given m 0 , if the inner trajectories {X m 0 ,m 1 } m 1 =1,...,M 1 are needed α times in the backward induction, we simulate α independent copies and use each copy once. This reduces further memory occupation as well as any superfluous dependence structure.

For each ordered couple (j < k) of indices that take their values in {1, ..., 2 L }, we introduce the stabilized regression basis

T m 0 t j ,t k ,M 1 : R d 1 x → t Γ m 0 t j ,t k ,M 1 x -X m 0 t k ∈ R d 1 (d 1 ≤ d 1 ), (2.3)
that performs a linear combination of the components of x -X m 0 t k using Γ m 0 t j ,t k ,M 1 that contains some eigenvectors from Γ m 0 t j ,t k ,M 1 obtained with the eigenvalue decomposition

Γ m 0 t j ,t k ,M 1 Λ m 0 t j ,t k ,M 1 t Γ m 0 t j ,t k ,M 1 (2.4) of the regression matrix 1 M 1 M 1 +M 1 m 1 =M 1 +1 X m 0 ,m 1 t j ,t k -X m 0 t k t X m 0 ,m 1 t j ,t k -X m 0 t k (2.5)
where t is the transpose operator. Once factorization (2.4)=(2.5) is performed, we obtain the diagonal matrix Λ m 0

t j ,t k ,M 1 = diag λ m 0 ,l t j ,t k ,M 1 l=1,...,d 1
of decreasing positive eigenvalues. Then, we define Λ m 0

t j ,t k ,M 1 = diag λ m 0 ,l t j ,t k ,M 1 l=1,...,d 1
as the truncation of Λ m 0 t j ,t k ,M 1 with d 1 defined by

d 1 = min    k ∈ {1, .., d 1 }, k l=1 λ m 0 ,l t j ,t k ,M 1 ≥ p% d 1 l=1 λ m 0 ,l t j ,t k ,M 1    , (2.6)
where p% ∈ [95%, 100%], d 1 keeps only eigenvalues that make the regression problem wellconditioned i.e. The ratio

λ m 0 ,l t j ,t k ,M 1 λ m 0 ,1 t j ,t k ,M 1 l=1,...,d 1
has to be bigger than 10 -6 in single precision or bigger than 10 -15 in double precision floating representation [START_REF] Press | Numerical Recipes in C++: The Art of Scientific Computing[END_REF]. In addition to ensuring a well-conditioned regression problem, equality (2.6) also performs a principal component analysis [START_REF] Press | Numerical Recipes in C++: The Art of Scientific Computing[END_REF]. At the same time that we set the components of Λ The regression with respect to t Γ m 0

t j ,t k ,M 1 X m 0 ,m 1 t j ,t k -X m 0 t k ∈ R d 1 with d 1 ≤ d 1 , instead of X m 0 ,m 1 t j ,t k -X m 0 t k ∈ R d 1
, involves the inversion of the diagonal matrix Λ m 0 t j ,t k ,M 1 which replaces the whole regression matrix (2.5). Since Λ is bounded below away from zero, its inverse is bounded and the same for the regression procedure. This stabilizes the computation of the regression estimator whose expression is detailed in Section 3. Besides, since we have a large number of regression matrices, we can batch compute these inversions like explained in [START_REF] Abbas-Turki | Resolving small random symmetric linear systems on graphics processing units[END_REF]. The latter reference presents a recipe to resolve efficiently and accurately large number of small symmetric linear systems on GPUs. In particular, the authors of [START_REF] Abbas-Turki | Resolving small random symmetric linear systems on graphics processing units[END_REF] present a batch parallelization strategy of Cuppen's divide & conquer algorithm (cf. [START_REF] Cuppen | A divide and conquer method for the symmetric tridiagonal eigenproblem[END_REF]) suited for CMCLD. As shown in the numerical examples of Section 5, performing a regression with re-

spect to t Γ m 0 t j ,t k ,M 1 X m 0 ,m 1 t j ,t k -X m 0 t k
for any couple of fine increments (t j < t k ) is mostly both unacessary and computationally heavy. Thus, only regressions with respect to

t Γ m 0 s j ,s k ,M 1 X m 0 ,m 1 s j ,s k -X m 0 s k
are considered, with (s j , s k ) being a couple that take their values in S × S ∩ [s j , s j ] . The set of the possible values S is a subset of the fine time discretization S. The couple (s j , s j ) ∈ (S∩]s j , T ]) 2 sets the minimal and the maximal time increment of the performed regressions conditionally on the realization X m 0 s j of the Markov process X. Given a coarse discretization set S, definitions 3.1 and 3.2 express the value of (s, s) s∈S .

The choice of S and of S s j → (s j , s j ) ∈ (S∩]s j , T ]) 2 have to be specific to each application. In some examples like the one presented in Section 5.6, S = S or S has to be fine enough because fine discretization is needed for the corresponding BSDE. Also for some situations as the Hamilton-Jacobi-Bellman simulation presented in Section 5.4, one has to cut the bias propagation and set s j < T for some values s j ∈ S. Consequently, S and S s → (s, s) ∈ (S∩]s, T ]) 2 are not predefined but rather add a degree of freedom to increase accuracy and stop the bias propagation. In Conditional Monte Carlo Learning for Diffusions part II, for i = 0, ..., L -L , we define (s j ) j=0,...,2 L to be iteratively equal to (T -s i j ) j=0,...,2 L starting with a homogeneously distributed sequence where each term is repeated 2 L-L times as follows

s 0 j = j2 L 2 L T 2 L where • is the ceiling function. (2.7)
Denoting S i the set of values taken by (T -s i j ) j=0,...,2 L , for example S 0 = {0, T /2 L , ..., (2 L -1)T /2 L , T }, it is possible to set a refinement strategy.

Given that (s j ) j∈{0,...,2 L } is a decreasing, and not strictly decreasing, sequence of coarse increments, we need to define on S a new operator δ S that associates to each s ∈ S the next increment in S. For a fixed index j ∈ {1, ...,

2 L }, we define δ S s j (•) on (s k ) k≤j , taken its values in S ∩ [s j , s j ], by δ S s j (s k ) = min (s j , min{s ∈ S; s k < s ≤ s j }) (2.8) with min(∅) = ∞.
When there is no confusion on the chosen set S, we use δ s j notation instead of δ S s j . When s k < s j , we use δ S notation instead of δ S s j . When there is no confusion on the chosen set S and s k < s j , we simplify both indices and use δ instead of δ S s j . This time operator will be largely used and for a given set S it has the following properties Pr1. s j = δ s j (s j ) = δ(s j ).

Pr2. As long as max(s

j 1 , s j 2 ) ≤ s k < min(s j 1 , s j 2 ), δ s j 1 (s k ) = δ s j 2 (s k ) = δ(s k ). Pr3. The nth composition of δ s j denoted δ n s j (•) is equal to s j when n ≥ |S∩]s j , s j ]| where | • | denotes the cardinal.

Illustration of coarse and fine approximations for conditional expectations

We consider here the following process

U t = E f (X T ) X t ,
with a deterministic function f . Thus, we assume that there is no path dependence through the sum on the realizations of X as done in (f ). In this path-independent situation for the fixed time set (2.7), it is clear that one can simulate U using 1NMC without any need of regression and thus without using our method. However, we choose to illustrate our method on this simple case and we will see at the end of this section what could be the benefits. For j, j ∈ {1, ..., 2 L } and time steps s j , s j ∈ S -{0} with s j < s j and for a fixed outer trajectory (X m 0 t k ) k=0,...,2 L , let us assume that we want to simulate one realization of U s j and of U s j . A straight way to do it is to draw the inner trajectories {X m 0 ,m 1 s j ,T } m 1 =1,...,M 1 and {X m 0 ,m 1 s j ,T } m 1 =1,...,M 1 , as in Figure 1, then average with respect to m 1 on the realizations f (X m 0 ,m 1 s j ,T ) and f (X m 0 ,m 1 s j ,T

) respectively. If s j and s j are close to each other in some sense 1 , our method makes possible to simulate U t k for any t k ∈ [s j , s j ) using

U m 0 t k = u m 0 ,S t k ,T = 1 M 1 M 1 m 1 =1 u m 0 ,S s j ,s j (X m 0 ,m 1 t k ,s j ).
Thus, u m 0 ,S s j ,s j (x) and

X m 0 ,m 1 t k ,s j replace respectively f (x) and X m 0 ,m 1 t k ,T
involved in a standard nested simulation that would approximate the m 0 realization of U t k by

M 1 m 1 =1 f (X m 0 ,m 1 t k ,T ) M 1 .
The introduced

• u m 0 ,S s,s (x) (0 ≤ s < s < s), called coarse approximation, is the regression representa- tion of E f (X T ) X s = x computed using the inner trajectories {X m 0 ,m 1 s,s } 1≤m 1 ≤M 1 . • u m 0 ,S s,s (0 ≤ s < T ), called fine approximation, is the Monte Carlo representa- tion of E f (X T ) X s computed

by averaging on the regression representations

{u m 0 ,S s,s (X m 0 ,m 1 s,s )} 1≤m 1 ≤M 1 with s = δ(s). The index s in u m 0 ,S s,s
provides the depth of the regression learning conditionally on the realization X m 0 s . Below, we explain how u m 0 ,S s,s (x) should be computed and how s should be set. First of all, since s j and s j are assumed to be "close enough" one can consider the same learning depth s j = s j = t 2 L = T and define

(u s j T ) & (u s j T ) u m 0 ,S s j ,T (x) = u m 0 ,S s j ,T (x) = f (x) and if s j = T then ( u s j ) U m 0 s j = u m 0 ,S s j ,T = 1 M 1 M 1 m 1 =1 u m 0 ,S s j ,T (X m 0 ,m 1 s j ,T
).

As will be presented in (3.1), we define

(u s j s j ) u m 0 ,S s j ,s j (x) = u m 0 ,S s j ,T + t T m 0 s j ,s j ,M 1 (x)A m 0 ,S s j ,s j ,
where the adaptation of (3.3) makes

(A T s j s j ) A m 0 ,S s j ,s j = Λ m 0 s j ,s j ,M 1 -1 M 1 M 1 m 1 =1 T m 0 s j ,s j ,M 1 (X m 0 ,m 1 s j ,s j ) u m 0 ,S s j ,T X m 0 ,m 1 s j ,T -u m 0 ,S s j ,T
.

1 Not necessary an Euclidean distance.

' T '' ''' 0 '''' Figure 1:
Given the realization of one outer trajectory (bold), we simulate inner trajectories to approximate U s j , U s j , U s j , U s j and U s j .

If we add a third increment s j (cf. Figure 1) such that s j , s j and s j are close enough, for any

t k ∈ [s j , s j ) one can set U m 0 t k = u m 0 ,S t k ,T
. The latter equality requires the definition of u m 0 ,S s j ,s j which can be obtained from (u s j s j ) (replace s j by s j and s j by s j in (u s j s j )) involving u m 0 ,S s j ,T and A m 0 ,S s j ,s j that can be computed using (A s j s j s j ). The calculations in (A s j s j s j ) use u m 0 ,S s j ,T and u m 0 ,S s j ,s j whose expression depends on u m 0 ,S s j ,T and A m 0 ,S s j ,s j . Finally, A m 0 ,S s j ,s j is the regression vector of u m 0 ,S s j ,T around u m 0 ,S s j ,T . Subsequently, the computations of u m 0 ,S s j ,T , u m 0 ,S s j ,T and u m 0 ,S s j ,T involve the dependence structure given in (2.9).

u m 0 ,S s j ,T → u m 0 ,S s j ,s j → u m 0 ,S s j ,T → u m 0 ,S s j ,s j → u m 0 ,S s j ,T → u m 0 ,S s j ,T = f ↑ ↑ A m 0 ,S s j ,s j A m 0 ,S s j ,s j → u m 0 ,S s j ,T = f u m 0 ,S s j ,s j → u m 0 ,S s j ,T → u m 0 ,S s j ,T = f ↑ A m 0 ,S s j ,s j → u m 0 ,S s j ,T = f (2.9)
By adding other increments s j and s j (cf. Figure 1), it can happen that s j , s j and s j can no longer be considered close enough to each other. In this situation, a regression on a linear basis around X m 0 s j would not be considered sufficient for inner trajectories that start at X m 0 s j or X m 0 s j . To deal with this situation, one should introduce (s j , s j ) and (s j , s j ) that set the beginning and the ending of the family of successive regressions. For instance if s j = s j and s j = s j , one starts the backward induction associated to the increment s j by the final condition u m 0 ,S s j ,s j (x) = u m 0 ,S s j ,s j (x) instead of u m 0 ,S s j ,T (x) = f (x) and the dependence tree (2.9) becomes

u m 0 ,S s j ,s j → u m 0 ,S s j ,s j → u m 0 ,S s j ,s j → u m 0 ,S s j ,s j → u m 0 ,S s j ,s j → u m 0 ,S s j ,s j → u m 0 ,S s j ,T ...f ↑ ↑ ↑ A m 0 ,S s j ,s j A m 0 ,S s j ,s j A m 0 ,S s j ,s j ...f u m 0 ,S s j ,s j → u m 0 ,S s j ,s j → u m 0 ,S s j ,s j → u m 0 ,S s j ,T ...f ↑ ↑ A m 0 ,S s j ,s j A m 0 ,S s j ,s j ...f u m 0 ,S s j ,s j = u m 0 ,S s j ,s j (2.10)
In Figure 1, we also set s j = s j as well as s j = s j and the tree (2.10) can be further changed to include the dependency structure induced by s j . Indeed, we urge the reader to check that (2.10) can be as easily completed as done for (2.9) to include the dependency structure induced by s j . Even with the simple example presented in this subsection, one can show the benefit of this method. Indeed, in addition to a fine simulation of U using u, this method defines a set of functions u that can be considered as coarse conditional approximation of U . These conditional approximations can be used as forward components of another functional. For instance, given the example presented above and illustrated in Figure 1, the simulation of an m 0 realization of V s j = E (U s j -U s j ) + X s j can be done with

V m 0 s j = 1 M 1 M 1 m 1 =1 u m 0 ,S s j ,s j (X m 0 ,m 1 s j ,s j ) -u m 0 ,S s j ,s j (X m 0 ,m 1 s j ,s j ) + .
These functions u can be also used for risk measures. For example, the conditional value at risk VaR α% U s j -U s j X s j of level α% can be computed after sorting the values u m 0 ,S s j ,s j (X m 0 ,m 1 s j ,s j ) -u m 0 ,S s j ,s j (X m 0 ,m 1 s j ,s j )

1≤m 1 ≤M 1 .
Remark 2.1. Referring to Figure 1, for any g, when E g(U s j ) X s j , E g(U s j ) X s j and E g(U s j ) X s j are well defined their simulation can be directly performed using u m 0 ,S s j ,s j , u m 0 ,S s j ,s j or u m 0 ,S s j ,s j . This is not the case for E g(U s j ) X s j since u m 0 ,S s j ,s j were not computed because s j = s j < s j .

The other benefit of our methodology is the possibility to have a parareal alike implementation [START_REF] Lions | A "parareal" in time discretization of PDE's[END_REF] and thus make the algorithm parallel in time in addition to have it parallel in paths. Indeed, refering to Figure 1, if we associate the final conditions u m 0 ,S s j ,s j and u m 0 ,S s j ,s j respectively to each subinterval [s j , s j ) and [s j , s j ), we can perform concurrent calculations on these intervals.

General methodology applied to BSDEs with a Markov forward process

Based on what was presented previously, we detail in Section 3.1 the simulation of approximations of U defined by (f ). Section 3.2 illustrates the adaptation of this new method to BSDEs.

Fine and coarse approximations

Considering the discretization sequence (s j ) j=0,...,2 L that takes its values in the set S ⊂ S, we use a learning procedure to associate to each scenario m 0 and to each discretization set S a couple of function families ( u m 0 ,S , u m 0 ,S ). Now, for given indices k, j ∈ {1, ..., 2 L } (k < j) that satisfy s j < s k < s j , for x ∈ R d 1 and s ∈ {s k , s k + ∆ t , ..., δ(s k ) -∆ t }, we define two approximation levels: A coarse approximation around X m 0 s k conditionally on X m 0 s j defined by

u m 0 ,S s j ,s k (x) = u m 0 ,S s k ,s k + t T m 0 s j ,s k ,M 1 (x)B m 0 ,S s j ,s k , (3.1)
and a fine approximation at X m 0 s defined by

u m 0 ,S s,s k = 1 M 1 M 1 m 1 =1   u m 0 ,S s k ,δ(s k ) (X m 0 ,m 1 s,δ(s k ) ) + δ(s k ) t l+1 >s f (t l , X m 0 ,m 1 s,t l , X m 0 ,m 1 s,t l+1 )   . (3.2)
The value of B m 0 ,S s j ,s k is given by

B m 0 ,S s j ,s k = ( Λ m 0 s j ,s k ,M 1 ) -1 1 M 1 M 1 m 1 =1 B m 0 ,S,δs j (s k ) s j ,s k ,M 1 (X m 0 ,m 1 s j ,s k :δs j (s k ) ) (3.3)
where

X m 0 ,m 1 s j ,s k :δs j (s k ) = X m 0 ,m 1 s j ,s k , X m 0 ,m 1 s j ,s k +∆t , ..., X m 0 ,m 1 s j ,δs j (s k )-∆t , X m 0 ,m 1 s j ,δs j (s k )
and B m 0 ,S,δs j (s k )

s j ,s k ,M 1 : Ω × R d 1 (δs j (s k )-s k )/∆t (ω, x 1 , ..., x (δs j (s k )-s k )/∆t ) → Ω × R d 1 is F δs j (s k ) ⊗ B(R d 1 (δs j (s k )-s k )/∆t
)measurable and defined by

B m 0 ,S,δs j (s k ) s j ,s k ,M 1 (x) = T m 0 s j ,s k ,M 1 (x 1 )        u m 0 ,S s j ,δs j (s k ) x δs j (s k )-s k ∆ t -u m 0 ,S s k ,δs k (s k ) + δs j (s k )-s k ∆ t -1 l=1 f (t ks k +l , x l , x l+1 )        B m 0 ,S,δs j (s k ) s j ,s k (x) (3.4)
where

k s k = s k /∆ t -1 and x = (x 1 , ..., x (δs j (s k )-s k )/∆t ).
T involved in (3.1) was already defined in (2.3). Regarding the regression vector B m 0 ,S s j ,s k , its value can be seen as an estimation of the vector a ∈ R d 1 that minimizes the quadratic error given by

E B m 0 ,S,δs j (s k ) s j ,s k (X m 0 ,m 1 s j ,s k :δs j (s k ) ) -t aT m 0 s j ,s k ,M 1 (X m 0 ,m 1 s j ,s k ) 2 . (3.5)
To complete this inductive interconnected backward definition of u and u, we set the final coarse approximation to

u m 0 ,S s j ,s j (x) = f (T, x) if s j = T, u m 0 ,S s j ,s j (x) = u m 0 ,S δ(s j ),s j (x) if s j < T, (3.6)
where the values s j , s j ∈ (S∩]s j , T ]) 2 with s j < s j , expressed in Definition 3.1, delimit the learning depth of the performed regressions conditionally on the realization X m 0 s j of the Markov process X. In numerical examples of Section 5, we see that s j and s j are needed when the bias becomes important because either T is sufficiently big or the variance produced by X is large enough. Otherwise, (3.6) can be replaced by u m 0 ,S s j ,T (x) = f (t 2 L , x) = f (T, x). In CMCLDII [START_REF] Abbas-Turki | Conditional Monte Carlo Learning for Diffusions II: extended methodology and application to risk measures and early stopping problems[END_REF], we suggest a method of iterative actualization of the couple {(s, s)} s∈S when we refine the discretization set S.

According to equations (3.1), (3.2), (3.3), (3.4) and (3.6), the functions u and u are defined backwardly. When u is a straight Monte Carlo involving u, the latter is defined using a regression around a point at which we expressed u. Consequently, u can be seen as a conditional first order Taylor expansion around the first layer of trajectories (X m 0 t k ) m 0 =1,...,M 0 k=1,...,2 L . The term of order zero in this expansion is played by u, where the term t T m 0 s j ,s k ,M 1 (x)B m 0 ,S s j ,s k , deduced from the minimization of (3.5), plays the order one.

Remark 3.1.

1. Since we do not want to increase further the algorithm complexity by considering higher order terms, the definition of u involves only linear terms for regressions around X m 0 s k .

2. When the dimension d 1 is not too high, it is possible to regress the residual of the first regression on higher order terms. These successive regressions do not increase drastically the complexity when compared to the standard procedure. Nevertheless, as it separates regression with respect to first order terms and regression with respect to higher order terms, it loses orthogonality between first and higher order terms.

3. In case X is a martingale, the linearity simplifies further computations since, for instance, (3.2) can be replaced by

u m 0 ,S s,s k = u m 0 ,S s k ,δ(s k ) (X m 0 s ) + 1 M 1 M 1 m 1 =1 δ(s k ) t l+1 >s f (t l , X m 0 ,m 1 s,t l , X m 0 ,m 1 s,t l+1 ).
Definition 3.1. Given a discretization set S ⊂ S

• For any s ∈ S, s is set to be equal to δ(s) and s is set backwardly to be the largest discretization time u ∈ S∩]s, δ(s)] that satisfies

1 M 0 M 0 m 0 =1   u m 0 ,S s,u -u m 0 ,S δ(s),δ(s) - δ(s) t l+1 >s f (t l , X m 0 t l , X m 0 t l+1 )   < S s (3.7)
where { S s } s∈S is a family of positive bias tuning parameters.

• For k, j ∈ {1, ..., 2 L } (k < j) that satisfy s j < s k ≤ s j < t 2 L = T , the simulation U m 0 ,m 1 s j ,s k of U around X m 0 s k conditionally on X m 0 s j is set to be equal to u m 0 ,S s j ,s k (X m 0 ,m 1 s j ,s k ) where u is given in (3.1) and (3.6).

• For k ∈ {1, ..., 2 L } and s ∈ {s k , s k + ∆ t , ..., δ(s k ) -∆ t } -{0}, the simulation U m 0 s of U at X m 0 s is set to be equal to u m 0 ,S s,s k with u expressed in (3.2).

• The average U lear 0 of learned values on U 0 is equal to

U lear 0 = 1 M 0 M 0 m 0 =1 u m 0 ,S 0,0 (3.8)
and the simulated value U sim 0 of U 0 is equal to

U sim 0 = 1 M 0 M 0 m 0 =1   u m 0 ,S δ(0),δ(0) + δ(0) t l+1 >0 f (t l , X m 0 t l , X m 0 t l+1 )   (3.9) with u expressed in (3.2).
The definition of s is straightforward since we take s = δ(s). To set the value of s, we need to know δ(s) whose value involves the one of δ 2 (s) and so on till T . In [START_REF] Abbas-Turki | Conditional Monte Carlo Learning for Diffusions II: extended methodology and application to risk measures and early stopping problems[END_REF], Section 3.1, we give an example of how to update the values of (s, s) s∈S when the coarse discretization S gets finer.

U m 0 ,m 1 can be seen as the inner or second layer approximation of U and U m 0 can be seen as the outer or first layer approximation of U . When U m 0 ,m 1 is only defined on S, it is remarkable to see that U m 0 is defined on the whole fine discretization set S. At t = 0, inner and outer trajectories coincide that yield U lear 0 and U sim 0 as two possible outer first layer approximations of U 0 . For any S, it is natural to have S s proportional to the value of the estimation 1

M 0 M 0 m 0 =1 u m 0 ,S
s,s . Used to control the bias, the choice of S s has also to take into account the confidence interval of the estimator of the left side of inequality (3.7). Based on equality (1.1), (3.7) is quite sufficient in the considered examples to have almost unbiased estimates. CMCLDII [START_REF] Abbas-Turki | Conditional Monte Carlo Learning for Diffusions II: extended methodology and application to risk measures and early stopping problems[END_REF] introduces a more stringent local bias control that can be used when (3.7) is not sufficient.

BSDEs with a Markov forward process

In Section 2.2, we saw the implementation of CMCLD on a simple problem and we showed its benefits when one has to simulate functionals of functionals of a Markov process. BSDEs with a Markov forward process are specific functionals of functionals of a Markov process. After [START_REF] Pardoux | Adapted solution of a backward stochastic differential equation[END_REF], BSDEs became very widely studied, especially in the quantitative finance community starting with [START_REF] El Karoui | Backward stochastic differential equations in Finance[END_REF]. Here we adapt CMCLD to the One step forward Dynamic Programming (ODP ) scheme for discrete BSDEs

(ODP ) Y T = ζ and for k < 2 L Y t k = E t k [Y t k+1 + ∆ t f (t k , Y t k+1 , Z t k )], Z t k = E t k [Y t k+1 (W t k+1 -W t k )/∆ t ].
(ODP ) was studied for instance in [START_REF] Gobet | Error expansion for the discretization of backward stochastic differential equations[END_REF][START_REF] Lemor | Rate of convergence of an empirical regression method for solving generalized backward stochastic differential equations[END_REF]. We consider ζ = f (T, X T ) to be some square integrable random variable that depends on X T . Given a discretization sequence (s j ) j=0,...,2 L ∈ S and referring to (2.1) and (2.2), the simulation of X involves the increments of an

R d 2 -Brownian motion W with ξ m 0 t k = W m 0 t k -W m 0 t k-1 and ξ m 0 ,m 1 s j ,t k = W m 0 ,m 1 s j ,t k - W m 0 ,m 1 s j ,t k-1 where W 1 , ..., W M 0 are independent realizations of W with W m 0 ,m 1 s j ,t k = W m 0 ,m 1 s j ,t k-1 + ∆W m 0 ,m 1 s j ,t k and W m 0 ,m 1 s j ,s j m 1 =1,...,M 1 +M 1 = W m 0 s j , (∆W m 0 ,m 1 s j ,t k ) (m 0 ,m 1 )∈{1,...,M 0 }×{1,...,M 1 +M 1 } k∈{j,...,2 L },j∈{1,...,2 L } are independent Brownian motion increments in- dependent from W 1 , ..., W M 0 with E([∆W m 0 ,m 1 s j ,t k ] 2 ) = ∆ t .
As pointed out in Section 2.1, if an inner trajectory {X m 0 ,m 1 } is needed several times in the backward induction, we simulate independent copies of it and thus independent copies of ξ m 0 ,m 1 and use each copy once.

For given indices k, j ∈ {1, ..., 2 L } with k < j that satisfy s j < s k < s j and using δ s j (s k ) defined in (2.8), we also set ∆W m 0 ,m 1 s j ,s k ,δs j (s k ) = W m 0 ,m 1 s j ,δs j (s k ) -W m 0 ,m 1 s j ,s k . For each k, the Borel B(R) ⊗ B(R d 2 )-measurable driver f (t k , •, •) is assumed to satisfy Lipschitz condition of Section 4.

Given the discretization set S, one can define two coarse approximations around X m 0 s k conditionally on X m 0 s j given by

y m 0 ,S s j ,s k (x) = y m 0 ,S s k ,s k + t T m 0 s j ,s k ,M 1 (x)C m 0 ,S s j ,s k , (3.10)
t z m 0 ,S s j ,s k (x) = t z m 0 ,S s k ,s k + t T m 0 s j ,s k ,M 1 (x)D m 0 ,S s j ,s k , (3.11)
as well as two fine approximations at X m 0 s , for s ∈ {s k , s k + ∆ t , ..., δ s j (s k ) -∆ t } with ∆ s = δ s j (s k ) -s and ∆ s k = δ s j (s k ) -s k , given by

y m 0 ,S s,s k = 1 M 1 M 1 m 1 =1   ∆ s f (s k , y m 0 ,S s k ,δs j (s k ) (X m 0 ,m 1 s,δs j (s k ) ), z m 0 ,S s,s k ) + y m 0 ,S s k ,δs j (s k ) (X m 0 ,m 1 s,δs j (s k ) )   , (3.12
)

z m 0 ,S s,s k = 1 M 1 ∆ s M 1 m 1 =1 y m 0 ,S s k ,δs j (s k ) (X m 0 ,m 1 s,δs j (s k ) ) W m 0 ,m 1 s,δs j (s k ) -W m 0 s (3.13)
and we set the final coarse approximation to

y m 0 ,S s j ,s j =    f (t 2 L , X m 0 ,m 1 s j ,t 2 L ) if s j = t 2 L ,
y m 0 ,S s j ,s j (X m 0 ,m 1 s j ,s j ) = y m 0 ,S δs j (s j ),s j (X m 0 ,m 1 s j ,s j ) if

s j < t 2 L , (3.14)
where the values s j , s j ∈ (S∩]s j , T ]) 2 with s j < s j , expressed in Definition 3.2, delimit the learning depth of the performed regressions conditionally on the realization X m 0 s j of the Markov process X.

Since T was already expressed in (2.3), to complete this inductive interconnected definition of (y, y, z, z), we set the vector C m 0 ,S s j ,s k and the matrix D m 0 ,S s j ,s k to be equal to

C m 0 ,S s j ,s k = ( Λ m 0 s j ,s k ,M 1 ) -1 M 1 M 1 m 1 =1 Y m 0 ,S,δs j (s k ) s j ,s k ,M 1 (X m 0 ,m 1 s j ,s k , X m 0 ,m 1 s j ,δs j (s k ) ), (3.15)
D m 0 ,S s j ,s k = ( Λ m 0 s j ,s k ,M 1 ) -1 M 1 M 1 m 1 =1 Z m 0 ,S,δs j (s k ) s j ,s k ,M 1
(X m 0 ,m 1 s j ,s k , ∆W m 0 ,m 1 s j ,s k ,δs j (s k ) , X m 0 ,m 1 s j ,δs j (s k ) ), (3.16) with Y m 0 ,S,δs j (s k )

s j ,s k ,M 1 (x , x) = T m 0 s j ,s k ,M 1 (x )Y m 0 ,S,δs j (s k ) s j ,s k (x , x) is F δs j (s k ) ⊗ B(R 2d 1 )-measurable and Z m 0 ,S,δs j (s k ) s j ,s k ,M 1 (x , w, x) = T m 0 s j ,s k ,M 1 (x ) t Z m 0 ,S,δs j (s k ) s j ,s k (w, x) is a vector function measurable with respect to F δs j (s k ) ⊗ B(R 2d 1 +d 2 ), where Y m 0 ,S,δs j (s k ) s j ,s k (x , x) =   ∆ s k f (s k , y m 0 ,S s j ,δs j (s k ) (x), z m 0 ,S s j ,s k (x )) + y m 0 ,S s j ,δs j (s k ) (x) -y m 0 ,S s k ,s k   (3.17) and Z m 0 ,S,δs j (s k ) s j ,s k (w, x) = y m 0 ,S s j ,δs j (s k ) (x) w ∆ s k -z m 0 ,S s k ,s k . (3.18)
From equations above, one can associate quadratic minimization problems to C m 0 ,S s j ,s k and to D m 0 ,S s j ,s k , as done for B m 0 ,S s j ,s k in (3.5). In the same fashion as in Definition 3.1, we define the double layer approximations (Y m 0 , Z m 0 ) and (Y m 0 ,m 1 , Z m 0 ,m 1 ) of functionals Y and Z. Definition 3.2. Given a discretization set S ⊂ S • For any s ∈ S, s is set to be equal to δ(s) and s is set backwardly to be the largest discretization time u ∈ S∩]s, δ(s)] that satisfies

1 M 0 M 0 m 0 =1 y m 0 ,S s,u -y m 0 ,S δ(s),δ(s) -(δ(s) -s)f (s, y m 0 ,S δ(s),δ(s) , z m 0 ,S s,s ) < S s (3.19)
where { S s } s∈S is a family of positive bias tuning parameters.

• For k, j ∈ {1, ..., 2 L } with k < j that satisfy s j < s k ≤ s j < t 2 L = T , the simulation

Y m 0 ,m 1 s j ,s k and Z m 0 ,m 1 s j ,s k
of Y and Z respectively around X m 0 s k conditionally on X m 0 s j are set to be equal to y m 0 ,S s j ,s k (X m 0 ,m 1 s j ,s k ) and z m 0 ,S s j ,s k (X m 0 ,m 1 s j ,s k ) where y and z are given in (3.10), (3.11) and (3.14). with y and z expressed in (3.12) and (3.13).

• For k ∈ {1, ..., 2 L } and s ∈ {s k , s k + ∆ t , ..., δ s k (s k ) -∆ t } -{0},
• The average Y lear 0 and Z lear 0 of learned values on Y 0 and Z 0 are respectively equal to Y m 0 ,m 1 and Z m 0 ,m 1 are the inner or second layer approximation of Y and Z respectively. Y m 0 and Z m 0 are the outer or first layer approximation of Y and Z respectively. Similar to what was presented for the generic example in Section 3.1, as inner and outer trajectories coincide at t = 0, we get two possible outer first layer approximations of Y 0 and Z 0 . It is also remarkable to have Y m 0 and Z m 0 defined on the whole fine discretization set S despite Y m 0 ,m 1 and Z m 0 ,m 1 are only defined on S. With BSDEs, it is possible to replace (3.19), using rather an M DP scheme (cf. [START_REF] Gobet | Linear regression MDP scheme for discrete backward stochastic differential equations under general conditions[END_REF]), by

Y lear 0 = 1 M 0 M 0 m 0 =1 y m 0 ,S 0,0 , Z lear 0 = 1 M 0 M 0 m 0 =1 z m 0 ,S 0,0 (3.20
1 M 0 M 0 m 0 =1 y m 0 ,S s,u -y m 0 ,S δ(r),δ(r) - r θ∈S,θ=s (δ(θ) -θ)f (s, y m 0 ,S δ(θ),δ(θ) , z m 0 ,S θ,θ ) < S s , ( 3

Regression-based NMC and bias propagation

Before presenting the main elements, we point out that we have intentionally considered only discrete functionals of a Markov process. The approximation due to discretization of the continuous version of BSDEs is not studied and we refer to [START_REF] Bouchard | Discrete time approximation and Monte Carlo simulation of backward stochastic differential equations[END_REF][START_REF] Gobet | Error expansion for the discretization of backward stochastic differential equations[END_REF][START_REF] Lemor | Rate of convergence of an empirical regression method for solving generalized backward stochastic differential equations[END_REF] among others that quantify well the resulting error. Moreover, we also consider the discretized version of the Markov process introduced in (2.1) and (2.2) where

E t k (x, ξ) = x + ∆ t b(t k , x) + σ(t k , x)ξ (4.1)
with the usual (cf. [START_REF] Pagès | Numerical probability: an introduction with applications to finance[END_REF]) Lipschitz continuity condition on the coefficients b(t, x) and σ(t, x) uniformly with respect to t ∈ [0, T ]. Similar to what was considered in Section 3.2, the noise ξ is given by increments of a vector of independent Brownian motions i.e.

ξ m 0 t k = W m 0 t k -W m 0 t k-1 and ξ m 0 ,m 1 s j ,t k = W m 0 ,m 1 s j ,t k -W m 0 ,m 1 s j ,t k-1 . (4.
1) can be read as an Euler scheme of a stochastic differential equation that admits a strong solution. In this paper, when the discretization is needed, we assume that L is sufficiently large to neglect the discretization error of the forward process X.

Given two arbitrary square integrable random variables χ 1 and χ 2 , consider {χ 3

m 1 } M 1 m 1 =1
to be the empirical regression of χ 1 with respect to χ 2 , the authors of [START_REF] Broadie | Risk estimation via regression[END_REF] established an upper bound error of the regression-based NMC estimator 1

M 1 M 1 m 1 =1 φ(χ 3 m 1 ) of E(φ(E(χ 1 |χ 2 )))
once we know the representation error κ = E(χ 1 |χ 2 ) -t RB(χ 2 ) induced by the projection of E(χ 1 |χ 2 ) on the basis B(χ 2 ). The fine approximations u and y presented earlier were computed by averaging on the empirical regressions u and y. It is then interesting to see how to control the error of the fine approximations through the representation error like in [START_REF] Broadie | Risk estimation via regression[END_REF]. First, for s j < s k < s j and Borel measurable Θ function of (X m 0 ,m 1 s j ,s k :δ(s k ) ) with Θ(X m 0 ,m 1 s j ,s k :δ(s k ) ) integrable, we denote E m 0 ,x s j ,s k , E m 0 ,x s j ,s k and E m 0 ,x s j ,s k the operators defined by

E m 0 ,x s j ,s k (Θ(X m 0 ,m 1 s j ,s k :δ(s k ) )) = E m 0 s j Θ(X m 0 ,m 1 s j ,s k :δ(s k ) )|X m 0 ,m 1 s j ,s k = x , E m 0 ,x s j ,s k (Θ(X m 0 ,m 1 s j ,s k :δ(s k ) )) = t T m 0 s j ,s k (x -X m 0 s k ) Λ m 0 s j ,s k -1 M 1 M 1 m 1 =1 T m 0 s j ,s k X m 0 ,m 1 s j ,s k -X m 0 s k Θ(X m 0 ,m 1 s j ,s k :δ(s k ) )
and

E m 0 ,x s j ,s k (Θ(X m 0 ,m 1 s j ,s k :δ(s k ) )) = t T m 0 s j ,s k (x -X m 0 s k )R m 0 s j ,s k Θ(X m 0 ,m 1 s j ,s k :δ(s k ) ) with R m 0 s j ,s k Θ(X m 0 ,m 1 s j ,s k :δ(s k ) ) ∈ argmin r∈R d 1 E m 0 s j   E m 0 s j Θ(X m 0 ,m 1 s j ,s k :δ(s k ) )|X m 0 ,m 1 s j ,s k -t T m 0 s j ,s k (X m 0 ,m 1 s j ,s k -X m 0 s k )r 2   .
When E s j is the conditional expectation knowing X m 0 s j , E m 0 s j is the conditional expectation knowing the trajectory of X m 0 starting from X m 0 s j . E m 0 s j is used as the regression basis depends on X m 0 . For a given s j ∈ S, in contrast to expressions presented in sections 2.2 and 3.2, we simplify the presentation here and we omit to center the regressions around u m 0 ,S δ(s j ),δ(s j ) or y m 0 ,S δ(s j ),δ(s j ) . Consequently, the value of u m 0 ,S s j ,s j and y m 0 ,S s j ,s j are obtained through respectively averaging on u m 0 ,S s j ,δ(s j ) (x) = E m 0 ,x s j ,δ(s j ) (Θ u (X m 0 ,m 1 s j ,δ(s j ):δ 2 (s j ) )) and on y m 0 ,S s j ,δ(s j ) (x) = E m 0 ,x s j ,δ(s j ) (Θ y (X m 0 ,m 1 s j ,δ(s j ):δ 2 (s j ) )), where Θ u (X m 0 ,m 1 s j ,δ(s j ):δ 2 (s j ) ) = u m 0 ,S s j ,δ 2 (s j ) (X m 0 ,m 1 s j ,δ 2 (s j ) ) +

δ 2 (s j ) t l ≥δ(s j ) f (t l , X m 0 ,m 1 s j ,t l , X m 0 ,m 1 s j ,t l+1 ), Θ y (X m 0 ,m 1 s j ,δ(s j ):δ 2 (s j ) ) = y m 0 ,S s j ,δ 2 (s j ) (X m 0 ,m 1 s j ,δ 2 (s j ) ) +∆ δ(s j ) f (δ(s j ), y m 0 ,S s j ,δ 2 (s j ) (X m 0 ,m 1 s j ,δ 2 (s j ) ), z m 0 ,S s j ,δ(s j ) (X m 0 ,m 1 s j ,δ(s j ) )).
We assume Lipschitz condition uniformly in time of the driver f involved in (ODP ) with respect to its Y and Z coordinates. Although this condition is not necessary to obtain good numerical results in Section 5, it is required to apply Theorem 2 of [START_REF] Broadie | Risk estimation via regression[END_REF] (cf. Assumption F2 in [START_REF] Broadie | Risk estimation via regression[END_REF]) that yields the following asymptotical result. Proposition 4.1. Given that assumptions A1, A2 and A3 of [START_REF] Broadie | Risk estimation via regression[END_REF] are fulfilled and that the driver involved in (ODP ) is [f ] Lip -Lipschitz we have the following asymptotical inequality

( ρ -ρ) 2 ≤ [ρ] Lip E m 0 s j (κ 2 (X m 0 ,m 1 s j ,δ(s j ) )) + O p (1/M 1 ) (4.2)
as M 1 -→ ∞ where ( ρ, ρ, [ρ] Lip , κ) is either equal to ( ρ u , ρ u , [ρ] u Lip , κ u ) for (f ) or equal to ( ρ y , ρ y , [ρ] y
Lip , κ y ) for (ODP ) with ρ u = u m 0 ,S s j ,s j , ρ y = y m 0 ,S s j ,s j ,

ρ u = E m 0 s j   E m 0 ,X m 0 ,m 1 s j ,δ(s j ) s j ,δ(s j )
Θ u (X m 0 ,m 1 s j ,δ(s j ):δ 2 (s j ) ) + δ(s j )

t l ≥s j f (t l , X m 0 ,m 1 s j ,t l , X m 0 ,m 1 s j ,t l+1 )   , ρ y = E m 0 s j E m 0 ,X m 0 ,m 1 s j ,δ(s j ) s j ,δ(s j ) Θ y (X m 0 ,m 1 s j ,δ(s j ):δ 2 (s j ) ) +∆ s j f s j , E m 0 ,X m 0 ,m 1 s j ,δ(s j ) s j ,δ(s j ) Θ y (X m 0 ,m 1 s j ,δ(s j ):δ 2 (s j ) ) , z m 0 ,S s j ,s j , [ρ] u Lip = 1, [ρ] y Lip = 1 + ∆ s j [f ] Lip and κ u (x) = E m 0 ,x s j ,δ(s j ) Θ u (X m 0 ,m 1 s j ,δ(s j ):δ 2 (s j ) ) -E m 0 ,x s j ,δ(s j ) Θ u (X m 0 ,m 1 s j ,δ(s j ):δ 2 (s j ) ) , κ y (x) = E m 0 ,x s j ,δ(s j ) Θ y (X m 0 ,m 1 s j ,δ(s j ):δ 2 (s j )
) -E m 0 ,x s j ,δ(s j ) Θ y (X m 0 ,m 1 s j ,δ(s j ):δ 2 (s j ) ) .

Proposition 4.1 results from Theorem 2 and Remark 2 of [START_REF] Broadie | Risk estimation via regression[END_REF]; we expressed [ρ] Lip associated to each problem and we replaced E by E m 0 s j as the regression basis depends on X m 0 . Assumptions A1, A2 and A3 of [START_REF] Broadie | Risk estimation via regression[END_REF] are standard assumptions for regressions (cf. [START_REF] White | Asymptotic theory for econometricians[END_REF]). Considering the regression basis presented in Section 2.1 with E(|X t | 2 ) < ∞ for any t ∈ [0, T ], these assumptions are fulfilled if: i) the conditional variance of each regressed quantity is integrable and bounded from below by v 0 > 0, ii) the regression value is unbiased and iii) each component of the regression basis as well as κ (denoted M in [START_REF] Broadie | Risk estimation via regression[END_REF]) admit a finite fourth moment. When the latter moment assumption iii is needed to establish error control and can be modified using truncation (cf. [START_REF] Gobet | A regression-based Monte Carlo method to solve backward stochastic differential equations The[END_REF]), the further i&ii are sufficient to ensure the existence and uniqueness of the regressed representation.

In Proposition 4.1, we provided a control on fine approximations u and y. In Proposition 4.2, we rather focus on coarse approximations and decompose the conditional mean square error E m0 sj [u m0,S sj ,s k (X m0,m1 sj ,s k ) -u s k (X m0,m1 sj ,s k )] 2 into a bias term W, a variance term V and a regression error term R. Proposition 4.2. Assuming i and iii introduced above, for s j < s < s k taking their values in the discretization set S, we define

W m 0 ,S s j ,s k (x) = E m 0 s j u m 0 ,S s j ,s k (x) -u m 0 ,S s j ,s k (x) , R m 0 ,S s j ,s k (x) = E m 0 s j u m 0 ,S s j ,s k (x) -u s k (x) , V m 0 ,S s j ,s k (x) = Var m 0 s j u m 0 ,S s j ,s k (x) , with u m 0 ,S s j ,s k (x) = E m 0 ,x s j ,s k   u δ(s k ) (X m 0 ,m 1 s j ,δ(s k ) ) + δ(s k ) t l+1 >s k f (t l , X m 0 ,m 1 s j ,t l , X m 0 ,m 1 s j ,t l+1 )   then E m0 sj [u m0,S sj ,s k (X m0,m1 sj ,s k ) -u s k (X m0,m1 sj ,s k )] 2 = E m0 sj (V m0,S sj ,s k (X m0,m1 sj ,s k )) +E m0 sj [R m0,S sj ,s k (X m0,m1 sj ,s k ) + W m0,S sj ,s k (X m0,m1 sj ,s k )] 2
and there exists a positive constant K m 0 1,s j ,s k depending on the regression basis such that

E m0 sj [W m0,S sj ,s k (X m0,m1 sj ,s k )] 2 ≤ K m0 1,sj ,s k E m0 sj [u m0,S sj ,δ(s k ) (X m0,m1 sj ,δ(s k ) ) -u δ(s k ) (X m0,m1 sj ,δ(s k ) )] 2 .
Proof. As we simulate several independent copies of X m 0 ,m 1 , we make sure that the approximations u are independent from X m 0 ,m 1 conditionally on X m 0 . Then, the expansion of

E m0 sj [u m0,S sj ,s k (X m0,m1 sj ,s k ) -u s k (X m0,m1 sj ,s k )]
2 is gotten when we notice that

u s k (X m0,m1 sj ,s k ) = E m0,X m 0 ,m 1 s j ,s k sj ,s k   u δ(s k ) (X m0,m1 sj ,δ(s k ) ) + δ(s k ) t l+1 >s k f (t l , X m0,m1 sj ,t l , X m0,m1 sj ,t l+1 )   .
An expression for the constant K m 0 1,s j ,s k can be obtained after expanding

E m 0 s j [W m 0 ,S s j ,s k (X m 0 ,m 1 s j ,s k )] 2 using u m 0 ,S s j ,s k (x) = E m 0 ,x s j ,s k   u m 0 ,S s j ,δ(s k ) (X m 0 ,m 1 s j ,δ(s k ) ) + δ(s k ) t l+1 >s k f (t l , X m 0 ,m 1 s j ,t l , X m 0 ,m 1 s j ,t l+1 )   .
Finally, we should point out that one could establish a similar result for (ODP ). Indeed, for instance, using the following coarse discretization to approximate (ODP )

     Ŷs k = E s k f s k ( Ŷδ(s k ) , Ẑs k ) , Ẑs k = 1 ∆s k E s k Ŷδ(s k ) (W δ(s k ) -W s k ) , (4.3)
with f s k (y, z) = y + ∆ s k f s k (y, z), the bias is then controlled as follows

E m0 sj [ W m0,S sj ,s k (X m0,m1 sj ,s k )] 2 ≤ K m0 1,sj ,s k E m0 sj      f s k ( Y sj ,δ(s k ) (X m0,m1 sj ,δ(s k ) ), Z sj ,s k (X m0,m1 sj ,s k )) -f s k (y m0,S sj ,δ(s k ) (X m0,m1 sj ,δ(s k ) ), z m0,S sj ,s k (X m0,m1 sj ,s k ))   2    ,
for some positive constant K m 0 1,s j ,s k depending on the regression basis. Therefore, the bias upper bound depends heavily on the driver choice.

As shown in Proposition 4.2, the bias W at time step s k is controlled by the mean square error at time step δ(s k ) decomposed into a variance term V, a regression error term R and a bias term at time step δ(s k ). Thus, increasing the number of time steps loosens the bias control as it involves more and more terms. In some situations, this accumulation of errors is a source of a significant bias backpropagation. A standard method to reduce this bias backpropagation is to use a coarser set S. In this paper, we rather proposed a new trick to cut this bias backpropagation using the couple (s, s) for any s ∈ S. Like for the example of Section 5.4, the latter method is quite effective and will be further studied in CMCLDII (cf. [START_REF] Abbas-Turki | Conditional Monte Carlo Learning for Diffusions II: extended methodology and application to risk measures and early stopping problems[END_REF], Section 4).

Some numerical results

In this section we test the presented conditional MC learning procedure on various BSDE examples. The fact that the driver f depends also on X is not a burden to the use of our method. All simulations are run on a laptop that has an Intel i7-7700HQ CPU and a single GeForce GTX 1060 GPU programmed with the CUDA/C application programming interface. We refer the reader to [START_REF]Cuda C PROGRAMMING GUIDE[END_REF] for an introduction to CUDA programming.

The proposed 1NMC simulation combined with local regressions leads to a simulation with batch parallel regressions that can be successfully implemented using the contribution [START_REF] Abbas-Turki | Resolving small random symmetric linear systems on graphics processing units[END_REF]. We had even to extend the work presented in [START_REF] Abbas-Turki | Resolving small random symmetric linear systems on graphics processing units[END_REF] to examples with larger dimensions. It is worth mentioning that the very competitive execution times that we present hereafter are possible not only because of the GPU use but also due to the adaptation of parallelization strategies presented in [START_REF] Abbas-Turki | Resolving small random symmetric linear systems on graphics processing units[END_REF].

Allen-Cahn equation

We consider (ODP ) simulation as presented in Section 3.2, we use the following functions We would like to approximate the solution y(t, x) of the Allen-Cahn PDE defined as follows, y(T, x) = f (T, x), ∂y ∂t (t, x) + y(t, x) -[y(t, x)] 3 + (∆ x y)(t, x) = 0.

f (t, x, y, z) = y -y 3 , f (T, x) = 2 + 2 5 |x| 2 d 1 -1 and E t k (x, w) = x + √ 2w, X t 0 = 0.
(5.1)

A benchmark approximation y b (0, x) for the solution y(0, x) of the PDE (5.1) is given in [Section 4.2; [START_REF] Han | Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations[END_REF]].

Table 1 provides the y(0, 0) simulated value of equation (5.1), approximated by Y learn 0 and Y sim 0 expressions, with respect to the number of inner trajectories M 1 . The benchmark solution y b (0, 0) is equal to 0.0528 for T = 0.3 and d 1 = 100. The standard deviation of each expression and the runtime in seconds are also given. We reduce the bias by increasing the number of inner trajectories. Table 1 shows that a relative small number of outer and inner trajectories is sufficient to observe a small variance and bias for both options. In fact, we show that the standard deviation is already acceptable for M 0 = 2 4 outer trajectories and the bias is acceptable for M 1 = 2 6 inner trajectories with an execution time of 56 millisecond.

The average learned Y learn 0 and the simulated Y sim 0 values in Table 2 are presented for the longer time horizon T = 1. The benchmark solution is equal to 0.0338 for T = 1 and d 1 = 100. To achieve a similar level of variance and bias, we need larger number of outer and inner trajectories than in Table 1. In fact, for M 0 = 2 5 of outer trajectories and M 1 = 2 6 of inner trajectories, we obtained an acceptable bias and standard deviation within 4 seconds of execution.

Multidimensional Burgers-type PDEs

We assume the (ODP ) setting presented in Section 3.2, we use the following functions

f (t, x, y, z) = y - 2 + d 1 2d 1 d 1 i=1 z i , f (T, x) = exp T + 1 d 1 d 1 i=1 x i 1 + exp T + 1 d 1 d 1 i=1 x i and E t k (x, w) = x + d 1 √ 2 w, X t 0 = 0.
We simulate the solution y(t, x) of the multidimensional Burgers-type PDE (cf [START_REF] Chassagneux | Linear multistep schemes for BSDEs[END_REF], Example 4.6) defined as follows, y(T, x) = f (T, x),

∂y ∂t (t, x) + d 1 2 2 (∆ x y)(t, x) + y(t, x) - 2 + d 1 2d 1 d 1 d 1 i=1 ∂y ∂x i (t, x) = 0. (5.2) PDE (5.
2) admits an explicit solution, we refer the reader to [Lemma 4.3, [START_REF] Han | Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations[END_REF]] for more details. The value of the solution y(0, 0) is 0.5000 for T = 0.2 and d 1 = 100. 4 is augmented till reaching a sufficiently small standard deviation within less than 7 minutes of execution time.

Time-dependent reaction-diffusion-type example

Let κ = 0.6, λ = 1 √ d 1 , we use the following functions

f (t, x, y, z) = min 1, y -κ -1 -sin λ d 1 i=1 x i e λ 2 d(t-T ) 2 , f (T, x) = 1 + κ + sin λ d 1 i=1 x i and E t k (x, w) = x + w, X t 0 = 0.
We simulate the solution y(t, x) of the time dependent reaction-diffusion-type PDE (cf [START_REF] Gobet | Adaptive importance sampling in leastsquares Monte Carlo algorithms for backward stochastic differential equations[END_REF], Section 6) defined as follows, y(T, x) = f (T, x),

∂y ∂t (t, x) + min 1, y -κ -1 -sin λ d 1 i=1 x i + 1 2 (∆ x y)(t, x) = 0. ( 5.3)
The explicit solution of the PDE (5.3) is given in [Lemma 4.4; [START_REF] Han | Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations[END_REF]].

Table 5 shows the approximated solution of the equation (5. Y sim 0 std 2 5 1.8197 (± 0.0386) 1.7587 (± 0.0287) 244 2 6 1.7125 (± 0.0104) 1.6799 (± 0.0116) 311 2 7 1.6605 (± 0.0037) 1.6376 (± 0.0091) 466 2 8 1.6458 (± 0.0023) 1.6290 (± 0.0089) 817 2 9 1.6439 (± 0.0019) 1.6283 (± 0.0061) 1526

Due in part to small number of time steps (L = 3), we obtain very accurate solutions within less than 2 seconds of execution.

A Hamilton-Jacobi-Bellman (HJB) equation

We assume here the driver to be equal to

f (t, x, y, z) = -|z| 2 d 1 , f (T, x) = ln 1 + |x| 2 d 1 /2 and E t k (x, w) = x + √ 2w, X t 0 = 0.
We calculate the solution y(t, x) of the HJB equation (cf [START_REF] Chassagneux | Numerical simulation of quadratic BSDEs[END_REF] Section 4.2) defined by y(T, x) = f (T, x),

∂y ∂t (t, x) + (∆ x y)(t, x) -|(∇ x y)(t, x)| 2 d 1 = 0.
(5.4) PDE (5.4) admits a benchmark solution. We refer the reader to [Lemma 4.2; [START_REF] Han | Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations[END_REF]] for more details.

Related to inequality (3.19), Figure 2 shows between two approximations 1 M 0

M 0 m 0 =1
y m 0 ,S s,s and 1 M 0

M 0 m 0 =1
y m 0 ,S δ(s),δ(s) + (δ(s) -s)f (s, y m 0 ,S δ(s),δ(s) , z m 0 ,S s,s ) of the same quantity with respect to the time discretization. The purpose here is to show the benefits of using the couple (s, s) = (δ(s), s) to stop the bias backpropagation. On Figure 2 left, we perform the conditional MC procedure taking s = T for all s ∈ S = {0, Figure 3 shows the convergence of the expressions (3.20) and (3.21) to the benchmark value with respect to the number of inner trajectories. In particular, we observe that both expressions converge to the benchmark solution with a small variance when M 1 = 2 17 . x i ∂y ∂x i (t, x) -y(t, x) , R l d1 i=1

x i ∂y ∂x i (t, x) -y(t, x) } = 0. PDE (5.5) has a benchmark solution given in [Section 4.4; [START_REF] Han | Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations[END_REF]]. This benchmark solution is equal to 21.299 for T = 0.5 and d 1 = 100.

Figure 4 shows the approximation of the solution of PDE (5.5) with respect to the number of inner trajectories. As before, we illustrate the values given by both (3.20) and (3.21) expressions. We show that 2 7 outer trajectories and 2 11 inner trajectories are sufficient to get an accurate approximation of the solution as the obtained values are in the corridor of the standard deviation of the benchmark solution given in [Section 4.4; [START_REF] Han | Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations[END_REF]]. No bias cut is needed here and s k = T for any s k ∈ S. The runtime with 2 7 outer trajectories and 2 11 inner trajectories is 53 seconds. Figure 6 shows the numerical solution y(0, 0) of (5.6) with respect to various choices of coarse time steps. L = 8 is sufficient to discretize the problem when the time horizon T = 1. The latter simulation is achieved in 620 seconds of runtime.

.

  the simulation Y m 0 s and Z m 0 s of Y and Z respectively at X m 0 s are set to be equal to y m 0 ,S s,s k and to z m 0 ,S s,δs k (s k )

) and the simulated values Y sim 0 and Z sim 0 of Y 0

 00 and Z 0 are respectively equal to

1

 1 3), calculated by the average learned (3.20) and the simulated (3.21) expressions, with respect to the number of inner trajectories M 1 . The standard deviation of each expression and the runtime in milliseconds are also given. The benchmark solution is equal to 1.6000 for T = 1, d 1 = 100. Table 5: Numerical simulations for PDE (5.3): T = 0.5, d 1 = 100, M 0 = 2 10 , L = 3; [Benchmark solution] y b (0, 0) = 1.6000. M Average learned Simulated Runtime in sec. (10 -3 ) Y learn 0 std

  Figure 2 right, we take s = s + 3 8 ∧ T when s ∈ S. The latter choice of s is then clearly justified, it can be even set automatically during the execution of the simulation using inequality (3.19) without any manual intervention.

Figure 2

 2 Figure 2: 1 M 0 M 0

Figure 3 :

 3 Figure 3: Numerical solution of PDE (5.4) calculated by (3.20) and (3.21) expressions: T = 1, d 1 = 100, M 0 = 2 7 , L = 3.

5. 5 1 i=1zσ d 1 i=1z 1 x 1 x

 51111 Pricing of European financial derivatives with different interest rates for borrowing and lendingAssuming µ = 0.06, σ = 0.2, R l = 0.04 and R b = 0.06, we introduce the following functionsf (t, x, y, z) = -R l y -(µ -R l ) σ d i + (R b -R l ) max {0, 1 i -y}, f (T, x) = max { max 1≤i≤d i -120, 0} -2 max { max 1≤i≤d i -150, 0}andE t k (x, w) = x exp (µ -σ22 )∆ t + σw , X t 0 = 100.Let y defined as the solution of the following PDE, y(T, x) = f (T, x),

Figure 4 :

 4 Figure 4: Numerical solution of PDE (5.5) calculated by the average learned and the simulated expressions: T = 0.5, d 1 = 100, M 0 = 2 7 , L = 2.

5. 6 A 1 α 1 , 1 .Figure 5 [

 61115 Figure 5 [Right], the distribution of y 253 256 , X m 0 253 256 is either simulated using Y m 0 253 256that in-

Table 1 :

 1 Numerical simulations for PDE (5.1): T = 0.3, M 0 = 2 4 , d 1 = 100, L = 4; [Benchmarck solution] y b (0, 0) = 0.0528.

	M 1	Average learned	Simulated	Runtime in sec. (10 -3 )
		Y learn 0	std	Y sim 0	std
	2 4 0.0454 (± 0.0093) 0.0455 (± 0.0073)	13
	2 5 0.0513 (± 0.0011) 0.0517 (± 0.0008)	23
	2 6 0.0523 (± 0.0004) 0.0518 (± 0.0006)	56
	2 7 0.0526 (± 0.0003) 0.0515 (± 0.0001)	119
	2 8 0.0525 (± 0.0002) 0.0517 (± 0.0002)	227
	2 9 0.0527 (± 0.0002) 0.0515 (± 0.0002)	414

Table 2 :

 2 Numerical simulations for PDE (5.1): T = 1, d 1 = 100, M 0 = 2 5 , L = 6; [Benchmarck solution] y b (0, 0) = 0.0338.

	M 1	Average learned	Simulated	Runtime in sec.
		Y learn 0	std	Y sim 0	std
	2 5 0.0345 (± 0.0008) 0.0350 (± 0.0021)	2
	2 6 0.0333 (± 0.0003) 0.0326 (± 0.0004)	4
	2 7 0.0334 (± 0.0002) 0.0330 (± 0.0003)	7
	2 8 0.0336 (± 0.0002) 0.0332 (± 0.0002)	12
	2 9 0.0336 (± 0.0001) 0.0331 (± 0.0001)	27

Table 3 :

 3 Numerical simulations for PDE (5.2): T = 0.2, d 1 = 100, M 0 = 2 6 , L = 5; [Explicit solution] y(0, 0) = 0.5000.

	M 1	Average learned	Simulated	Runtime in sec.
		Y learn 0	std	Y sim 0	std
	2 8 0.4785 (± 0.0428) 0.0517 (± 0.0431)	7
	2 9 0.5113 (± 0.0450) 0.5108 (± 0.0450)	16
	2 10 0.4966 (± 0.0448) 0.4912 (± 0.0447)	27
	2 11 0.5022 (± 0.0421) 0.5012 (± 0.0435)	49

Table 4 :

 4 Numerical simulations for PDE (5.2): T = 0.2, d 1 = 100, M 1 = 2 11 , L = 5; [Explicit solution] y(0, 0) = 0.5000.

	M 0	Average learned	Simulated	Runtime in sec.
		Y learn 0	std	Y sim 0	std
	2 5 0.4953 (± 0.0618) 0.4941 (± 0.0615)	24
	2 6 0.5022 (± 0.0424) 0.5013 (± 0.0435)	49
	2 7 0.5079 (± 0.0346) 0.5066 (± 0.0342)	103
	2 8 0.5158 (± 0.0221) 0.5151 (± 0.0221)	194
	2 9 0.5023 (± 0.0164) 0.5029 (± 0.0164)	408

Table 3

 3 presents the approximations (3.20) and (3.21) of the equation (5.2), with respect to the number of inner trajectories M 1 . It shows that the standard deviation of both results should be reduced by increasing the number of outer trajectories. The number M 0 of outer trajectories in Table

called "Trained" vs. its "Tested" counterpart used for the ex-

The obtained distributions are quite similar which strengthen the possibility of using (3.12) not only for s = s k but for various intermediary discretization times s ∈ {s k , s k + ∆ t , ..., δ s j (s k ) -∆ t } (cf. CMCLDII [START_REF] Abbas-Turki | Conditional Monte Carlo Learning for Diffusions II: extended methodology and application to risk measures and early stopping problems[END_REF]).