Multilevel Monte Carlo algorithms for Lévy-driven SDEs with Gaussian correction

0.33 MB·English

by  Steffen Dereich

TheAnnalsofAppliedProbability
2011,Vol.21,No.1,283–311
DOI:10.1214/10-AAP695
(cid:13)c InstituteofMathematicalStatistics,2011
MULTILEVEL MONTE CARLO ALGORITHMS FOR LE´VY-DRIVEN
SDES WITH GAUSSIAN CORRECTION
1
1
By Steffen Dereich
0
2
Philipps-Universit¨at Marburg

n
a WeintroduceandanalyzemultilevelMonteCarlo algorithms for
J the computation of Ef(Y), where Y =(Y ) is the solution of
t t∈[0,1]
7 a multidimensional L´evy-driven stochastic differential equation and
f is a real-valued function on the path space. The algorithm relies

]
on approximations obtained by simulating large jumps of the L´evy
R
process individually and applying a Gaussian approximation for the
P
smalljumppart.Upperboundsareprovidedfortheworstcaseerror
.
h over the class of all measurable real functions f that are Lipschitz
t continuouswithrespecttothesupremumnorm.Theseupperbounds
a
m areeasilytractableonceoneknowsthebehavioroftheL´evymeasure
around zero.
[
Inparticular,onecanderiveupperboundsfromtheBlumenthal–

1 GetoorindexoftheL´evyprocess.InthecasewheretheBlumenthal–
v Getoor index is larger than one, this approach is superior to algo-
9 rithms that do not apply a Gaussian approximation. If the L´evy
6 processdoesnotincorporate aWienerprocessoriftheBlumenthal–
3 Getoor index β is larger than 4, then the upper bound is of order
1 3
τ−(4−β)/(6β) when the runtime τ tends to infinity. Whereas in the
.
1 case, where β is in [1,4] and the L´evy process has a Gaussian com-
0 ponent,weobtainboun3dsoforderτ−β/(6β−4).Inparticular,theerror
1
is at most of order τ−1/6.
1
:
v
i 1. Introduction. Let dY N and denote by D[0,1] the Skorokhod space
X ∈
of functions mapping [0,1] to RdY endowed with its Borel-σ-field. In this
r
a article, we analyze numerical schemes for the evaluation of
S(f):=E[f(Y)],
where
Received August 2009; revised January 2010.
AMS 2000 subject classifications. Primary 60H35; secondary 60H05, 60H10, 60J75.
Key words and phrases. Multilevel Monte Carlo, Koml´os–Major–Tusna´dy coupling,
weak approximation, numerical integration, L´evy-drivenstochastic differential equation.
This is an electronic reprint of the original article published by the
Institute of Mathematical Statistics in The Annals of Applied Probability,
2011,Vol. 21, No. 1, 283–311. This reprint differs from the original in pagination
and typographic detail.
1
2 S.DEREICH
Y =(Y ) is a solution to a multivariate stochastic differential equa-
t t [0,1]
• tion drive∈n by a multidimensional L´evy process (with state space RdY),
and
f:D[0,1] R is a Borel measurable function that is Lipschitz continuous
• →
with respect to the supremum norm.
This is a classical problem which appears for instance in finance, where Y
models the risk neutral stock price and f denotes the payoff of a (possibly
pathdependent)option,andinthepastseveralconceptshavebeenemployed
for dealing with it.
A common stochastic approach is to perform aMonte Carlo simulation of
numericalapproximations tothesolution Y.Typically,theEulerorMilstein
schemesareusedtoobtainapproximations.Alsohigherorderschemescanbe
applied provided that samples of iterated Itˆo integrals are supplied and the
coefficients of the equation are sufficiently regular. In general, the problem
is tightly related to weak approximation which is, for instance, extensively
studied in the monograph by Kloeden and Platen [12] for diffusions.
Essentially, one distinguishes between two cases. Either f(Y) depends
only on the state of Y at a fixed time or alternatively it depends on the
wholetrajectoryof Y.Intheformercase,extrapolation techniquescanoften
be applied to increase the order of convergence, see [21]. For L´evy-driven
stochasticdifferentialequations,theEulerschemewasanalyzedin[17]under
the assumption that the increments of the L´evy process are simulatable.
Approximate simulations of the L´evy increments are considered in [11].
In this article, we consider functionals f that depend on the whole tra-
jectory. Concerning results for diffusions, we refer the reader to the mono-
graph [12]. For L´evy-driven stochastic differential equations, limit theorems
in distribution are provided in [10] and [18] for the discrepancy between the
genuine solution and Euler approximations.
Recently, Giles [7, 8] (see also [9]) introduced the so-called multilevel
Monte Carlo method to compute S(f). It is very efficient when Y is a dif-
fusion. Indeed, it even can be shown that it is—in some sense—optimal,
see [5]. For L´evy-driven stochastic differential equations, multilevel Monte
Carlo algorithms are first introduced and studied in [6]. Let us explain their
findings in terms of the Blumenthal–Getoor index (BG-index) of the driv-
ing L´evy process which is an index in [0,2]. It measures the frequency of
small jumps, see (3), where a large index corresponds to a process which
has small jumps at high frequencies. In particular, all L´evy processes which
have a finite number of jumps has BG-index zero. Whenever the BG-index
is smaller or equal to one, the algorithms of [6] have worst case errors at
most of order τ 1/2, when the runtime τ tends to infinity. Unfortunately,
−
the efficiency decreases significantly for larger Blumenthal–Getoor indices.
MLMC FOR LE´VY SDEWITH GAUSSIANCORRECTION 3
Fig. 1. Order of convergence in dependence on the Blumenthal–Getoor index.
Typically, it is not feasible to simulate the increments of the L´evy process
perfectly, and one needs to work with approximations. This necessity typi-
cally worsens the performance of an algorithm, when the BG-index is larger
than one due to the higher frequency of small jumps. It represents the main
bottleneck inthesimulation. In thisarticle, weconsider approximative L´evy
increments thatsimulatethelargejumpsandapproximatethesmallonesby
a normal distribution (Gaussian approximation) in the spirit of Asmussen
and Rosin´ski [2] (see also [4]). Whenever the BG-index is larger than one,
this approach is superior to the approach taken in [6], which neglects small
jumps in the simulation of L´evy increments.
To be more precise, we establish a new estimate for the Wasserstein met-
ricbetweenanapproximativesolutionwithGaussianapproximationandthe
genuine solution, see Theorem 3.1. It is based on a consequence of Zaitsev’s
generalization [22] of the Komlo´s–Major–Tusn´ady coupling [13, 14] which
might be of its own interest itself, see Theorem 6.1. With these new esti-
mates,weanalyzeaclassofmultilevelMonteCarloalgorithmstogetherwith
acost functionwhichmeasuresthecomputationalcomplexity oftheindivid-
ualalgorithms.Weprovideuppererrorboundsforindividualalgorithmsand
optimize the error over the parameters under a given cost constraint. When
the BG-index is larger than one, appropriately adjusted algorithms lead to
significantly smaller worst case errors over the class of Lipschitz functionals
than the ones analyzed so far, see Theorem 1.1, Corollary 1.2 and Figure 1.
In particular, one always obtains numerical schemes with errors at most of
order τ 1/6 when the runtime τ of the algorithm tends to infinity.
−
4 S.DEREICH
Notation and universal assumptions. We denote by the Euclidean
|·|
norm for vectors as well as the Frobenius norm for matrices and let
k·k
denote the supremum norm over the interval [0,1]. X =(X ) denotes an
t t 0
d -dimensional L2-integrable L´evy process. By the L´evy–K≥hintchine for-
X
mula, it is characterized by a square integrable L´evy-measure ν [a Borel
measure on RdX 0 with x 2ν(dx)< ], a positive semi-definite matrix
\{ } | | ∞
ΣΣ∗ (Σ being a dX dX-matrix), and a drift b RdX via
× R ∈
Eei θ,Xt =etψ(θ),
h i
where
1
ψ(θ)= Σ θ 2+ b,θ + (ei θ,x 1 i θ,x )ν(dx).
∗ h i
2| | h i − − h i
ZRdX
Briefly, we call X a (ν,ΣΣ ,b)-L´evy process, and when b=0, a (ν,ΣΣ )-
∗ ∗
L´evy martingale. All L´evy processes under consideration are assumed to be
c`adla`g. As is well known, we can represent X as sum of three independent
processes
X =ΣW +L +bt,
t t t
where W =(W ) is a d -dimensional Wiener process and L=(L ) is
t t 0 X t t 0
a L2-martingale t≥hat comprises the compensated jumps of X. We con≥sider
the integral equation
t
(1) Y =y + a(Y )dX ,
t 0 t t
Z0 −
wherey0 RdY isafixeddeterministicinitialvalue.Weimposethestandard
∈
Lipschitz assumption on the function a:RdY RdY×dX: for a fixed K< ,
→ ∞
and all y,y′ RdY, one has
∈
a(y) a(y ) K y y and a(y ) K.
′ ′ 0
| − |≤ | − | | |≤
Furthermore, we assume without further mentioning that
x 2ν(dx) K2, Σ K and b K.
| | ≤ | |≤ | |≤
Z
We refer to the monographs [3] and [20] for details concerning L´evy pro-
cesses. Moreover, a comprehensive introduction to the stochastic calculus
for discontinuous semimartingales and, in particular, L´evy processes can be
found in [16] and [1].
In order to approximate the small jumps of the L´evy process, we need to
impose a uniform ellipticity assumption.
MLMC FOR LE´VY SDEWITH GAUSSIANCORRECTION 5
Assumption UE. There are h (0,1], ϑ 1 and a linear subspace
∈ ≥ H
of RdX such that for all h (0,h] the L´evy measure ν is supported on
B(0,h)
∈ |
and satisfies
H
1
y,x 2ν(dx) y ,x 2ν(dx) ϑ y,x 2ν(dx)
′
ϑ h i ≤ h i ≤ h i
ZB(0,h) ZB(0,h) ZB(0,h)
for all y,y with y = y .
′ ′
∈H | | | |
Mainresults. WeconsideraclassofmultilevelMonteCarloalgorithms
A
together with a cost function cost: [0, ) that are introduced explicitly
A→ ∞
in Section 2. For each algorithm S , we denote by S(f) a real-valued
∈A
random variable representing the random output of the algorithm when
applied to a given measurable funcbtion f:D[0,1] R. Wbe work in the real
→
numbermodelof computation,whichmeansthatweassumethatarithmetic
operations with real numbers and comparisons can be done in one time
unit, see also [15]. Our cost function represents the runtimeof the algorithm
reasonably well when supposing that
onecansamplefromthedistributionν /ν(B(0,h)c)andtheuniform
B(0,h)c
• |
distribution on [0,1] in constant time,
one can evaluate a at any point y RdY in constant time, and
• ∈
f can beevaluated forpiecewiseconstant functionsinlessthanaconstant
•
multiple of its breakpoints plus one time units.
As pointed out below, in that case, the average runtime to evaluate S(f) is
less than a constant multiple of cost(S). We analyze the minimal worst case
error b
b
err(τ)= inf sup E[S(f) S(f)2]1/2, τ 1.
b | − | ≥
S : f Lip(1)
∈bA ∈
cost(S) τ
≤ b
Hereandelsewhere, Lip(1) denotestheclass ofmeasurablefunctions f:D[0,
1] R that are Lipschitz continuous with respect to supremum norm with
→
coefficient one.
In this article, we use asymptotic comparisons. We write f g for 0<
liminf f limsupf < , and f -g or, equivalently g%f, for≈limsupf <
g ≤ g ∞ g
. Our main findings are summarized in the following theorem.
∞
Theorem 1.1. Assume that Assumption UE isvalid and let g:(0, )
∞ →
(0, ) be a decreasing and invertible function such that for all h>0
∞
x 2
| | 1ν(dx) g(h)
h2 ∧ ≤
Z
6 S.DEREICH
and, for a fixed γ>1,
γ
(2) g h 2g(h)
2 ≥
(cid:18) (cid:19)
for all sufficiently small h>0.
(I) If Σ=0 or
g 1(x)%x 3/4 as x ,
− −
→∞
then
err(τ)-g 1((τlogτ)2/3)τ1/6(logτ)2/3 as τ .
−
→∞
(II) If
g 1(x)-x 3/4 as x ,
− −
→∞
then
logτ
err(τ)- as τ ,
sg∗(τ) →∞
where g (τ)=inf x>1:x3g 1(x)2(logx) 1 τ .
∗ − −
{ ≥ }
The class of algorithms together with appropriate parameters which
A
establish the error estimates above are stated explicitly in Section 2.
In terms of the Blumenthal–Getoor index
(3) β:=inf p>0: x pν(dx)< [0,2]
| | ∞ ∈
(cid:26) ZB(0,1) (cid:27)
we get the following corollary.
Corollary1.2. Assume that Assumption UE isvalid and that the BG-
index satisfies β 1. If Σ=0 or β 4, then
≥ ≥ 3
4 β
sup γ 0:err(τ)-τ−γ − ,
{ ≥ }≥ 6β
and, if Σ=0 and β< 4,
6 3
β
sup γ 0:err(τ)-τ γ .
−
{ ≥ }≥ 6β 4
−
MLMC FOR LE´VY SDEWITH GAUSSIANCORRECTION 7
Visualization of the results and relationship to other work. Figure1illus-
trates our findings and related results. The x-axis and y-axis represent the
Blumenthal–Getoor index and the order of convergence, respectively. Note
that MLMC 0 stands for the multilevel Monte Carlo algorithm which does
not apply a Gaussian approximation, see [6]. Both lines marked as MLMC
1 illustrate Corollary 1.2, where the additional (G) refers to the case where
the SDE comprises a Wiener process.
These results are to be compared with the results of Jacod et al. [11].
Here an approximate Euler method is analyzed by means of weak approx-
imation. In contrast to our investigation, the object of that article is to
compute Ef(X ) for a fixed time T >0. Under quite strong assumptions
T
(for instance, a and f have to be four times continuously differentiable and
theeights momentof theL´evy processneedstobefinite),they provideerror
bounds for a numerical scheme which is based on Monte Carlo simulation
of one approximative solution. In the figure, the two lines quoted as JKMP
represent the order of convergence for general, respectively pseudo symmet-
rical, L´evy processes. Additionally to the illustrated schemes, [11] provide
an expansion which admits a Romberg extrapolation under additional as-
sumptions.
Westressthefactthatouranalysisisapplicabletogeneralpathdependent
functionals and that our error criterion is the worst case error over the
class of Lipschitz continuous functionals with respect to supremum norm.
In particular, our class contains most of the continuous payoffs appearing in
finance.
We remark that our results provide upper bounds for the inferred error
and sofar nolower boundsareknown.Theworst exponentappearingin our
estimates is 1 which we obtain for L´evy processes with Blumenthal–Getoor
6
index 2. Interestingly, this is also the worst exponent appearing in [19] in
the context of strong approximation of SDEs driven by subordinated L´evy
processes.
Agenda. The article is organized as follows. In Section 2, we introduce
a class of multilevel Monte Carlo algorithms together with a cost function.
Here, we also provide the crucial estimate for the mean squared error which
motivates theconsideration oftheWasserstein distancebetween an approxi-
mativeandthegenuinesolution,see(6).Section3statesthecentralestimate
for the former Wasserstein distance, see Theorem 3.1. In this section, we ex-
plain the strategy of the proof and the structure of the remaining article
in detail. For the proof, we couple the driving L´evy process with a L´evy
process constituted by the large jumps plus a Gaussian compensation of the
small jumps and we write the difference between the approximative and the
genuine solution as a telescoping sum including further auxiliary processes,
see (9) and (10). The individual errors are then controlled in Sections 4 and
8 S.DEREICH
5 for the terms which donot dependon the particular choice of the coupling
and in Section 7 for the error terms that do dependon the particular choice.
In between, in Section 6, we establish the crucial KMT like coupling result
for the L´evy process. Finally, in Section 8, we combine the approximation
result for the Wasserstein metric (Theorem 3.1) with estimates for strong
approximation of stochastic differential equations from [6] toprove themain
results stated above.
2. Multilevel Monte Carlo. Based on a number of parameters, we define
a multilevel Monte Carlo algorithm S: We denote by m and n ,...,n nat-
1 m
ural numbers and let ε ,...,ε and h ,...,h denote decreasing sequences
1 m 1 m
of positive reals. Formally, the algorbithm S can be represented as a tuple
constituted by these parameters, and we denote by the set of all possi-
A
ble choices for S. We continue with defininbg processes that depend on the
latter parameters. For ease of notation, the parameters are omitted in the
definitions belowb.
WechooseasquarematrixΣ(m) suchthat(Σ(m)(Σ(m)) ) = x x
∗ i,j B(0,hm) i j×
ν(dx).Moreover,fork=1,...,m,weletL(k)=(L(k)) denotethe(ν ,
t t≥0 R |B(0,hk)c
0)-L´evy martingale which comprises the compensated jumps of L that are
larger than h , that is
k
(k)
L = ∆L t xν(dx).
t s t1{|∆Ls|≥hk} s− ZB(0,hk)c
X≤
Here and elsewhere, we denote ∆L =L L . We let B =(B ) be an
t t t t t 0
independent Wiener process (independent−of W−and L(k)), and con≥sider, for
k=1,...,m, the processes (k)=(ΣW +Σ(m)B +L(k)+bt) as driving
processes. Let Υ(k) denote tXhe solutiontto t t t≥0
t
(k) (k)
Υ =y + a(Υ )d ,
t 0 Z0 s− Xι(k)(s)
where (ι(k)(t)) is given via ι(k)(t)=max(I(k) [0,t]) and the set I(k) is
t 0
≥ (k) ∩
constituted by the random times (Tj )j Z+ that are inductively defined via
∈
(k)
T =0 and
0
(k) (k) (k)
T =inf t (T , ): ∆L h or t=T +ε .
j+1 { ∈ j ∞ | t|≥ k j k}
Clearly, Υ(k) is constant on each interval [T(k),T(k)) and one has
j j+1
(k) (k) (k)
(4) Υ =Υ +a(Υ )( ).
Tj(+k)1 Tj(k) Tj(k) XTj(+k)1 −XTj(k)
MLMC FOR LE´VY SDEWITH GAUSSIANCORRECTION 9
Note that we can write
m
E[f(Υ(m))]= E[f(Υ(k)) f(Υ(k 1))]+E[f(Υ(1))].
−
−
k=2
X
The multilevel Monte Carlo algorithm—identified with S—estimates each
expectation E[f(Υ(k)) f(Υ(k 1))] (resp., E[f(Υ(1))]) individually by sam-
−
−
pling independently nk (resp., n1) versions of f(Υ(k)) fb(Υ(k−1)) [f(Υ(1))]
−
and taking the average. The output of the algorithm is then the sum of the
individual estimates. We denote by S(f) a random variable that models the
random output of the algorithm when applied to f.
b
The mean squared error of an algorithm. The Monte Carlo algorithm
introduced above induces the mean squared error
m
1
mse(S,f)= E[f(Y)] E[f(Υ(m))] 2+ var(f(Υ(k)) f(Υ(k 1)))
−
| − | n −
k
k=2
X
b 1
+ var(f(Υ(1))),
n
1
whenappliedtof.FortwoD[0,1]-valuedrandomelementsZ(1) andZ(2),we
denoteby (Z(1),Z(2))theWasserstein metricof second-orderwithrespect
W
to supremum norm, that is
1/2
(5) (Z(1),Z(2))=inf z(1) z(2) 2dξ(z(1),z(2)) ,
W ξ k − k
(cid:18)Z (cid:19)
where the infimum is taken over all probability measures ξ on D[0,1]
×
D[0,1] having first marginal P and second marginal P . Clearly, the
Z(1) Z(2)
Wasserstein distance depends only on the distributions of Z(1) and Z(2).
Now, we get for f Lip(1), that
∈
m
1
mse(S,f) (Y,Υ(m))2+ E[ Υ(k) Υ(k 1) 2]
−
≤W n k − k
k
k=2
(6) X
b
1
+ E[ Υ(1) y 2].
0
n k − k
1
We set
mse(S)= sup mse(S,f),
f Lip(1)
∈
b b
and remark that estimate (6) remains valid for the worst case error mse(S).
The main task of this article is to provide good estimates for the Wasser-
stein metric (Y,Υ(m)). The remaining terms on the right-hand side of (b6)
W
are controlled with estimates from [6].
10 S.DEREICH
The cost function. In order to simulate one pair (Υ(k 1),Υ(k)), we need
−
to simulate all displacements of L of size larger or equal to h on the time
k
interval[0,1].Moreover,weneedtheincrementsoftheWienerprocessonthe
timeskeleton (I(k 1) I(k)) [0,1].Thenwecanconstructourapproximation
−
∪ ∩
via (4). In the real number model of computation (under the assumptions
described in the Introduction), this can be performed with runtime less
than a multiple of the number of entries in I(k) [0,1], see [6] for a detailed
∩
description of an implementation of a similar scheme. Since
1 1
E[#(I(k) [0,1])] 1+ +E =ν(B(0,h )c)+ +1,
∩ ≤ εk (cid:20)t [0,1]1{|∆Lt|≥hk}(cid:21) k εk
∈X
we define, for S ,
∈A
m
b cost(S)= n ν(B(0,h )c)+ 1 +1 .
k k
ε
k=1 (cid:20) k (cid:21)
X
b
Then supposing that ε 1 and ν(B(0,h )c) 1 for k =1,...,m, yields
1 ≤ k ≤ εk
that
m
1
(7) cost(S) 3 n .
k
≤ ε
k
k=1
X
b
Algorithms achieving the error rates of Theorem 1.1. Let us now quote
the choice of parameters which establish the error rates of Theorem 1.1. In
general, one chooses ε =2 k and h =g 1(2k) for k Z . Moreover, in
k − k − +
∈
case (I), for sufficiently large τ, one picks
g 1(2k)
m= log C (τlogτ)2/3 and n = C τ1/3(logτ) 2/3 −
⌊ 2 1 ⌋ k 2 − g 1(2m)
(cid:22) − (cid:23)
for k=1,...,m,
whereC and C are appropriateconstants that donotdependon τ.Incase
1 2
(II), one chooses
g (τ)2 g 1(2k)
∗ −
m= log C g (τ) and n = C
⌊ 2 1 ∗ ⌋ k 2logg (τ)g 1(2m)
(cid:22) ∗ − (cid:23)
for k=1,...,m,
where again C and C are appropriate constants. We refer the reader to
1 2
the proof of Theorem 1.1 for the error estimates of this choice.