Variations on branching methods for non linear PDEs
Xavier Warin













 


callettersOMScmsymn
calletters





.equation

    TheoremTheorem[section]
    Definition[Theorem]DefinitionProposition[Theorem]PropositionPropertyPropertyAssumption[Theorem]AssumptionAlgorithmAlgorithmLemma[Theorem]LemmaCorollary[Theorem]CorollaryRemark[Theorem]RemarkExample[Theorem]ExampleProofProofIntroductionIntroductionaddtoresetequationsection.equation
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    -0.1 0.2 -0.1 0.2 -0.02   
    
    
    
    𝔻𝔼𝔽Ł𝕃ℙℚℝ𝕊𝕄ℕ
     A B C D E F G H I J K L M N O P Q R S T U V W X Y ZLD̅CCOBÑB̃p̃ỸZ̃τ̅θ̅ψψ M I
    αασαΩωΩτ̂ε𝐱0h̋XΛa̅Σ F𝔽ψφ
    
    ×⊗ 1 Proof. ⟨⟩FLp̂IXWT^1^-⌊⌋Variations on branching methods for non linear PDEs
         Xavier Warin EDF R&D & FiME, Laboratoire de Finance des Marchés de l'Energie, xavier.warin@edf.fr
        December 30, 2023
    ======================================================================================================
    
     The branching methods developed in <cit.>, <cit.> are effective methods to solve some  semi linear PDEs and are shown
    numerically to be  able to solve some full non linear PDEs.
    These  methods are however restricted to some small coefficients in the PDE and small maturities.
    This article shows numerically that these methods can be adapted to solve the problems with longer maturities in the semi-linear case by using a new derivation scheme and some nested method.
    As for the case of full non linear PDEs, we introduce new schemes  and we show numerically that they provide an effective alternative to the schemes previously developed.
    
    
    § INTRODUCTION
    
    The resolution of low dimensional non linear PDEs  is often achieved by some deterministic methods such as finite difference schemes, finite elements and finite volume.
    Due the curse of dimensionality, these methods cannot be used in dimension greater than three : both the computer time and the memory required are too large even for supercomputers.
    In the recent years the probabilistic community has developed some representation of semi linear PDE:
    
    
    
     -∂_tu- u =f(u,Du),             u_T=g,    t<T, x∈^d,
    
    
    by means of backward stochastic differential equations (BSDE), as introduced by  <cit.>.
    Numerical Monte Carlo algorithms have been developed to solve efficiently these BSDE by  <cit.>,  <cit.>.
    The representation  of  the following full non linear PDE: 
    
     
    -∂_tu- u =f(u,Du,D^2u),             u_T=g,    t<T, x∈^d,
    
    
    has been given by the mean of second order backward stochastic differential  equation (SOBSDE) by  <cit.>.
    A numerical algorithm  developed by <cit.> has been derived  to solve these  full non linear PDE by the mean of SOBSEs.
    
    The BSDE and SOBSDE schemes developed rely on the approximation of conditional expectation and the most effective implementation is based on regression methods as developed in  <cit.>, <cit.>.
    These regression methods develop an approximation of   conditional expectations based on an expansion  on basis functions. The size of this expansion has to grow exponentially with the dimension of the problem so we have to face again the curse of dimensionality. 
    Notice that the BSDE methodology could be used in dimension 4 or 5 as regressions has been successfully used in dimension 6 in   <cit.> using some local regression function.
    
    Recently a new representation of semi linear equations (<ref>) for a polynomial function f of u and Du has been given by <cit.> : this representation
    uses the automatic differentiation approximation  as used in <cit.>, <cit.>, <cit.>, and <cit.>. 
    The authors have shown that the representation gives a finite variance estimator only for small maturities or small non linearities and numerical  examples until dimension 10 are given. Besides, they have shown that the given scheme using Malliavin weights cannot be used to solve the full non linear equation (<ref>).
    
    <cit.> have introduced a re-normalization technique improving numerically the convergence of the scheme diminishing the variance observed for the semi linear case. Besides, the authors haved introduced a scheme to solve the full non linear equation (<ref>). Without proof of convergence they numerically have shown that the developed scheme is effective.
    
    The aim of the paper is to provide some numerical variation on the algorithm developed in  <cit.>.
    In a first part we will show, with simple ideas, that it is possible to deal with longer maturities than the ones possible with the initial algorithm.
    
    In a second part we give some alternative schemes to the one proposed in <cit.> and, testing  them on some numerical examples, we show that they are superior than the scheme previously developed.
    
    In the numerical results presented in the article, all errors are estimated  as  the log of the  standard deviation observed divided  by the square root of the number of particles used and these errors are plotted as a function of the log of the number of particles used. As our methods are pure Monte Carlo methods we expect to have lines with slope -1/2 when the numerical variance is bounded. 
    
     
    
    
    § THE SEMI LINEAR CASE
    
    Let σ_0 ∈^d be some constant non-degenerate matrix,  μ∈^d be some constant vector  f: [0,T] ^d ^d →,
    and g: ^d → bounded Lipschitz functions, 
    we consider the semi linear parabolic PDE:
    ∂_t u + 1/2σ_0 σ_0^⊤ : D^2 u + μ.Du + f(·, u ,Du) = 0,
         [0,T) ^d,
    
    with terminal condition u(T, ·) = g(·) where 	A:B :=  (AB^⊤) for two matrices A, B ∈^d.
    
    When f is a polynomial in (u, Du) in the form
    
    f(t,x,y,z,γ)
     =∑_ℓ = (ℓ_0, ℓ_1, ·, ℓ_m) ∈ L c_ℓ(t,x) y^ℓ_0∏_i=1^m (b_i · z)^ℓ_i ,
    
    for some m ≥ 1, L⊂^1 + m, 
    where (b_i)_i=1, m is a sequence of ^d-valued bounded continuous functions defined on [0,T] ^d, 
    and (c_ℓ)_ℓ∈ L is a sequence of bounded continuous functions defined on [0,T] ^d.
     <cit.> obtained a probabilistic representation to the above PDE by branching diffusion processes under some technical conditions.
    In the sequel, we simplify the setting by taking f  as a constant (in u, Du) plus a monomial  in u, (b_i.Du) , i=1,m : 
    
    f(t,x,y,z)
     =
     h(t,x)+ c(t,x) y^ℓ_0∏_i=1^m (b_i · z)^ℓ_i ,
    
    for some m ≥ 1,
    where (b_i)_i=1, m  is a sequence of ^d-valued bounded continuous function defined on [0,T] ^d, (ℓ_i)_i=0,m∈^m+1 supposing that ∑_i=0,mℓ_i >0,
    and c is a  bounded continuous function defined on [0,T] ^d. We note L= ∑_i=0^m ℓ_i.
    
    The case with f a general polynomial only complexifies the notation : it can be simply treated as in <cit.> by introducing some probability mass function (p_ℓ)_ℓ∈ L (i.e. p_ℓ≥ 0 and ∑_ℓ∈ L p_ℓ = 1) that are used to select with monomial to consider during the branching procedure. Another approach can be used : instead of sampling the monomial to use, it is possible to consider successively  all terms of the f but this doesn't give a representation as nice as the one in  <cit.>.
    
    
     §.§ Variation on the original scheme of <cit.>
    
    In this section we present the original scheme of <cit.> and explain how to diminish the variance increase the maturities of the problem.
    
    
      §.§.§ The branching process
    
    Let us first introduce a branching process with arrival time of distribution density function ρ.  At the arrival time, the particle branches into |ℓ | offsprings.
    We introduce a sequence  of   i.i.d. positive random variables (τ^k)_k = (k_1, ⋯, k_n-1, k_n) ∈^n, n>1 with all the values k_i ∈ [1,L], for  i >0.
    
    We construct an age-dependent branching process using the following procedure :
    
      *   We start from a particle marked by 0, indexed by (1), of generation 1,
    		whose arrival time is given by T_(1) := τ^(1)∧ T.
    
      *   Let k = (k_1, ⋯, k_n-1, k_n) ∈^n be a particle of generation n, with arrival time T_k that branches into L offspring particles noted (k_1, ⋯, k_n-1, k_n,i) for i=1,...,L. We  define the set of its offspring particles by
    	
        S(k) := {(k_1, ⋯, k_n, 1), ⋯, (k_1, ⋯, k_n, L) },
    
      We first mark  the ℓ_0 particles by 0,  the ℓ_1 next by 1 , and so on,  so that each particle has a mark i for i = 0, ⋯, m.
    
      *  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- + τ^k) ∧ T. Let us denote Δ T_k = T_k -T_k-.
    
      *  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.
    
    We denote further
    
        θ_k :=  k,
           ^n_t := { k   n  T_k-≤ t < T_k },      t ∈ [0,T),
        {k    n   T_k = T},       t = T,
    
    and also
    
        ^n_t := ∪_s ≤ t^n_s,
            _t := ∪_n ≥ 1^n_t
               _t := ∪_n ≥ 1^n_t.
    
    Clearly, _t (resp. ^n_t) denotes the set of all living particles (resp. of generation n) in the system at time t,
    and _t (resp. ^n_t) denotes the set of all particles (resp. of generation n) being alive at or before time t.
    
    
    
    We next equip each particle with a Brownian motion in order to define a branching Brownian motion.
    Let (Ŵ^k)_k = (k_1, ⋯, k_n-1, k_n) ∈^n, n>1 be a sequence of independent d-dimensional Brownian motion, which is also independent of (τ^k)_k = (k_1, ⋯, k_n-1, k_n) ∈^n, n>1.
    Define W^(1)_t =  Ŵ^(1)_t for all t ∈[0, T_(1)] and then for each k = (k_1, ⋯, k_n) ∈_T ∖{(1)},
    define
    
    W^k_t  :=  W^k-_T_k- + Ŵ^k_t - T_k-,     t ∈ [T_k-, T_k].
    
    Then (W^k_·)_k ∈_T is a branching Brownian motion.
    
    
      §.§.§ The original algorithm
    
    Let us denote  F̅(t):=∫_t^∞ρ(s)ds.
    Denoting X^k_t := x + μ t + σ_0 W^k_t for all k ∈_T and t ∈ [T_k-, T_k] and
    by _t,x the expectation operator conditional on the starting data X_t=x at time t, 
     we obtain from the Feynman-Kac formula the representation of the solution u of equation (<ref>) as:
     
     u(0,x)
     =
     _0,x[F̅(T)g(X_T)/F̅(T)+∫_0^T f(u,Du)(t,X_t)/ρ(t)ρ(t)dt]
     =
     _0,x[ϕ(T_(1),X^(1)_T_(1))],
     
    where T_(1):=τ^(1)∧ T, and
      ϕ(t,y)
      := _{t≥ T}/F̅(T)g(y)
                        +_{t<T}/ρ(t)( h+c u^ℓ_0∏_i=1^m (b_i · Du)^ℓ_i)(t,y).
     
    On the event {_{T_(1)<T}}, using the independence of the (τ^k,W^k)   we are left to calculate
    
    
    [c u^ℓ_0   ∏_i=1^m (b_i · Du)^ℓ_i](T_(1),X_T_(1))  = c ∏_j=1^ℓ_0_T_(1),X_T_(1)[ ϕ(T_(1,j),X^(1)_T_(1,j))]
       ∏_i=1^m ( b_i(T_(1),X_T_(1)).D _T_(1),X_T_(1)[ϕ(T_(1,p),X^(1,p)_T_(1,p))])^ℓ_i
    Using differentiation with respect to the heat kernel, i.e. the marginal density of the Brownian motion we get :
    
    [c u^ℓ_0∏_i=1^m (b_i · Du)^ℓ_i](T_(1),X_T_(1))  = c ∏_j=1^ℓ_0_T_(1),X_T_(1)[ ϕ(T_(1,j),X^(1,j)_T_(1,j))]  
        ∏_i=1^m ( b_i(T_(1),X_T_(1))._T_(1),X_T_(1)[ (σ_0^⊤)^-1Ŵ^(1,p)_Δ T_(1,p)/Δ T_(1,p)ϕ(T_(1,p),X^(1,p)_T_(1,p))])^ℓ_i 
    Using equations (<ref>) and (<ref>) recursively and the tower property , we get the following
    representation 
    
    u(0,x)
     = _0,x[  _(1)]
    
    where _(1) is given by the backward recursion :
    let  _k := g(X^k_T) -g(X^k_T_k-)_{θ_k ≠ 0}/(Δ T_k) for every k ∈_T, then let
    _k
     := 1/ρ(Δ T_k)( h(T_k,X^k_T_k) + c(T_k, X^k_T_k)
    ∏_k̃∈ S(k)_k̃_k̃),
         k ∈_T ∖_T.
    
    
    where 
    _k
    		  = _{θ_k = 0} + _{θ_k ≠ 0}  b_θ_k(T_k-, X^k_T_k-)
    		· (σ_0^⊤)^-1Ŵ^k_Δ T_k/Δ T_k.
    
    and we have used that  _0,x[ g(X^k_T_k-)  b_θ_k(T_k-, X^k_T_k-)
    		·σ_0^⊤)^-1Ŵ^k_Δ T_k]  =0.
    This backward representation is slightly different from the  elegant representation  introduced  in <cit.>.
    Clearly on our case the variance of the method used will be lower than with the representation   in <cit.> for a similar computational cost.
    
    
    In the case where the operator f is linear  and a function of the gradient (ℓ_0=0, m=1 and ℓ_1=1) using the arguments in <cit.>  it can be easily seen by conditioning with respect to the number of branching that equation (<ref>) is of finite variance if 1/x ρ(x)^2 = O(x^α) as x ⟶ 0 with α≥ 0.
    
    When τ follows for example a gamma law with parameters κ and θ, the finite variance is proved as soon as κ≤ 0.5 for PDE coefficients and maturities small enough.
    
    In the non linear case, <cit.> have shown that the variance is in fact finite  for maturities small enough and small coefficients  as soon as κ <0.5 but numerical results show that κ=0.5 is optimal in term of efficiency: for  a given θ the numerical variance is nearly the same for the values of κ between 0.4 and 0.5 but a higher κ value limits the number of branching thus meaning a smaller computational cost.
    
    
      §.§.§ Variation on the original  scheme
    
    As indicated in the introduction, the method is restricted to small maturities or small non linearities.
    Having a given non linearity we are interested in adapting the methodology in  order to be able to treat longer maturities. A simple idea consists in noting that the Monte Carlo  method is applied by sampling the conditional expectation _t,x for t>0 appearing in equation (<ref>) only once. Using nested Monte Carlo, so by sampling each term of equation (<ref>) more that one time
    one can expect a reduction in the variance observed. A nested method of order n is defined as a method using n sampling to estimate each function u or Du at each branching.
    Of course the computational time will grow exponentially with the number of samples taken and for example 
    trying to use a gamma law  with a non linearity of Burger's type  u(b.Du)  with κ =0.5  is very costly: due to the  high values of the density ρ near 0, trajectories can have many branching.
    
    Some different strategies have been tested to  be able to use this technique :
    
      *  A first possibility consists in trying  to re-sample more at the beginning of the resolution and decreasing the number of samples as time goes by or as the number of branching increases.  The methodology works  slightly better than  a re-sampling with a constant number of particles but  has to be adapted to each  maturity and each case so it has been given up.
    
      *   Another observation is that the gamma law is only necessary  to treat the gradient term: so it is possible to use two laws: a first one, an exponential law,  will be used to estimate the u function while an gamma law will be used for the Du terms.  This second technique is the most effective and is used for the results obtained in the section.
    
    For a given dimension d , we take σ_0 = 1/√(d)_d, μ =,
    
     f(t,x,y,z)= d(t,x) +  y (b · z),
    
    where b := 0.2/d (1+1/d, 1+2/d, ⋯, 2) and 
    
    h(t,x) 
    :=
    cos( x_1 +⋯ + x_d) 
    ( α + σ_0^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 semi linear PDE (<ref>) is given by
    
        u(t,x) =  cos(x_1+ ⋯ + x_d) e^α (T-t).
    
    Our goal is to estimate u at t=0, x=0.5.
    This test case will be noted test A in the sequel. 
    
    We use the nested algorithm with two distributions for τ:
    
      *  an exponential law with density ρ(s)= λ e^-λ s with λ=0.4 to calculate the u terms,
    
      *  a gamma distribution ρ(s)  = 1/Γ(κ) θ^κ s^κ -1exp(- s/θ) _{s > 0} with Γ(κ) := ∫_0^∞ s^κ-1 e^-s ds and the parameters κ=0.5, 1/θ=0.4 to calculate the Du terms.
    
     
    We first give on figures <ref>,  <ref> and  <ref>  the results obtained for test A for different maturities and a dimension d=4 so the analytical solution is -0.508283. We plot for each maturity :
    
      *  the solution obtained by increasing the number of Monte Carlo scenarios used,
    
      *  the error calculated as explained in the introduction.
    
    Nested n  curves stand for the curves using  the nested method of order  n, so the  Nested 1 curve stands for the original method.
    
    
    
    On figure <ref>, for maturity 2.5  the error observed with the orignal method (Nested 1) is around 1000 so it has not been plotted.
    
    
    Because of the number of branching due to the gamma law, it seems difficult to use a nested method of order n>2 for long maturities : the time needed explodes.
    But clearly the nested method permits to have accurate solution for longer maturities.
    
    For a maturity of 2 we also give the results obtained in dimension 6 on figure <ref> giving an analytical solution  -1.4769 : once again the original method fails to converge while the nested one give good results.
    
    
    
     §.§ Adaptation of the original branching to the re-normalization technique
    
    As introduced in <cit.>, we introduce a modification of the original branching process that let us use exponential laws  for the branching dates 
    to treat the Du terms in the method previously described.
     	Recall that ^1_T = {(1)}, we introduce an associated ghost particle, denoted by (1^1),
    	and denote ^1_T := {(1), (1^1) }.
    	Next, given the collection ^n_T of all particles (as well as ghost particles) of generation n,
    	we define the collection ^n+1_T as follows.
    	For every k = (k_1, ⋯, k_n) ∈^n, we denote by o(k) = (k̂_1, ⋯, k̂_n) its original particles,
    	where k̂_i := j when k_i = j    j^1.
    	Further, when k = (k_1, ⋯, k_n) is such that k_n ∈, we denote k^1 := (k_1, ⋯, k_n-1, k_n^1).
    	The mark of k ∈^n will be the same as its original particle o(k), i.e. θ_k := θ_o(k);
    	and T_k := T_o(k), Δ T_k := Δ T_o(k) and τ^k = τ^o(k).
    	Define also ^n_T:= {k ∈^n_T  : o(k) ∈_T}.	
    	For every k = (k_1, ⋯, k_n) ∈^n_T ∖^n_T, we still define the set of its offspring particles by
    	
        S(k) := {(k_1, ⋯, k_n, 1), ⋯, (k_1, ⋯, k_n, L) },
    
    	and the set of ghost offspring particles by
    	
        S^1(k)
        		  :=  { (k_1, ⋯, k_n, 1^1), ⋯, (k_1, ⋯, k_n, L^1) }.
    
     	Then the collection ^n+1_T of all particles (and ghost particles) of generation n+1 is
    	
        ^n+1_T
        		 := ∪_k ∈^n_T ∖^n_T( S(k) ∪ S^1(k) ).
    
    	Define also
    	
        _T := ∪_n ≥ 1^n_T,
        		    _T := ∪_n ≥ 1^n_T.
    
    
      §.§.§ The original re-normalization technique
    
            We next equip each particle with a Brownian motion in order to define a branching Brownian motion.
    	Further, let W^∅_0 := 0, and for every k = (k_1, ⋯, k_n) ∈^n_T, let
    	
    		W^k_s
    		 := 
    		W^k-_T_k- + _k_n ∈Ŵ^o(k)_s - T_k-,
    		     
    		X^k_s := μ s +σ_0  W^k_s,
    		   ∀ s ∈ [T_k-, T_k].
    	
    
    On figure <ref>, we give the original Galton-Watson tree and the ghost particles associated.
    
    
    
    	
    
    The initial equation (<ref>) remains unchanged (first  step of the algorithm) but equation (<ref>) is modified
    by replacing the term 
    _T_(1),X_T_(1)[ (σ_0^⊤)^-1Ŵ^(1,p)_Δ T_(1,p)/Δ T_(1,p)ϕ(T_(1,p),X^(1,p)_T_(1,p))  ]
    
    by
    _T_(1),X_T_(1)[ (σ_0^⊤)^-1Ŵ^(1,p)_Δ T_(1,p)/Δ T_(1,p)(ϕ(T_(1,p),X^(1,p)_T_(1,p))  -ϕ(T_(1,p),X^(1,p^1)_T_(1,p)) )  ].
    
    Notice that since W^(1,p^1) has been obtained  by (<ref>), Ŵ^(1,p)_Δ T_(1,p) and ϕ(T_(1,p),W^(1,p^1)_T_(1,p)) are orthogonal so that adding the second term acts as a control variate.
    Recursively using the modified version of equation (<ref>) induced by the use of (<ref>), <cit.> gave defined the re-normalized estimator by a backward induction:
    let _k := g(X^k_T)/(Δ T_k) for every k ∈_T, then let
    _k
    :=
    1/ρ(Δ T_k) ( h(T_k, X^k_T_k)+c(T_k, X^k_T_k)
    ∏_k̃∈ S(k)(_k̃ -  _k̃^1_{θ(k̃) ≠ 0}) _k̃),
    		   k ∈_T ∖_T.
    
    where the weights  are given by equation (<ref>),
    so we have 
    
    u(0,x)
     = _0,x[  _(1)].
    
    
    As explained in section <ref>, equation (<ref>) used in representation (<ref>) force us to take laws for branching dates  with a high probability  of low values that leads to a high number of recursions defined by equation  (<ref>).
    Besides such laws using some rejection algorithm, as gamma laws, are very costly to generate.
    The use of  (<ref>) permits us to use exponential laws very cheap to simulate and with a low probability of small values. 
    
    Indeed it can be easily seen in the linear case  (f  function of the gradient with ℓ_0=0, m=1 and ℓ_1=1) by conditioning with respect to the number of branching that the variance is bounded for small maturities
    and coefficients if
    _0,x[ (_k -  _k^1_{θ(k) ≠ 0})^2 (b_θ_k(T_k-, X^k_T_k-)
    		· (σ_0^⊤)^-1Ŵ^o(k)_Δ T_k)^2/(Δ T_k)^2] < ∞.
    
    By X^k^1_t construction using g regularity, it is easily seen that for small time steps Δ T_k,  _0,x,Δ T_k[ (_k -  _k^1)^2 ] = O(Δ T_k) as Δ T_k ⟶ 0 and (<ref>) is satisfied for every ρ densities.
    
    
      §.§.§ Re-normalization techniques and antithetic
    
    We give a version of the re-normalization technique using antithetic variables.
    Equation (<ref>) is modified  by :
    
    W^k_s
     := 
    W^k-_T_k- + _k_n ∈Ŵ^o(k)_s - T_k- - _k_n ∉Ŵ^o(k)_s - T_k-, 
         
    		X^k_s := μ s +σ_0  W^k_s,
    		   ∀ s ∈ [T_k-, T_k],
                    
     for every k = (k_1, ⋯, k_n) ∈^n_T.
    
    Then equation (<ref>) is modified by  :
    
      *  First , replacing the term tacking into account the power of uϕ(T_(1,j),X^(1,j)_T_(1,j))
    
     by
    1/2( ϕ(T_(1,j),X^(1,j)_T_(1,j))+
    ϕ(T_(1,j),X^(1,j^1)_T_(1,j)) ),
    
      *  and the term taking into account the gradient 
    _T_(1),X_T_(1)[ (σ_0^⊤)^-1Ŵ^(1,p)_Δ T_(1,p)/Δ T_(1,p)ϕ(T_(1,p),X^(1,p)_T_(1,p))  ]
    
    by 
    _T_(1),X_T_(1)[ (σ_0^⊤)^-1Ŵ^(1,p)_Δ T_(1,p)/Δ T_(1,p)1/2(ϕ(T_(1,p),X^(1,p)_T_(1,p))  -ϕ(T_(1,p),X^(1,p^1)_T_(1,p)) )  ].
    
    
    Notice that with this version the variance of the gradient term is finite with the same argument as in the original re-normalization version in subsection <ref>.
    
    By backward induction we get the re-normalized antithetic estimator  modifying (<ref>) by:
    _k
        := 
    1/ρ(Δ T_k) ( h(T_k, X^k_T_k)+c(T_k, X^k_T_k)
    ∏_k̃∈ S(k)1/2(_k̃ -  _k̃^1_{θ(k̃) ≠ 0} + _k̃^1_{θ(k̃) = 0}) _k̃),
    	
        k ∈_T ∖_T.
    
    where the weights  are given by equation (<ref>).
    Then  we have 
    
    u(0,x)
     = _0,x[  _(1)].
    
      §.§.§ Numerical result for semi linear with re-normalization
    
    
    We apply our nested algorithm on the original re-normalized technique  and on the re-normalization technique with antithetic variables on  two test cases.
    
    First we give some results for test case A in dimension 4.
    We give the Monte Carlo error  obtained by the nested method on figure <ref>.
    For the maturity T=3, without nesting the error of the original re-normalization technique  has an order of magnitude of 2000 so the curve has not been given.
    For the maturity T=4, the nested  original re-normalization technique with  an order 2 doesn't seem to converge.
    
    
    As the maturity increases, nesting with a higher order becomes necessary.
    Notice that with the re-normalization it is possible to use the nested method of a high order because of the small number of branching used.
    For example, for T=2,  for an accuracy of 0.0004, in dimension d=4:
    
      *  the original method in section <ref> with a nested method of order 2 achieves an accuracy of 0.0004 for a CPU time of 1500 seconds using 28 cores,
    
      *  the re-normalized version of section <ref> with a nested method of order 4 reaches the same accuracy in 1800 seconds,
    
      *  the re-normalized version with antithetic of section <ref> without nesting  reaches the same accuracy in 11 seconds.
    
    
    For the same test case A we plot in dimension 6 the error on figure <ref> to show that the method converges in high dimension.
    
    
    Besides on figure <ref>, we show that the derivative is accurately calculated.
    
    We then use a second test case B :
    For a given dimension d , we take σ_0 = 1/√(d)_d, μ=,
    
     f(t,x,y,z)=  0.1/d  (· z)^2
    
     with a terminal condition g(x) = cos( x_1+ ⋯ + x_d).
    This test case cannot be solve by the nested method without re-normalization due to the high cost involved by the potential high number of branching. We give the results obtained for case B by the re-normalization methods of section  <ref> and <ref>  in dimension 4 on figure <ref>.
    
    
    At last we give the results obtained in dimension 6 pour T=1.5 and T=3 on figure <ref>.
    
    
    The nested method with re-normalization and antithetic appears to be the most effective and permits to solve semi-linear equations with quite long maturities.
    The re-normalization technique is however far more memory consuming than the original scheme of section <ref>. This memory cost  explodes with very high maturities.
    The nested version of the   original scheme of section <ref> isn't affected by these  memory problems but is affected by an explosion of the computational time with longer maturities.
    
    
     §.§ Extension to variable coefficients
    
    In the case of time and space dependent  coefficients μ and σ_0 of the PDE, it is possible to use the method consisting in “freezing” the coefficients first proposed  in <cit.> for non fixed μ  and extended in  the general case in <cit.>.
    This method increases the variance of the estimator, therefore it is more efficient for treating log maturities to use an Euler scheme to take into account the variation of the coefficients.
    Introducing an Euler time step δ t, between the dates T_k- and  T_k, the SDE is discretized as :
    
        X^k_T_k-+i δ t =    X^k_T_k-+(i-1) δ t + μ(T_k-+(i-1) δ t , X^k_T_k-+(i-1) δ t) δ t + 
           σ_0(T_k-+(i-1) δ t, X^k_T_k-+(i-1) δ t) Ŵ^k,i_δ t,
         i=1, ..,N,
        
        X^k_T_k =    X^k_T_k-+ N δ t+μ(T_k-+N δ t, X^k_T_k-+ N δ t) (Δ T_k - Nδ t) +  
           σ_0(T_k-+N δ t, X^k_T_k-+N δ t) Ŵ^k,i_Δ T_k - Nδ t,
    
    
    where N = ⌊Δ T_k/δ t⌋,  and 
    (Ŵ^k,i)_k = (k_1, ⋯, k_n-1, k_n) ∈^n, n>1, i≥ 1  is a sequence of independent d-dimensional Brownian motion.
    
    Using an integration by part on the first time step, in the original scheme of section <ref>, the gradient term in equation (<ref>) is replaced
    
        _T_(1),X_T_(1)[ (σ_0(T_(1),X^(1)_T_(1))^⊤)^-1Ŵ^(1,p),1_min(δ t,Δ T_(1,p))/min(δ t,Δ T_(1,p))ϕ(T_(1,p),X^(1,p)_T_(1,p))]
    
    In the case of the renormalization technique of section <ref>, the ghost is obtained from the original particule by removing the part associated to the first brownian.
    Then for every k = (k_1, ⋯, k_n) ∈^n_T, the particule dynamic is given by  
    
        X^k_T_k-+  δ t : =    X^k_T_k- + μ(T_k-, X^k_T_k-) δ t + 
           _k_n ∈σ_0(T_k-, X^k_T_k-) Ŵ^k,1_δ t,
        
        X^k_T_k-+i δ t =    X^k_T_k-+(i-1) δ t + μ(T_k-+(i-1) δ t , X^k_T_k-+(i-1) δ t) δ t + 
           σ_0(T_k-+(i-1) δ t, X^k_T_k-+(i-1) δ t) Ŵ^k,i_δ t,
         i=2, ..,N,
        
        X^k_T_k =    X^k_T_k-+ N δ t+μ(T_k-+N δ t, X^k_T_k-+ N δ t) (Δ T_k - Nδ t) +  
           σ_0(T_k-+N δ t, X^k_T_k-+N δ t) Ŵ^k,i_Δ T_k - Nδ t,
    
    if N>0 and
    
        X^k_T_k: =    X^k_T_k- + μ(T_k-, X^k_T_k-) Δ T_k + 
           _k_n ∈σ_0(T_k-, X^k_T_k-) Ŵ^k,1_Δ T_k
    
    otherwise. 
    
    The renormalization technique of section <ref> leads to the following estimation of the gradient  in equation (<ref>):
    _T_(1),X_T_(1)[ (σ_0(T_(1),X^(1)_T_(1))^⊤)^-1Ŵ^(1,p),1_min( δ t,Δ T_(1,p)/min( δ t,Δ T_(1,p))(ϕ(T_(1,p),X^(1,p)_T_(1,p))  -ϕ(T_(1,p),X^(1,p^1)_T_(1,p)) )  ].
    
    The renormalization technique described for the renormalization technique of section <ref> can be straightforwardly adapted  to the renormalization scheme with antithetics of 
    section <ref>.
    
    Of course using equation (<ref>) we expect that variance of the scheme will degrade with the diminution of the time step and we expect the scheme
    (<ref>) to correct this behaviour.
    On figure <ref> we give the error estimations given by the original scheme and the renormalization technique  (with anithetics of section <ref>) depending on the time step for a case with
     burgers non linearity in dimension 4 with 1e6 particles: as we refine the time step the scheme (<ref>) becomes unusable while the scheme (<ref>)
    gives stable results.
    
    
    § THE FULL NON LINEAR CASE
    
    In order to  treat some full non linear case, so with a second order derivative D^2u in f, the re-normalization technique is necessary as no distribution can meet the finite variance requirement
     even when f is linear in D^2u (see <cit.>).
    
    Suppose that the f function is as follows :
    
    f(t,x,y,z,γ)
     :=
    h(t,x)+
     c(t,x) y^ℓ_0∏_i=1^m ( (b_i · z)^ℓ^1_i)
    ∏_i=m+1^2m( (a_i : γ)^ℓ_i),
    
    for a given (ℓ_0, ℓ_1, ·, ℓ_m, ℓ_m+1, ⋯, ℓ_2m) ∈^1+2m,
    m ≥ 1, b_i:[0,T] ^d →^d for i=1, ⋯, m are  bounded continuous,
     h: [0,T] ^d →  is a  bounded continuous function, and a_i : [0,T] ^d →^d, for   i=m+1, ⋯, 2m are  bounded continuous functions.
    We note L= ∑_i=0^2mℓ_i.
    
    We use a similar algorithm to the one proposed in section <ref>. 
    Instead of approximating f using representation (<ref>), we have to take into account  the D^2u term :
    
    [c u^ℓ_0   ∏_i=1^m (b_i · Du)^ℓ^1_i∏_i=m+1^2m (a_i :D^2u)^ℓ_i ](T_(1),X_T_(1))  = 
       c ∏_j=1^ℓ_0_T_(1),X_T_(1)[ ϕ(T_(1,j),X^(1)_T_(1,j))]  
       ∏_i=1^m ( b_i(T_(1),X_T_(1)).D _T_(1),X_T_(1)[ϕ(T_(1,p),X^(1,p)_T_(1,p))])^ℓ_p^1
       ∏_i=m+1^2m (a_i :D^2 _T_(1),X_T_(1)[ϕ(T_(1,p),X^(1,p)_T_(1,p))])^ℓ_i.
    
    The terms 
        _T_(1),X_T_(1)[ ϕ(T_(1,j),X^(1)_T_(1,j))]
     and 
        ( b_i(T_(1),X_T_(1)).D _T_(1),X_T_(1)[ϕ(T_(1,p),X^(1,p)_T_(1,p))])
     are approximated by the different schemes previously seen.
    It remains to give  an approximation of the (a_i :D^2 _T_(1),X_T_(1)[ϕ(T_(1,p),X^(1,p)_T_(1,p))] u) term.
    
    
    
     §.§ Ghost particles of dimension q
    
    We extend the definition given in <cit.> of ghost tree for the full non linear case.
    For a particle in dimension (1) of generation n=1, we introduce  q associated ghost particles denoted (1^i) for i=1,...,q.
    Let ^1_T := {(1), (1^1), ...,  (1^q) }
    Then given the collection ^n_T of all particles and ghost particles of generation n, 
    we define ^n+1_T as follows.
    Given k=(k_1, ⋯, k_n) ∈^n_T, we denote by o(k) its original particle; and when k_n ∈, we denote k^i := (k_1, ⋯, k_n-1, k_n^i)  for i ∈ [1,q] and i is noted the order of k^i. The function κ allows us to give the order of a particle  for k=(k_1, ⋯, k_n) ∈^n_T :
    κ(k)     =  i,   k_n = p^i  p ∈, 
    κ(k)     =   0,  k_n = p  p ∈,
    
    The variables T_k as well as the mark θ_k inherits that of the original particle o(k). Similarly Δ T_k = Δ T_o(k).
    Denote also ^n_T := { k ∈^n_T  : o(k) ∈^n_T }.
    For every k = (k_1, ⋯, k_n) ∈^n_T ∖^n_T, 
    we define the collection of its offspring particles by
    
        h(k) := {(k_1, ⋯, k_i, 1), ⋯, (k_1, ⋯, k_i, L) },
    
    and generalizing the definition in section <ref>, we introduce q collections of all offspring ghost particles:
    
        S^i(k)
          :=  { (k_1, ⋯, k_n, 1^i), ⋯, (k_1, ⋯, k_n, L^i) },   i=1,...,q
     
    Then the collection ^n+1_T of all particles and ghost particles of generation n+1 is given by
    
        ^n+1_T
         := ∪_k ∈^n_T ∖^n_T ( S(k) ∪ S^1(k) ∪ ... ∪  S^q(k)  ).
    
    
     §.§ D^2u approximations
    
    In this section, we give some different schemes that can be used  to approximate the D^2u term and that we will compared on some numerical test cases.
    
    
    
      §.§.§ The original D^2u approximation
    
    The approximation developed in this paragraph was first proposed  in <cit.> and uses some ghost particle of dimension q=2.
    To obtained the position of a particle, we freeze its position if its order is 2 and  inverse its increment if its order is 1 , so  for every k = (k_1, ⋯, k_n) ∈^n_T
    W^k_s
     := 
    W^k-_T_k- + _κ(k)=0Ŵ^o(k)_s - T_k- - _κ(k)=1Ŵ^o(k)_s - T_k-, 
         
    		X^k_s := μ s +σ_0  W^k_s,
    		   ∀ s ∈ [T_k-, T_k].
    
    
    Then we use the following representation for the D^2u term in equation (<ref>) :
    
     D^2     _T_(1),X_T_(1)[ϕ(T_(1,p),X^(1,p)_T_(1,p))]  = 
       _T_(1),X_T_(1)[ (σ_0^⊤)^-1Ŵ^(1,p)_Δ T_(1,p)(Ŵ^(1,p)_Δ T_(1,p))^⊤ - Δ T_(1,p) I_d/(Δ T_(1,p))^2σ_0^-1ψ],  
    
    where
    ψ= 1/2[ ϕ(T_(1,p),X^(1,p)_T_(1,p)) + ϕ(T_(1,p),X^(1,p^1)_T_(1,p)) - 2 ϕ(T_(1,p),X^(1,p^2)_T_(1,p)) ].
    
    Using for example the equation (<ref>) for the first derivative Du, <cit.> gave the following re-normalized estimator defined by a backward induction:
    let _k := g(X^k_T)/(Δ T_k) for every k ∈_T, then let
    _k
        := 1/ρ(Δ T_k) ( h(T_k, X^k_T_k)+c(T_k, X^k_T_k)
    ∏_k̃∈ S(k)(_k̃_θ(k̃)=0 + (_k̃ - _k̃^2) _ 1 ≤θ(k̃) ≤ m  +  
       1/2(_k̃+_k̃^1 - _k̃^2) _m+1 ≤θ(k̃) ≤ 2m}) _k̃),
    		     k ∈_T ∖_T.
    
    where 
    _k 
       := 
    _{θ_k = 0} 
    + _{θ_k  ∈{1, ⋯, m}}b_θ_k(T_k-, X^k_T_k-) · (σ_0^⊤)^-1Ŵ^o(k)_Δ T_k/Δ T_k
       + _{θ_k ∈{m+1, ⋯, 2m}} a_θ_k : (σ_0^⊤)^-1Ŵ^o(k)_Δ T_kŴ^o(k)_Δ T_k - Δ T_k I_d/(Δ T_k)^2σ_0^-1.
    
    Then we have 
    
    u(0,x)
    = _0,x[  _(1)].
    
    
     §.§ A second representation
    
    
    This second representation  uses some ghost particle of dimension q=3.
    Let 
        (Ŵ^k,i)_k = (k_1, ⋯, k_n-1, k_n) ∈^n, n>1, i=1, 2
     be a sequence of independent d-dimensional Brownian motion, which is also independent of (Δ T_k)_k = (k_1, ⋯, k_n-1, k_n) ∈^n, n>1.
    The dynamic of the original particles and the ghosts is given by :
    
    W^k_s
     := 
    W^k-_T_k- + _κ(k)=0Ŵ^o(k),1_s - T_k- +Ŵ^o(k),2_s - T_k-/√(2) + _κ(k)=1Ŵ^o(k),1_s - T_k-/√(2) + _κ(k)=2Ŵ^o(k),2_s - T_k-/√(2)
         
    		X^k_s := μ s +σ_0  W^k_s,
    		   ∀ s ∈ [T_k-, T_k].
    
    
    
    
    We then replace (<ref>) by
    
     D^2  _T_(1),X_T_(1)[ϕ(T_(1,p),X^(1,p)_T_(1,p))]  = _T_(1),X_T_(1)[ 2 (σ_0^⊤)^-1Ŵ^(1,p),1_Δ T_k(Ŵ^(1,p),2_Δ T_k)^⊤/(Δ T_(1,p))^2σ_0^-1ψ) ],  
    
    where
    ψ = ϕ(T_(1,p),X^(1,p)_T_(1,p)) + ϕ(T_(1,p),X^(1,p^3)_T_(1,p)) -  ϕ(T_(1,p),X^(1,p^1)_T_(1,p))-  ϕ(T_(1,p),X^(1,p^2)_T_(1,p)).
    
    This scheme can be can be easily obtained by applying the differentiation rule used for semi linear equations  on two successive steps with size Δ T_(1,p)/2.
    A simple calculation  shows that the original scheme has a variance bounded by |D^2u|_∞^2 39/2 while this one has  variance bounded by |D^2u|_∞^2 9 so we expect a diminution of the variance observed with this new  scheme.
    
    This derivation on two consecutive time steps has already been used implicitly for example in <cit.> and already was  numerically superior to a scheme directly using second order Malliavin weight.
    
    Recursively the re-normalized estimator is defined by a backward induction:
    let _k := g(X^k_T)/(Δ T_k) for every k ∈_T, then let
    _k
        := 
    1/ρ(Δ T_k) ( h(T_k, X^k_T_k)+c(T_k, X^k_T_k)
    ∏_k̃∈ S(k)(_k̃_θ(k̃)=0 + (_k̃ - _k̃^3) _ 1 ≤θ(k̃) ≤ m  + 
            (_k̃+_k̃^3 - _k̃^1 - _k̃^2) _m+1 ≤θ(k̃) ≤ 2m}) _k̃),
    		     k ∈_T ∖_T.
    
    where 
    _k 
       := 
    _{θ_k = 0} 
    + _{θ_k  ∈{1, ⋯, m}}b_θ_k(T_k-, X^k_T_k-) · (σ_0^⊤)^-1Ŵ^o(k),1_Δ T_k/Δ T_k
       + _{θ_k ∈{m+1, ⋯, 2m}} a_θ_k : 2  (σ_0^⊤)^-1Ŵ^o(k),1_Δ T_kŴ^o(k),2_Δ T_k/(Δ T_k)^2σ_0^-1.
    
    Then we have 
    
    u(0,x)
    = _0,x[  _(1)].
    
    
     §.§ A third representation
    
    This representation is only the antithetic version  of the second one and uses some ghost particle of dimension q=6.
    The dynamic of the original particles and the ghosts is given by :
    
    W^k_s
        :=
    W^k-_T_k- + _κ(k)=0Ŵ^o(k),1_s - T_k- +Ŵ^o(k),2_s - T_k-/√(2) + _κ(k)=1Ŵ^o(k),1_s - T_k-/√(2) + _κ(k)=2Ŵ^o(k),2_s - T_k-/√(2) - 
       _κ(k)=4Ŵ^o(k),1_s - T_k- +Ŵ^o(k),2_s - T_k-/√(2) - _κ(k)=5Ŵ^o(k),1_s - T_k-/√(2) - _κ(k)=6Ŵ^o(k),2_s - T_k-/√(2)
            
    		X^k_s := μ s +σ_0  W^k_s,
    		   ∀ s ∈ [T_k-, T_k].
    
    
    
    
    We then replace (<ref>) by
    
     D^2    _T_(1),X_T_(1)[ϕ(T_(1,p),X^(1,p)_T_(1,p))]  = 
       _T_(1),X_T_(1)[  (σ_0^⊤)^-1Ŵ^(1,p),1(Ŵ^(1,p),2)^⊤/(Δ T_(1,p))^2σ_0^-1ψ) ],  
    
    where
    ψ     =   ϕ(T_(1,p),X^(1,p)_T_(1,p)) + 2 ϕ(T_(1,p),X^(1,p^3)_T_(1,p)) -  ϕ(T_(1,p),X^(1,p^1)_T_(1,p))-  ϕ(T_(1,p),X^(1,p^2)_T_(1,p)) + 
       ϕ(T_(1,p),X^(1,p)_T_(1,p^4))- ϕ(T_(1,p),X^(1,p^5)_T_(1,p))- ϕ(T_(1,p),X^(1,p^6)_T_(1,p)) ),
    
    and the weights are still given by equation (<ref>).
    The backward induction is defined as follows:
    let _k := g(X^k_T)/(Δ T_k) for every k ∈_T, then let
    _k
        :=
    1/ρ(Δ T_k) ( h(T_k, X^k_T_k)+ c(T_k, X^k_T_k)/2∏_k̃∈ S(k)((_k̃+_k̃^4)_θ(k̃)=0 + (_k̃ - _k̃^4) _ 1 ≤θ(k̃) ≤ m  + 
       1/2 (_k̃+ 2_k̃^3 - _k̃^1 - _k̃^2 +_k̃^4-  _k̃^5 - _k̃^6) _m+1 ≤θ(k̃) ≤ 2m}) _k̃),
    		     k ∈_T ∖_T.
    
    where the weights are given by equation (<ref>).
    And as usual we have 
    
    u(0,x)
    = _0,x[  _(1)].
    
    Extension to schemes for derivatives of order more than 3 is obvious with the two last schemes. 
    
    
     §.§ Numerical results
    
    For all test cases in this section we take μ = 0.2, σ_0 = 0.5 and we want to evaluate u(0,0.5 ).
    We test the 3 schemes previously described :
    
      *  Version 1 stands for the original version of the scheme using backward recursion (<ref>),
    
      *  Version 2 stands for the second representation using second backward recursion (<ref>),
    
      *  Version 3 stands  for the third representation corresponding to the antithetic version of the second representation and using backward recursion 
    (<ref>). Notice that in this case all terms in u in f are treated with antithetic ghosts.
     
    We give results for the non nested version as the nested version doesn't improve the results very much.
    
      *   We first choose a non linearity  
        f(u,Du,D^2u) = h(t,x)+0.1/d u(:D^2u),
     where  μ = 0.2, σ_0 = 0.5 and
    
    h(t,x)=    (α+ σ_0^2/2) cos(x_1+..+x_d) e^α(T-t)+ 
         0.1 cos(x_1+..+x_d)^2 e^2α(T-t)+ μsin(x_1+..+x_d) e^α(T-t),
    
     with α=0.2.
    We suppose that the final solution is given by g(x)= cos(x_1+..+x_d) such that the analytical solution is 
    
    u(t,x)= cos(x_1+..+x_d) e^α(T-t).
    
    This test case will be noted test C.
    In the example we want evaluate u(0,0.5 ).
     First we take d=4 and give   the results obtained for different maturities on figures <ref> and <ref>.
    
    
    We then test in dimension 6 the different schemes on figure <ref>. Besides on figure <ref> we show that
    the schemes provide a good accuracy for the computation of the derivatives by plotting (.Du) for the three versions : as expected the
    accuracy is however slightly less good than for  the function evaluation.
    
      *  At last we consider the test D where d=4,  and 
    
    f(u,Du,D^2u) = 0.0125 (.DU) (:D^2u).
    
    We give the solution and error obtained for the 3 methods on figure <ref>.
    
    On all the test cases, the last representation using antithetic variables gives the best result in term of variance reduction but at a price of memory consumption increase: as order of the ghost representation increase so does the memory needed.
     
    
    
    
    § CONCLUSION
    
    As the scheme and methods developped here let us extend the maturities than can be used to evaluate the solution of some semi linear and full non linear equation.
    This is achieved by an increase of the computational time and the memory consumption.
    
    
    plain
