Central limit theorem for the multilevel Monte Carlo Euler method

0.26 MB·English

by  Mohamed Ben Alaya

TheAnnalsofAppliedProbability
2015,Vol.25,No.1,211–234
DOI:10.1214/13-AAP993
(cid:13)c InstituteofMathematicalStatistics,2015
CENTRAL LIMIT THEOREM FOR THE MULTILEVEL
5 MONTE CARLO EULER METHOD
1
0
By Mohamed Ben Alaya1 and Ahmed Kebaier2
2

n Universit´e Paris 13
a
J ThispaperfocusesonstudyingthemultilevelMonteCarlomethod

6 recentlyintroducedbyGiles [Oper. Res. 56(2008) 607–617] which is
2 significantly more efficient than the classical Monte Carlo one. Our
aim is to prove a central limit theorem of Lindeberg–Feller type for

] the multilevel Monte Carlo method associated with the Euler dis-
R
cretization scheme.Todoso, weprovefirstastable lawconvergence
P
theorem, in the spirit of Jacod and Protter [Ann. Probab. 26 (1998)
.
h 267–307], for theEulerschemeerror on twoconsecutivelevelsof the
t algorithm.Thisleadstoanaccuratedescriptionoftheoptimalchoice
a
m ofparametersandtoanexplicitcharacterizationofthelimitingvari-
ance in the central limit theorem of the algorithm. A complexity of
[
themultilevel Monte Carlo algorithm is carried out.

1
v
5
1. Introduction. In many applications, in particular in the pricing of fi-
6
3 nancial securities, we are interested in the effective computation by Monte
6 Carlo methods of the quantity Ef(X ), where X :=(X ) is a diffusion
T t 0≤t≤T
0
process and f a given function. The Monte Carlo Euler method consists of
.
1
two steps. First, approximate the diffusion process (X ) by the Eu-
0 t 0≤t≤T
5 ler scheme (Xn) with time step T/n. Then approximate Ef(Xn) by
t 0≤t≤T T
1 1 N f(Xn ), where f(Xn ) is a sample of N independent copies
: N i=1 T,i T,i 1≤i≤N
v of f(Xn). This approximation is affected, respectively, by a discretization
i P T
X
r
a
Received March 2013; revised September2013.
1Supportedby Laboratory of Excellence MME-DII.
URL:http://labex-mme-dii.u-cergy.fr/.
2This research benefited from the support of the chair “Risques Financiers,” Fon-
dation du Risque. Also supported by Laboratory of Excellence MME-DII. URL:
http://labex-mme-dii.u-cergy.fr/.
AMS 2000 subject classifications. 60F05, 62F12, 65C05, 60H35.
Key words and phrases. Centrallimittheorem,multilevelMonteCarlomethods,Euler
scheme, finance.
This is an electronic reprint of the original article published by the
Institute of Mathematical Statistics in The Annals of Applied Probability,
2015,Vol. 25, No. 1, 211–234. This reprint differs from the original in pagination
and typographic detail.
1
2 M. BEN ALAYAANDA. KEBAIER
error and a statistical error
N
1
ε :=E(f(Xn) f(X )) and f(Xn ) Ef(Xn).
n T − T N T,i − T
i=1
X
On one hand, Talay and Tubaro [21] prove that if f is sufficiently smooth,
thenε c/nwithcagivenconstantandinamoregeneralcontext,Kebaier
n
∼
[17] proves that the rate of convergence of the discretization error ε can
n
be 1/nα for all values of α [1/2,1] (see, e.g., Kloeden and Platen [18] for
∈
more details on discretization schemes). On the other hand, the statistical
error is controlled by the central limit theorem with order 1/√N. Further,
the optimal choice of thesample size N in the classical Monte Carlo method
mainly depends on the order of the discretization error. More precisely, it
turns out that for ε =1/nα the optimal choice of N is n2α. This leads to
n
a total complexity in the Monte Carlo method of order C =n2α+1 (see
MC
Duffie and Glynn [5] for related results). Let us recall that the complexity of
an algorithm isproportionaltothemaximum numberof basiccomputations
performed by this one. Hence, expressing this complexity in terms of the
−2−1/α
discretization error ε , we get C =ε .
n MC n
Inordertoimprovetheperformanceofthismethod,Kebaierintroduceda
two-level Monte Carlo method [17] (called the statistical Romberg method)
reducingthecomplexity C whilemaintainingtheconvergence ofthealgo-
MC
rithm. This method uses two Euler schemes with time steps T/n and T/nβ,
β (0,1) and approximates Ef(X ) by
T
∈
1 N1 1 N2
f(Xnβ)+ f(Xn ) f(Xnβ),
N T,i N T,i − T,i
1 2
i=1 i=1
X X
b
where Xnβ is a second Euler scheme with time step T/nβ and such that the
T
BrownianpathsusedforXn andXnβ hastobeindependentoftheBrownian
T T
paths ubsedto simulate Xnβ.Itturnsoutthat for agiven discretization error
T
ε =1/nα (α [1/2,1]), the optimal choice is obtained for β =1/2, N =
n 1
∈
n2α and N =n2α−(1/2)b. With this choice, the complexity of the statistical
2
Romberg method is of order C =n2α+(1/2) =ε−2−1/(2α), which is lower
SR n
than the classical complexity in the Monte Carlo method.
More recently, Giles [8] generalized the statistical Romberg method of
Kebaier[17]andproposedthemultilevel MonteCarloalgorithm,inasimilar
approachtoHeinrich’smultilevelmethodforparametricintegration[12](see
also Creutzig et al. [3], Dereich [4], Giles [7], Giles, Higham and Mao [9],
Giles and Szpruch [10], Heinrich [11], Heinrich and Sindambiwe [13] and
Hutzenthaler, Jentzen and Kloeden [14] for related results). The multilevel
THE MULTILEVEL MONTE CARLO METHOD 3
MonteCarlomethodusesinformationfromasequenceofcomputationswith
decreasing step sizes and approximates the quantity Ef(X ) by
T
Q = 1 N0 f(X1 )+ L 1 Nℓ (f(Xℓ,mℓ) f(Xℓ,mℓ−1)), m N 0,1 ,
n N T,k N T,k − T,k ∈ \{ }
0 ℓ
k=1 ℓ=1 k=1
X X X
where the fine discretization step is equal to T/n thereby L= logn. For
logm
ℓ,mℓ ℓ,mℓ−1
ℓ 1,...,L , processes (X ,X ) , k 1,...,N , are inde-
∈ { } t,k t,k 0≤t≤T ∈ { ℓ}
ℓ,mℓ ℓ,mℓ−1
pendent copies of (X ,X ) whose components denote the Eu-
t t 0≤t≤T
ler schemes with time steps m−ℓT and m−(ℓ−1)T. However, for fixed ℓ, the
ℓ,mℓ ℓ,mℓ−1
simulation of (X ) and (X ) has to bebased on the same
t 0≤t≤T t 0≤t≤T
Brownian path. Concerning the first empirical mean, processes (X1 ) ,
t,k 0≤t≤T
k 1,...,N , are independent copies of (X1) which denotes the Eu-
∈{ 0} t 0≤t≤T
ler scheme with time step T. Here, it is important to point out that all
these L+1 Monte Carlo estimators have to be based on different indepen-
dent samples. Due to the above independence assumption on the paths, the
variance of the multilevel estimator is given by
L
σ2:=Var(Q )=N−1Var(f(X1))+ N−1σ2,
n 0 T ℓ ℓ
ℓ=1
X
where σ2=Var(f(Xℓ,mℓ) f(Xℓ,mℓ−1)). Assuming that the diffusion coeffi-
ℓ T − T
cients of X and the function f are Lipschitz continuous, then it is easy to
check, using properties of the Euler scheme that
L
σ2 c N−1m−ℓ
≤ 2 ℓ
ℓ=0
X
for some positive constant c (see Proposition 1 for more details). Giles
2
[8] uses this computation in order to find the optimal choice of the multi-
level Monte Carlo parameters. More precisely, to obtain adesired root mean
squared error (RMSE), say of order 1/nα, for the multilevel estimator, Giles
[8] uses the above computation on σ2 to minimize the total complexity of
the algorithm. It turns out that the optimal choice is obtained for (see The-
orem 3.1 of [8])
logn T logn
(1) N =2c n2α +1 for ℓ 0,...,L and L= .
ℓ 2 logm mℓ ∈{ } logm
(cid:18) (cid:19)
Hence, for an error ε =1/nα, this optimal choice leads to a complexity
n
for the multilevel Monte Carlo Euler method proportional to n2α(logn)2=
ε−2(logε )2. Interesting numerical tests, comparing three methods (crude
n n
4 M. BEN ALAYAANDA. KEBAIER
Monte Carlo, statistical Romberg and the multilevel Monte Carlo), were
processed in Korn, Korn and Kroisandt [19].
In the present paper, we focus on central limit theorems for the inferred
error; a question which has not been addressed in previous research. To
do so, we use techniques adapted to this setting, based on a central limit
theorem for triangular array (see Theorem 2)together with Toeplitz lemma.
It is worth to note that our approach improves techniques developed by
Kebaier [17] in his study of the statistical Romberg method (see Remark 2
for more details). Hence, our main result is a Lindeberg–Feller central limit
theorem for the multilevel Monte Carlo Euler algorithm (see Theorem 4).
Further, this allows us to prove a Berry–Esseen-type bound on our central
limit theorem.
In order to show this central limit theorem, we first prove a stable law
convergence theorem, for the Euler scheme error on two consecutive levels
mℓ−1 and mℓ, of the type obtained in Jacod and Protter [16]. Indeed, we
prove the following functional result (see Theorem 3):
mℓ
(Xℓ,mℓ Xℓ,mℓ−1) stablyU as ℓ ,
s(m 1)T − ⇒ →∞
−
whereU isthesamelimitprocessgiven inTheorem3.2ofJacodandProtter
[16].Ourresultusesstandardtoolsdevelopedintheirpaperbutitcannotbe
deducedwithoutaspecificandlaboriousstudy.Further,theirresult,namely
mℓ
(Xℓ,mℓ X) stablyU as ℓ ,
T − ⇒ →∞
r
is neither sufficient nor appropriate to prove our Theorem 4, since the mul-
tilevel Monte CarloEulermethodinvolves theerrorprocess Xℓ,mℓ Xℓ,mℓ−1
rather than Xℓ,mℓ X. −
−
ThankstoTheorem4,weobtain aprecisedescriptionforthechoice ofthe
parameters to run the multilevel Monte Carlo Euler method. Afterward, by
a complexity analysis we obtain the optimal choice for the multilevel Monte
Carlo Euler method. It turns out that for a total error of order ε =1/nα
n
the optimal parameters are given by
(m 1)T logn
(2) N = − n2αlogn for ℓ 0,...,L and L= .
ℓ mℓlogm ∈{ } logm
This leads us to a complexity proportional to n2α(logn)2 = ε−2(logε )2
n n
which is the same order obtained by Giles [8]. By comparing relations (1)
and (2), we note that our optimal sequence of sample sizes (N ) does
ℓ 0≤ℓ≤L
not depend on any given constant, since our approach is based on proving
a central limit theorem and not on obtaining an upper bound for the vari-
anceofthealgorithm.However,somenumericaltestscomparingtheruntime
THE MULTILEVEL MONTE CARLO METHOD 5
with respecttotherootmeansquareerror,show thatweareinlinewiththe
original work of Giles [8]. Nevertheless, the major advantage of our central
limittheorem is thatitfillsthegap intheliteratureforthemultilevel Monte
Carlo Euler method and allows to construct a more accurate confidence in-
terval compared to the one obtained using Chebyshev’s inequality. All these
results are stated and proved in Section 3. The next section is devoted to
recall some useful stochastic limit theorems and to introduce our notation.
2. General framework.
2.1. Preliminaries. Let (X ) be a sequence of random variables with
n
values in a Polish space E defined on a probability space (Ω, ,P). Let
F
(Ω, ,P) be an extension of (Ω, ,P), and let X be an E-valued random
F F
variable on the extension. We say that (X ) converges in law to X stably
n
aendewreite Xn stablyX, if
⇒
E(Uh(X )) E(Uh(X))
n
→
for all h:E R bounded continuous and all bounded random variable U
→ e
on (Ω, ). This convergence is obviously stronger than convergence in law
F
that we will denote here by “ .” According to Section 2 of Jacod [15] and
⇒
Lemma 2.1 of Jacod and Protter [16], we have the following result.
Lemma 1. Let V and V be defined on (Ω, ) with values in another
n
F
metric space E′.
P
If V V,X stablyX then (V ,X ) stably(V,X).
n n n n
→ ⇒ ⇒
Conversely, if (V,X ) (V,X) and V generates the σ-field , we can re-
n
alize this limit as (V,X⇒) with X defined on an extension ofF(Ω, ,P) and
F
X stablyX.
n
⇒
Now, we recall a result on the convergence of stochastic integrals for-
mulated from Theorem 2.3 in Jacod and Protter [16]. This is a simplified
version butit is sufficient for our study.Let Xn=(Xn,i) bea sequence
1≤i≤d
of Rd-valued continuous semimartingales with the decomposition
n,i n,i n,i n,i
X =X +A +M , 0 t T,
t 0 t t ≤ ≤
where, for each n N and 1 i d, An,i is a predictable process with finite
∈ ≤ ≤
variation, null at 0 and Mn,i is a martingale null at 0.
Theorem 1. Assume that the sequence (Xn) is such that
T
Mn,i + dAn,i
h iT | s |
Z0
6 M. BEN ALAYAANDA. KEBAIER
is tight. Let Hn and H be a sequence of adapted, right-continuous and left-
hand side limited processes all defined on the same filtered probability space.
If (Hn,Xn) (H,X) then X is a semimartingale with respect to the filtra-
⇒
tiongeneratedbythelimitprocess (H,X),andwehave(Hn,Xn, HndXn)
⇒
(H,X, HdX).
R
R
We recall also the following Lindeberg–Feller central limit theorem that
will be used in the sequel (see, e.g., Theorems 7.2 and 7.3 in [1]).
Theorem 2 (Central limit theorem for triangular array). Let (kn)n∈N
be a sequence such that k as n . For each n, let X ,...,X be
n→∞ →∞ n,1 n,kn
k independent random variables with finite variance such that E(X )=0
n n,k
for all k 1,...,k . Suppose that the following conditions hold:
n
∈{ }
(A1) lim kn E X 2=σ2,σ>0.
n→∞ k=1 | n,k|
(A2) Lindeberg’s condition: for all ε > 0, lim kn E(X 2
P n→∞ k=1 | n,k| ×
1 )=0. Then
{|Xn,k|>ε} P
kn
X (0,σ2) as n .
n,k
⇒N →∞
k=1
X
Moreover, if the X have moments of order p>2, then the Lindeberg’s
n,k
condition can be obtained by the following one:
(A3) Lyapunov’s condition: lim kn E X p=0.
n→∞ k=1 | n,k|
P
2.2. The Euler scheme. Let X :=(X ) be the process with values
t 0≤t≤T
in Rd, solution to
(3) dX =b(X )dt+σ(X )dW , X =x Rd,
t t t t 0
∈
whereW =(W1,...,Wq)isaq-dimensionalBrownianmotiononsomegiven
filtered probability space =(Ω, ,( ) ,P) with ( ) is the standard
t t≥0 t t≥0
Brownian filtration, b andBσ are, rFespeFctively, Rd and RFd×q valuedfunctions.
WeconsiderthecontinuousEulerapproximationXn withstepδ=T/ngiven
by
dXn=b(X )dt+σ(X )dW , η (t)=[t/δ]δ.
t ηn(t) ηn(t) t n
It is well known that under the global Lipschitz condition
C >0, such that, b(x) b(y) + σ(x) σ(y) C y x ,
T T
( ) ∃ | − | | − |≤ | − |
Hb,σ x,y Rd,
∈
THE MULTILEVEL MONTE CARLO METHOD 7
the Euler scheme satisfies the following property (see, e.g., Bouleau and
L´epingle [2]):
p 1, sup X , sup Xn Lp and
∀ ≥ | t| | t |∈
0≤t≤T 0≤t≤T
( )
P K (T)
E sup X Xn p p , K (T)>0.
| t− t | ≤ np/2 p
0≤t≤T
h i
Note that according to Theorem 3.1 of Jacod and Protter [16], under the
weaker condition
( ) b and σ are locally Lipschitz with linear growth,
b,σ
H
we have only the uniform convergence in probability, namely the property
e
P
( ) sup X Xn 0.
P | t− t |→
0≤t≤T
Foellowing the notation of Jacod and Protter [16], we rewrite diffusion (3) as
follows:
q
j
dX =ϕ(X )dY = ϕ (X )dY ,
t t t j t t
j=0
X
where ϕ is the jth column of the matrix σ, for 1 j q, ϕ = b and
j 0
Y :=(t,W1,...,Wq)′. Then the continuous Euler app≤roxi≤mation Xn with
t t t
time step δ=T/n becomes
q
(4) dXn=ϕ(Xn )dY = ϕ (Xn )dYj, η (t)=[t/δ]δ.
t ηn(t) t j ηn(t) t n
j=0
X
3. The multilevel Monte Carlo Euler method. Let (Xmℓ) denotes
t 0≤t≤T
the Euler scheme with time step m−ℓT for ℓ 0,...,L , where L=logn/
∈{ }
logm. Noting that
L
(5) Ef(Xn)=Ef(X1)+ E(f(Xmℓ) f(Xmℓ−1)),
T T T − T
ℓ=1
X
the multilevel method is to estimate independently by the Monte Carlo
methodeachoftheexpectations ontheright-handsideoftheaboverelation.
Hence, we approximate Ef(Xn) by
T
(6) Q = 1 N0 f(X1 )+ L 1 Nℓ (f(Xℓ,mℓ) f(Xℓ,mℓ−1)).
n N T,k N T,k − T,k
0 ℓ
k=1 ℓ=1 k=1
X X X
Here, it is important to point out that all these L+1 Monte Carlo esti-
mators have to be based on different, independent samples. For each ℓ
∈
8 M. BEN ALAYAANDA. KEBAIER
ℓ,mℓ ℓ,mℓ−1
1,...,L the samples (X ,X ) are independent copies of
{ } T,k T,k 1≤k≤Nℓ
ℓ,mℓ ℓ,mℓ−1
(X ,X ) whose components denote the Euler schemes with time
T T
steps m−ℓT and m−(ℓ−1)T and simulated with the same Brownian path.
Concerning the first empirical mean, the samples (X1 ) are inde-
T,k 1≤k≤N0
pendent copies of X1. The following result gives us a first description of
T
the asymptotic behavior of the variance in the multilevel Monte Carlo Euler
method.
Proposition 1. Assume that b and σ satisfy condition ( ). For a
b,σ
H
Lipschitz continuous function f:Rd R, we have
−→
L
(7) Var(Q )=O N−1m−ℓ .
n ℓ
!
ℓ=0
X
Proof. We have
L
Var(Q )=N−1Var(f(X1))+ N−1Var(f(Xℓ,mℓ) f(Xℓ,mℓ−1))
n 0 T ℓ T − T
ℓ=1
X
N−1Var(f(X1))
≤ 0 T
L
+2 N−1(Var(f(Xmℓ) f(X ))+Var(f(Xmℓ−1) f(X )))
ℓ T − T T − T
ℓ=1
X
N−1Var(f(X1))
≤ 0 T
L
+2[f] N−1E sup Xmℓ X 2+ sup Xmℓ−1 X 2 ,
lip ℓ | t − t| | t − t|
0≤t≤T 0≤t≤T
Xℓ=1 h i
|f(u)−f(v)|
where [f] :=sup . We complete the proof by using property
lip u6=v |u−v|
( ) on the strong convergence of the Euler scheme. (cid:3)
P
Inequality(7)indicatesthedependenceofthevarianceofQ onthechoice
n
oftheparametersN ,...,N .Thisvariancecanbesmallerthanthevariance
0 L
off(Xn),sothatQ appearsasagoodcandidateforthevariancereduction.
T n
Themain resultof thissection is aLindeberg–Feller central limit theorem
(see Theorem 4 below). In order to prove this result, we need to prove first
a new stable law convergence theorem for the Euler scheme error adapted
to the setting of multilevel Monte Carlo algorithm. This is crucial and is the
aim of the following subsection.
THE MULTILEVEL MONTE CARLO METHOD 9
3.1. Stable convergence. In what follows, we prove a stable law conver-
gence theorem, for the Euler scheme error on two consecutive levels mℓ−1
and mℓ, of the type obtained in Jacod and Protter [16]. Our result in Theo-
rem 3 below is an innovative contribution on the Euler scheme error that is
different and more tricky than the original work by Jacod and Protter [16]
since it involves the error process Xℓ,mℓ Xℓ,mℓ−1 rather than Xℓ,mℓ X.
Note that the study of the error Xℓ,mℓ −Xℓ,mℓ−1 as ℓ can be red−uced
− →∞
to the study of the error Xmn Xn as n where Xmn and Xn stand
− →∞
for the Euler schemes with time steps T/(mn) and T/n constructed on the
same Brownian path.
Theorem 3. Assume that b and σ are 1 with linear growth then the
C
following result holds:
mn
For all m N 0,1 (Xmn Xn) stablyU as n ,
∈ \{ } (m 1)T − ⇒ →∞
r −
with (U ) the d-dimensional process satisfying
t 0≤t≤T
1 q t
(8) U = Z Hi,jdBij, t [0,T],
t √2 t s s ∈
i,j=1 Z0
X
where
(9) Hi,j =(Z )−1ϕ˙ ϕ¯ with ϕ˙ := ϕ (X ) and ϕ¯ :=ϕ (X ),
s s s,j s,i s,j ∇ j s s,i i s
and (Z ) is the Rd×d valued process solution of the linear equation
t 0≤t≤T
q t
Z =I + ϕ˙ dYjZ , t [0,T].
t d s,j s s ∈
j=0Z0
X
Here, ϕ is a d d matrix with ( ϕ ) is the partial derivative of ϕ with
j j ik ij
∇ × ∇
respect to the kth coordinate, and (Bij) is a standard q2-dimensional
1≤i,j≤q
Brownian motion independent of W. This process is defined on an extension
(Ω, ,( ) ,P) of the space (Ω, ,( ) ,P).
t t≥0 t t≥0
F F F F
eNeoteethat byeletting formally m tendto infinity,we recover theJacod and
Protter’s result [16].
ProofofTheorem3. ConsidertheerrorprocessUmn,n=(Umn,n) ,
t 0≤t≤T
defined by
Umn,n:=Xmn Xn, t [0,T].
t t − t ∈
10 M. BEN ALAYAANDA. KEBAIER
Combining relation (4), for both processes Xmn and Xn, together with a
Taylor expansion
q
dUmn,n= ϕ˙n (Xmn Xn )dYj,
t t,j ηmn(t)− ηn(t) t
j=0
X
whereϕ˙n isthed dmatrixwhoseithrowisthegradientofthereal-valued
t,j ×
function ϕ at a point between Xn and Xmn . Therefore, the equation
ij ηn(t) ηmn(t)
satisfied by Un can be written as
t q
Umn,n= ϕ˙n Umn,ndYj +Gmn,n,
t s,j s s t
Z0 j=0
X
with
t q t q
Gmn,n= ϕ˙n (Xn Xn )dYj ϕ˙n (Xmn Xmn )dYj.
t s,j s − ηn(s) s − s,j s − ηmn(s) s
Z0 j=0 Z0 j=0
X X
In the following, let (Zmn,n) be the Rd×d valued solution of
t 0≤t≤T
t q
Zmn,n=I + ϕ˙n dYj Zmn,n.
t d s,j s s
Z0 j=0 !
X
Theorem48,page326in[20],ensuresexistenceoftheprocess((Zmn,n)−1)
t 0≤t≤T
defined as the solution of
t q t q
(Zmn,n)−1=I + (Zmn,n)−1 (ϕ˙n )2ds (Zmn,n)−1 ϕ˙n dYj.
t d s s,j − s s,j s
Z0 j=1 Z0 j=0
X X
Thanks to Theorem 56, page 333 in the same reference [20], we get
t
Umn,n=Zmn,n (Zmn,n)−1dGmn,n
t t s s
(Z0
t q
(Zmn,n)−1 (ϕ˙n )2(Xn Xn )ds
− s s,j s − ηn(s)
Z0 j=1
X
t q
+ (Zmn,n)−1 (ϕ˙n )2(Xmn Xmn )ds .
s s,j s − ηmn(s)
Z0 j=1 )
X
Since the increments of the Euler scheme satisfy
q
Xn Xn = ϕ¯n (Yi Yi )
s − ηn(s) s,i s − ηn(s)
i=0
X