Table Of ContentVariations on branching methods for non linear PDEs
Xavier Warin ∗
January 27, 2017
7
1
0
2
n Abstract
a
J The branching methods developed in [9], [11] are effective methods to solve some
6
semi linear PDEs and are shown numerically to be able to solve some full non linear
2
PDEs. ThesemethodsarehoweverrestrictedtosomesmallcoefficientsinthePDEand
]
small maturities. This article shows numerically that these methods can be adapted
R
P to solve the problems with longer maturities in the semi-linear case by using a new
h. derivation scheme and some nested method. As for the case of full non linear PDEs,
t we introduce new schemes and we show numerically that they provide an effective
a
m alternative to the schemes previously developed.
[
1
1 Introduction
v
0
6
The resolution of low dimensional non linear PDEs is often achieved by some deter-
6
7 ministic methods such as finite difference schemes, finite elements and finite volume.
0
Due the curse of dimensionality, these methods cannot be used in dimension greater
.
1
than three : both the computer time and the memory required are too large even for
0
7 supercomputers. In the recent years the probabilistic community has developed some
1
representation of semi linear PDE:
:
v
i
X −∂ u−Lu = f(u,Du), u = g, t < T, x ∈ Rd, (1.1)
t T
r
a
by means of backward stochastic differential equations (BSDE), as introduced by [13].
NumericalMonteCarloalgorithmshavebeendevelopedtosolveefficientlytheseBSDE
by [2], [14]. The representation of the following full non linear PDE:
−∂ u−Lu = f(u,Du,D2u), u = g, t < T, x ∈ Rd, (1.2)
t T
has been given by the mean of second order backward stochastic differential equation
(SOBSDE) by [4]. A numerical algorithm developed by [6] 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
∗EDF R&D & FiME, Laboratoire de Finance des March´es de l’Energie, [email protected]
1
expectation and the most effective implementation is based on regression methods as
developed in [8], [12]. These regression methods develop an approximation of condi-
tional expectations based on an expansion on basis functions. The size of this expan-
sion 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
indimension4or5asregressionshasbeensuccessfullyusedindimension6in[3]using
some local regression function.
Recently a new representation of semi linear equations (1.1) for a polynomial function
f of u and Du has been given by [9] : this representation uses the automatic differenti-
ation approximation as used in [7], [1], [10], and [5]. 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 (1.2).
[11] have introduced a re-normalization technique improving numerically the conver-
gence of the scheme diminishing the variance observed for the semi linear case. Be-
sides,theauthorshavedintroducedaschemetosolvethefullnonlinearequation(1.2).
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 devel-
oped in [9, 11]. 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 [11] 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 when the numerical variance is bounded.
2
2 The Semi Linear case
Let σ ∈ Sd be some constant non-degenerate matrix, µ ∈ Rd be some constant vector
0
f : [0,T]×Rd×R×Rd → R, and g : Rd → R bounded Lipschitz functions, we consider
the semi linear parabolic PDE:
1
∂ u+ σ σ(cid:62) : D2u+µ.Du+f(·,u,Du) = 0, on [0,T)×Rd, (2.1)
t 2 0 0
with terminal condition u(T,·) = g(·) where A : B := Trace(AB(cid:62)) for two matrices
A, B ∈ Md.
2
When f is a polynomial in (u,Du) in the form
m
(cid:88) (cid:89)
f(t,x,y,z,γ) = c (t,x)y(cid:96)0 (b ·z)(cid:96)i,
(cid:96) i
(cid:96)=((cid:96)0,(cid:96)1,·,(cid:96)m)∈L i=1
for some m ≥ 1, L ⊂ N1+m, where (b ) is a sequence of Rd−valued bounded
i i=1,m
continuous functions defined on [0,T] × Rd, and (c ) is a sequence of bounded
(cid:96) (cid:96)∈L
continuous functions defined on [0,T]×Rd. [9] 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 .Du) , i = 1,m :
i
m
(cid:89)
f(t,x,y,z) =h(t,x)+c(t,x)y(cid:96)0 (b ·z)(cid:96)i, (2.2)
i
i=1
for some m ≥ 1, where (b ) is a sequence of Rd−valued bounded continuous
i i=1,m
(cid:88)
function defined on [0,T]×Rd, ((cid:96) ) ∈ Nm+1 supposing that (cid:96) > 0, and c is
i i=0,m i
i=0,m
a bounded continuous function defined on [0,T]×Rd. We note L = (cid:80)m (cid:96) .
i=0 i
Remark 2.1. The case with f a general polynomial only complexifies the notation :
it can be simply treated as in [9] by introducing some probability mass function (p )
(cid:96) (cid:96)∈L
(cid:80)
(i.e. p ≥ 0 and p = 1) that are used to select with monomial to consider
(cid:96) (cid:96)∈L (cid:96)
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 [9].
2.1 Variation on the original scheme of [9]
In this section we present the original scheme of [9] and explain how to diminish the
variance increase the maturities of the problem.
2.1.1 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 |(cid:96)| offsprings. We introduce
a sequence of i.i.d. positive random variables (τk) with all the
k=(k1,···,kn−1,kn)∈Nn,n>1
values k ∈ [1,L], for i > 0.
i
We construct an age-dependent branching process using the following procedure :
1. We start from a particle marked by 0, indexed by (1), of generation 1, whose
arrival time is given by T := τ(1)∧T.
(1)
2. Letk = (k ,··· ,k ,k ) ∈ Nn beaparticleofgenerationn,witharrivaltimeT
1 n−1 n k
that branches into L offspring particles noted (k ,··· ,k ,k ,i) for i = 1,...,L.
1 n−1 n
We define the set of its offspring particles by
S(k) := {(k ,··· ,k ,1),··· ,(k ,··· ,k ,L)},
1 n 1 n
3
We first mark the (cid:96) particles by 0, the (cid:96) next by 1 , and so on, so that each
0 1
particle has a mark i for i = 0,··· ,m.
3. For a particle k = (k ,··· ,k ,k ) of generation n+1, we denote by k− :=
1 n n+1
(k ,··· ,k ) the “parent” particle of k, and the arrival time of k is given by
1 n
T := (cid:0)T +τk(cid:1)∧T. Let us denote ∆T = T −T .
k k− k k k−
4. In particular, for a particle k = (k ,··· ,k ) of generation n, and T is its birth
1 n k−
time and also the arrival time of k−. Moreover, for the initial particle k = (1),
one has k− = ∅, and T = 0.
∅
We denote further
(cid:40)(cid:8) (cid:9)
k of generation n s.t. T ≤ t < T , when t ∈ [0,T),
θ := mark of k, Kn := k− k
k t
{k of generation n s.t. T = T}, when t = T,
k
and also
Kn := ∪ Kn, K := ∪ Kn and K := ∪ Kn.
t s≤t s t n≥1 t t n≥1 t
Clearly, K (resp. Kn) denotes the set of all living particles (resp. of generation n)
t t
n
in the system at time t, and K (resp. K ) denotes the set of all particles (resp. of
t t
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 (Wˆ k) be a sequence of independent d-
k=(k1,···,kn−1,kn)∈Nn,n>1
dimensionalBrownianmotion,whichisalsoindependentof(τk) .
k=(k1,···,kn−1,kn)∈Nn,n>1
DefineW(1) = Wˆ (1) forallt ∈ (cid:2)0,T (cid:3)andthenforeachk = (k ,··· ,k ) ∈ K \{(1)},
t t (1) 1 n T
define
Wk := Wk− +Wˆ k , for all t ∈ [T ,T ]. (2.3)
t Tk− t−Tk− k− k
Then (Wk) is a branching Brownian motion.
· k∈KT
2.1.2 The original algorithm
Let us denote F¯(t) := (cid:82)∞ρ(s)ds. Denoting Xk := x+µt+σ Wk for all k ∈ K and
t t 0 t T
t ∈ [T ,T ] and by E the expectation operator conditional on the starting data
k− k t,x
X = x at time t, we obtain from the Feynman-Kac formula the representation of the
t
solution u of equation (2.1) as:
u(0,x) = E (cid:104)F¯(T)g(XT) +(cid:90) T f(u,Du)(t,Xt)ρ(t)dt(cid:105) = E (cid:2)φ(cid:0)T ,X(1) (cid:1)(cid:3),
0,x F¯(T) ρ(t) 0,x (1) T(1)
0
(2.4)
where T := τ(1)∧T, and
(1)
φ(t,y) := 1{t≥T}g(y)+1{t<T}(h+cu(cid:96)0(cid:89)m (b ·Du)(cid:96)i)(t,y). (2.5)
F¯(T) ρ(t) i
i=1
4
On the event {1 }, using the independence of the (τk,Wk) we are left to calcu-
{T <T}
(1)
late
[cu(cid:96)0(cid:89)m (b ·Du)(cid:96)i](T ,X ) = c(cid:89)(cid:96)0 E (cid:2)φ(cid:0)T ,X(1) (cid:1)(cid:3)
i (1) T(1) T(1),XT(1) (1,j) T(1,j)
i=1 j=1
m
(cid:89)(b (T ,X ).DE (cid:2)φ(cid:0)T ,X(1,p) (cid:1)(cid:3))(cid:96)i (2.6)
i (1) T(1) T(1),XT(1) (1,p) T(1,p)
i=1
Using differentiation with respect to the heat kernel, i.e. the marginal density of the
Brownian motion we get :
[cu(cid:96)0(cid:89)m (b ·Du)(cid:96)i](T ,X ) = c(cid:89)(cid:96)0 E (cid:2)φ(cid:0)T ,X(1,j) (cid:1)(cid:3)
i (1) T(1) T(1),XT(1) (1,j) T(1,j)
i=1 j=1
m Wˆ (1,p)
(cid:89)(b (T ,X ).E (cid:2)(σ(cid:62))−1 ∆T(1,p)φ(cid:0)T ,X(1,p) (cid:1)(cid:3))(cid:96)i (2.7)
i (1) T(1) T(1),XT(1) 0 ∆T(1,p) (1,p) T(1,p)
i=1
Usingequations(2.4)and(2.7)recursivelyandthetowerproperty,wegetthefollowing
representation
(cid:104) (cid:105)
u(0,x) = E0,x ψ(cid:98)(1) (2.8)
g(Xk)−g(Xk )1
where ψ(cid:98)(1) is given by the backward recursion : let ψ(cid:98)k := T F(∆TTk−) {θk(cid:54)=0} for
k
every k ∈ K , then let
T
ψ(cid:98)k := ρ(∆1T )(cid:0)h(Tk,XTkk)+c(Tk,XTkk) (cid:89)ψ(cid:98)k˜Wk˜(cid:1), for k ∈ KT \KT. (2.9)
k
k˜∈S(k)
where
b (T ,Xk )·(σ(cid:62))−1Wˆ k
W = 1 + 1 θk k− Tk− 0 ∆Tk. (2.10)
k {θk=0} {θk(cid:54)=0} ∆T
k
and we have used that E (cid:2)g(Xk )b (T ,Xk )·σ(cid:62))−1Wˆ k (cid:3) = 0. This backward
0,x Tk− θk k− Tk− 0 ∆Tk
representation is slightly different from the elegant representation introduced in [9].
Clearly on our case the variance of the method used will be lower than with the rep-
resentation in [9] for a similar computational cost.
In the case where the operator f is linear and a function of the gradient ((cid:96) = 0,
0
m = 1 and (cid:96) = 1) using the arguments in [5] it can be easily seen by conditioning
1
with respect to the number of branching that equation (2.8) is of finite variance if
1 = O(xα) as x −→ 0 with α ≥ 0.
xρ(x)2
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, [9] have shown that the variance is in fact finite for maturities
5
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
thesameforthevaluesofκbetween0.4and0.5butahigherκvaluelimitsthenumber
of branching thus meaning a smaller computational cost.
2.1.3 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 method-
ology 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 E
t,x
for t > 0 appearing in equation (2.7) only once. Using nested Monte Carlo, so by
sampling each term of equation (2.7) 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
samplingtoestimateeachfunctionuorDuateachbranching. Ofcoursethecomputa-
tional 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
resolutionanddecreasingthenumberofsamplesastimegoesbyorasthenumber
ofbranchingincreases. Themethodologyworksslightlybetterthanare-sampling
with a constant number of particles but has to be adapted to each maturity and
each case so it has been given up.
• Anotherobservationisthatthegammalawisonlynecessarytotreatthegradient
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 Id, µ = 0,
d
f(t,x,y,z) = d(t,x)+y(b·z),
where b := 0.2(1+ 1,1+ 2,··· ,2) and
d d d
(cid:16) σ2 3d+1 (cid:17)
h(t,x) := cos(x +···+x ) α+ 0 +csin(x +···+x ) eα(T−t) eα(T−t).
1 d 1 d
2 2d
With terminal condition g(x) = cos(x +···+x ), the explicit solution of semi linear
1 d
PDE (2.1) is given by
u(t,x) = cos(x +···+x )eα(T−t).
1 d
6
Our goal is to estimate u at t = 0, x = 0.51. 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/θ)1 with
Γ(κ)θκ {s>0}
Γ(κ) := (cid:82)∞sκ−1e−sds and the parameters κ = 0.5, 1 = 0.4 to calculate the Du
0 θ
terms.
We first give on figures 1, 2 and 3 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.
Nestedncurvesstandforthecurvesusingthenestedmethodofordern, sotheNested
1 curve stands for the original method.
Figure 1: Estimation and error in d = 4 on case test A. Maturity T = 1.
Figure 2: Estimation and error in d = 4 on case test A. Maturity T = 1.5
7
On figure 3, for maturity 2.5 the error observed with the orignal method (Nested
1) is around 1000 so it has not been plotted.
Figure 3: Error in d = 4 on case test A. Maturity T = 2., T = 2.5
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 4
givingananalyticalsolution−1.4769: onceagaintheoriginalmethodfailstoconverge
while the nested one give good results.
Figure 4: Estimation and error in d = 6 on case test A. Maturity T = 2.
2.2 Adaptation of the original branching to the re-normalization
technique
Asintroducedin[11],weintroduceamodificationoftheoriginalbranchingprocessthat
letususeexponentiallawsforthebranchingdatestotreattheDutermsinthemethod
1
previously described. Recall that K = {(1)}, we introduce an associated ghost parti-
T
cle, denoted by (11), and denote K(cid:101)1 := {(1),(11)}. Next, given the collection K(cid:101)n of all
T T
particles (as well as ghost particles) of generation n, we define the collection K(cid:101)n+1 as
T
8
follows. For every k = (k1,··· ,kn) ∈ K(cid:101)n, we denote by o(k) = (kˆ1,··· ,kˆn) its original
particles, where kˆ := j when k = j or j1. Further, when k = (k ,··· ,k ) is such that
i i 1 n
kn ∈ N, wedenotek1 := (k1,··· ,kn−1,kn1). Themarkofk ∈ K(cid:101)n willbethesameasits
original particle o(k), i.e. θ := θ ; and T := T , ∆T := ∆T and τk = τo(k).
k o(k) k o(k) k o(k)
Define also K(cid:98)Tn := {k ∈ K(cid:101)Tn : o(k) ∈ KT}. For every k = (k1,··· ,kn) ∈ K(cid:101)Tn \K(cid:98)Tn, we
still define the set of its offspring particles by
S(k) := {(k ,··· ,k ,1),··· ,(k ,··· ,k ,L)},
1 n 1 n
and the set of ghost offspring particles by
S1(k) := (cid:8)(k ,··· ,k ,11),··· ,(k ,··· ,k ,L1)(cid:9).
1 n 1 n
Then the collection K(cid:101)n+1 of all particles (and ghost particles) of generation n+1 is
T
K(cid:101)n+1 := ∪ (cid:0)S(k)∪S1(k)(cid:1).
T k∈K(cid:101)n\K(cid:98)n
T T
Define also
K(cid:101)T := ∪n≥1K(cid:101)Tn, and K(cid:98)T := ∪n≥1K(cid:98)Tn.
2.2.1 The original re-normalization technique
We next equip each particle with a Brownian motion in order to define a branching
Brownian motion. Further, let W0∅ := 0, and for every k = (k1,··· ,kn) ∈ K(cid:101)Tn, let
Wsk := WTkk−− + 1kn∈NWˆso−(kT)k−, and Xsk := µs+σ0Wsk, ∀s ∈ [Tk−,Tk]. (2.11)
On figure 5, we give the original Galton-Watson tree and the ghost particles asso-
ciated.
9
(b) Tree with ghost particle k = (1,11)
(a) Original Galton-Watson tree
WW((11),1=)=WˆWˆ(1()1)+Wˆ(1,1) WW((11,)1=1)W=ˆ(W1ˆ)(1)
W(1,2)=Wˆ(1)+Wˆ(1,2) W(1,2)=Wˆ(1)+Wˆ(1,2)
W(1,1,1)=Wˆ(1)+Wˆ(1,1)+Wˆ(1,1,1) W(1,11,1)=Wˆ(1)+Wˆ(1,1,1)
W(1,1,2)=Wˆ(1)+Wˆ(1,1)+Wˆ(1,1,2) W(1,11,2)=Wˆ(1)+Wˆ(1,1,2)
(d) Tree with ghost particles k = (11) and
(c) Tree with ghost particle k = (11) k = (11,11)
W(11)=0 W(11)=0
W(11,1)=Wˆ(1,1) W(11,11)=0
W(11,2)=Wˆ(1,2) W(11,2)=Wˆ(1,2)
W(11,1,1)=Wˆ(1,1)+Wˆ(1,1,1) W(11,11,1)=Wˆ(1,1,1)
W(11,1,2)=Wˆ(1,1)+Wˆ(1,1,2) W(11,11,2)=Wˆ(1,1,2)
Figure 5: Original Galton-Watson tree, different trees with ghost particles (excluding ghost
particles at the extreme leaves) for a Brownian motion where Wk stands for Wk and Wˆ k
T
k
stands for Wˆ k .
∆T
k
Theinitialequation(2.4)remainsunchanged(firststepofthealgorithm)butequa-
tion (2.7) is modified by replacing the term
Wˆ (1,p)
E (cid:2)(σ(cid:62))−1 ∆T(1,p)φ(cid:0)T ,X(1,p) (cid:1)(cid:3)
T(1),XT(1) 0 ∆T(1,p) (1,p) T(1,p)
by
Wˆ (1,p)
E (cid:2)(σ(cid:62))−1 ∆T(1,p)(cid:0)φ(cid:0)T ,X(1,p) (cid:1)−φ(cid:0)T ,X(1,p1)(cid:1)(cid:1)(cid:3). (2.12)
T(1),XT(1) 0 ∆T(1,p) (1,p) T(1,p) (1,p) T(1,p)
Notice that since W(1,p1) has been obtained by (2.11), Wˆ (1,p) and φ(cid:0)T ,W(1,p1)(cid:1)
∆T (1,p) T
(1,p) (1,p)
are orthogonal so that adding the second term acts as a control variate. Recursively
using the modified version of equation (2.7) induced by the use of (2.12), [11] gave
10
