Table Of ContentOPTIMAL STRONG APPROXIMATION OF THE
ONE-DIMENSIONAL SQUARED BESSEL PROCESS
MARIO HEFTER AND ANDRE´ HERZWURM
Abstract. Weconsidertheone-dimensionalsquaredBesselprocessgivenbythe
stochastic differential equation (SDE)
(cid:112)
dX =1dt+2 X dW , X =x , t∈[0,1],
t t t 0 0
and study strong (pathwise) approximation of the solution X at the final time
point t=1. This SDE is a particular instance of a Cox-Ingersoll-Ross (CIR) pro-
cess where the boundary point zero is accessible. We consider numerical methods
thathaveaccesstovaluesofthedrivingBrownianmotionW atafinitenumberof
6 time points. We show that the polynomial convergence rate of the n-th minimal
1 errors for the class of adaptive algorithms as well as for the class of algorithms
0 that rely on equidistant grids are equal to infinity and 1/2, respectively. This
2
shows that adaption results in a tremendously improved convergence rate. As a
n by-product, we obtain that the parameters appearing in the CIR process affect
a
the convergence rate of strong approximation.
J
7
]
R
1. Introduction
P
h. In recent years, strong approximation of stochastic differential equations (SDEs)
t has intensively been studied for SDEs of the form
a
m (cid:112)
(1) dX = (a−bX )dt+σ X dW , X = x , t ≥ 0,
[ t t t t 0 0
1 withaone-dimensionalBrownianmotionW, anda,x0 ≥ 0, b ∈ R, andσ > 0. These
v SDEs are known to have a unique non-negative strong solution. Such SDEs were
5
proposed in [11] as a model for short-term interest rates. The solution is called Cox-
5
4 Ingersoll-Ross (CIR) process. Moreover, CIR processes are used as the volatility
1
process in the Heston model [18].
0
Strong approximation is of particular interest due to the multi-level Monte Carlo
.
1
technique, see [14, 15, 17]. For an optimality result of this technique applied to
0
6 quadrature problems of the form E(f(X)) with f: C([0,1]) → R, we refer to [12].
1
In mathematical finance, the functional f often represents a discounted payoff of
:
v some derivative and E(f(X)) is the corresponding price.
i
X In [1], various numerical schemes have been proposed and numerically tested
r for the SDE (1) with different choices of the corresponding parameters. These
a
numerical results indicate a convergence at a polynomial rate, which depends on the
parameters a,σ. More precisely, the empirical convergence rate is monotonically
decreasing in the quotient σ2/(2a) for all numerical schemes that have been tested.
Polynomial convergence rates for strong approximation of (1) have been proven by
[4, 13, 2, 25, 19], where either a global or the final time error w.r.t. the L -norm is
p
studied. Alltheseresultsonlyholdforsomeparameterrangewithinσ2/(2a) < 2and
share the same monotonicity, see Figure 1. For an overview of numerical schemes
Key words and phrases. Cox-Ingersoll-Rossprocess;strongapproximation;n-thminimalerror;
adaptive algorithm; reflected Brownian motion.
1
2 HEFTER AND HERZWURM
and results on strong convergence without a rate we refer to [13] and the references
therein.
1
p=1
p=2
p=4
0.8
e
t
a
r
e 0.6
c
n
e
g
r
e
v
n
o 0.4
c
0.2
0
0 0.5 1 1.5 2
<2=(2a)
Figure 1. Convergence rates of the known (best) upper bounds from
[13, 2, 25, 19] are shown for different values of p ∈ [1,∞[ with error
givenby(4). Alltheseupperboundsareobtainedforthedrift-implicit
Euler scheme. The dashed line shows the corresponding empirical
convergence rates from [1] for p = 1.
In this paper we consider the particular case of the SDE (1) with
σ2/(2a) = 2.
By rescaling with (2/σ)2, we thus may restrict ourselves to SDEs of the form
(cid:112)
(2) dX = (1−bX )dt+2 X dW , X = x , t ≥ 0,
t t t t 0 0
with x ≥ 0 and b ∈ R. Furthermore, we focus on the particular instance of (2)
0
with b = 0, i.e.,
(cid:112)
(3) dX = 1dt+2 X dW , X = x , t ≥ 0.
t t t 0 0
Its solution is called the square of a 1-dimensional Bessel process. For a detailed
study of (squared) Bessel processes we refer to [26, Chap. XI].
In the context of (strong) approximation of SDEs, the majority of numerical
methods in the literature are non-adaptive [22]. A non-adaptive algorithm uses a
fixeddiscretisationofthedrivingBrownianmotionwhereasadaptivealgorithmsmay
sequentially choose the evaluation points. The most frequently studied methods in
the class of non-adaptive algorithms are Euler or Milstein-type methods that are
based on values of the driving Brownian motion on an equidistant grid. For various
OPTIMAL STRONG APPROXIMATION OF THE 1-D SQUARED BESSEL PROCESS 3
strong approximation problems of SDEs satisfying standard assumptions, adaption
does not help up to a multiplicative constant, see, e.g., [24]. In particular, [23]
showed for strong approximation of scalar SDEs at the final time point that no
adaptive method can be better (up to a multiplicative constant) than the classical
Milstein scheme. Let us stress that these standard assumptions are not fulfilled
by (3) since the diffusion coefficient is not even locally Lipschitz continuous.
In contrast to that, the main result of this paper is that adaptive methods are far
superior to methods that are based on an equidistant grid for strong approximation
of the solution X of (3). For this, we determine the polynomial convergence rate
1
of the corresponding n-th minimal errors, which will be introduced below.
Let X be the solution of (3) at time t = 1, and let p ∈ [1,∞[. The error of an
1
approximation X(cid:98) of X is defined by
1 1
(cid:16) (cid:16)(cid:12) (cid:12)p(cid:17)(cid:17)1/p
(4) e (X(cid:98) ) = E (cid:12)X −X(cid:98) (cid:12) .
p 1 1 1
Atfirst,weconsidertheclassofmethodsthatonlyusevaluesofthedrivingBrownian
motion W on an equidistant grid with n points given by
(cid:110) (cid:111)
Ceq(n) = X(cid:98)1 = Φ(W1,W2,...,W1) : Φ: Rn → R Borel-measurable .
n n
The corresponding n-th minimal error for the approximation of X is given by
1
(cid:110) (cid:111)
(5) eeq(n) = inf e (X(cid:98) ) : X(cid:98) ∈ Ceq(n) .
p p 1 1
Roughly speaking, eeq(n) is the error of the best algorithm for the approxima-
p
tion of X w.r.t. the L -norm that only uses W ,W ,...,W . Clearly, Euler and
1 p 1 2 1
n n
Milstein-type schemes fit into this class of algorithms. In the case p = 2, the opti-
mal approximation is given by the conditional expectation of X given the σ-algebra
1
generated by W ,W ,...,W .
1 2 1
n n
The class of adaptive methods that use values of the driving Brownian motion W
at n sequentially chosen points is given by
(cid:110)
Cad(n) = X(cid:98) = Φ(W ,...,W ) : Φ: Rn → R Borel-measurable,
1 t1 tn
t ∈ [0,1],
1
t = ϕ (W ), ϕ : R → [0,1] Borel-measurable,
2 2 t1 2
.
.
.
(cid:111)
t = ϕ (W ,...,W ), ϕ : Rn−1 → [0,1] Borel-measurable .
n n t1 tn−1 n
Here, in contrast to the class Ceq, the k-th evaluation site t may depend on the
k
previous k − 1 observations of W. Moreover, considering the particular choice of
constantmappingsϕ = k/nyieldsCeq(n) ⊆ Cad(n)foralln ∈ N. Then-thminimal
k
error for the approximation of X for the class of adaptive methods is given by
1
(cid:110) (cid:111)
(6) ead(n) = inf e (X(cid:98) ) : X(cid:98) ∈ Cad(n) .
p p 1 1
We clearly have ead(n) ≤ eeq(n) for all n ∈ N.
p p
In the following we present our main results. We write a (cid:52) b for sequences of
n n
non-negative reals a and b if there exists a constant c > 0 such that a ≤ c·b for
n n n n
all n ∈ N. Moreover, we write a (cid:16) b if a (cid:52) b and b (cid:52) a .
n n n n n n
4 HEFTER AND HERZWURM
We show that the polynomial convergence rate of the n-th minimal error eeq is
p
equal to 1/2 for all p ∈ [1,∞[. More precisely, Corollary 1 yields
(7) eeq(n) (cid:16) n−1/2
p
for all p ∈ [1,∞[. Of course, the constants hidden in the “(cid:16)”-notation may depend
on p. Furthermore, the corresponding upper bound is attained by the drift-implicit
Euler scheme X(cid:98)n,imp, see Theorem 2 and (23). In the more general case of the
1
SDE (2), the drift-implicit Euler scheme is given by
(cid:114) 2
(cid:113) (cid:16) (cid:17) (cid:16)(cid:113) (cid:16) (cid:17)(cid:17)2
X(cid:98)n,imp + W −W + X(cid:98)n,imp + W −W
(8) X(cid:98)n,imp = nk k+n1 nk nk k+n1 nk
k+1 2+ b
n
n
for k = 0,...,n−1 and X(cid:98)n,imp = x . Let us mention that the drift-implicit Euler
0 0
scheme is actually proposed for the SDE (1) with parameters satisfying σ2/(2a) < 2,
see [1]. Nevertheless, it is still well defined in the limiting case σ2/(2a) = 2. Note
that the upper bound from (7) is the first strong convergence result with a positive
rate in the case σ2/(2a) = 2, cf. Figure 1.
For adaptive algorithms the situation is rather different. Corollary 2 shows that
(9) ead(n) (cid:52) n−q
p
for all p ∈ [1,∞[ and for all q ∈ [1,∞[. Hence the polynomial convergence rate of
the n-th minimal error ead is equal to infinity. More precisely, for every q ∈ [1,∞[
p
we construct an adaptive algorithm that converges (at least) at a polynomial rate q,
see Theorem 4. In fact, numerical experiments suggest an exponential decay, see
Figure 3. Moreover, such algorithms can be easily implemented on a computer with
number of operations of order n2. Combining (7) and (9) establishes our claim that
adaptive algorithms are far superior to non-adaptive algorithms that are based on
equidistant grids for strong approximation of (3). Let us stress that this is the first
result on SDEs where adaption results in an improved convergence rate compared
to methods that are based on equidistant grids.
Akeystepfortheproofsof (7)and(9)consistsofidentifyingthepathwisesolution
of (3), see Proposition 1, and link this problem to global optimization under the
Wiener measure. Let us mention that the analysis of the adaptive algorithm in
Theorem 4 heavily relies on results of [9].
Although we have shown that adaptive algorithms are far superior to methods
that are based on an equidistant grid for a particular choice of the parameters of
SDE (1), it is open whether this superiority also holds for more general parameter
constellations.
We now turn to the more general case of the SDE (2) with b ∈ R. Moreover, we
consider a stronger error criterion which is pathwise given by the supremum norm.
In this case we obtain
(cid:18) (cid:18) (cid:19)(cid:19)1/p
(10) E sup (cid:12)(cid:12)X −Xn(cid:12)(cid:12)p (cid:52) n−1/2 ·(cid:112)ln(1+n)
t t
0≤t≤1
n
for all p ∈ [1,∞[, where X denotes a projected equidistant Euler scheme, see
Remark 7. This scheme coincides with the drift-implicit Euler scheme for b = 0. Let
us stress that this error bound is the first strong convergence result with a positive
OPTIMAL STRONG APPROXIMATION OF THE 1-D SQUARED BESSEL PROCESS 5
rate in the case σ2/(2a) = 2 and arbitrary b ∈ R, cf. Figure 1. At present, we have
only shown the upper bound (10). Nevertheless, we expect this upper bound to be
sharp even for adaptive algorithms.
Let us briefly comment on some consequences of the results presented above for
strong approximation of CIR processes.
It is well-known that the parameters a, b, and σ in SDE (1) have an influence
on the behavior of its solution. For instance, the solution remains strictly positive
(the boundary point 0 is inaccessible) if and only if the so-called Feller condition
σ2/(2a) ≤ 1 is satisfied. As illustrated in Figure 1, the drift-implicit Euler scheme
converges at least with rate 1 if σ2/(2a) < min(2/(3p),1/2), see [2, 25], and so
does the corresponding n-th minimal error for methods using an equidistant grid.
Hence the quotient σ2/(2a) affects the convergence rate of the n-th minimal error for
equidistant methods since it drops down to 1/2 for σ2/(2a) = 2 and b = 0 according
to (7).
In contrast to the known upper bounds, cf. Figure 1, the convergence rate of the
drift-implicit Euler scheme for (3) does not depend on the L -norm appearing in the
p
error criterion, see (7).
Let us comment on lower bounds for strong approximation of SDEs at the final
time point based on the values of the driving Brownian motion. In [10], a two-
dimensional SDE is presented where the corresponding convergence rate is shown
to be 1/2. In contrast to rate 1 for smooth scalar SDEs [23], the difficulty in [10]
arises from the presence of L´evy areas. More recently, the existence of SDEs with
smooth coefficients has been shown where the corresponding n-th minimal error
converges arbitrarily slow to zero, see [16, 20]. It is crucial that these SDEs are
multi-dimensional. Apart from (7), we are not aware of any other lower bound with
convergence rate less than 1 for a scalar SDE.
Thispaperisorganizedasfollows. InSection2wederiveanexplicitrepresentation
of the solution of (3) and the more general case of (2). Using this representation we
show sharp upper and lower bounds of eeq in Section 3. In Section 4 we consider a
p
particularadaptivemethodthatachievesanarbitrarilyhighpolynomialconvergence
rate. Finally, we illustrate our results by numerical experiments.
2. Squared Bessel Process of Dimension One
In this section, we will derive an explicit expression for the strong solution of (3)
by using basic results about reflected SDEs. Subsequently, we will extend this
technique to the more general case of SDE (2).
In the following let (Ω,F,P) be a complete probability space and let (F ) be
t t≥0
a filtration on this space satisfying the usual conditions.
Given x ≥ 0 and a Brownian motion B w.r.t. (F ) , we define
0 t t≥0
(cid:90) t √
(11) W = sgn(B + x )dB
t s 0 s
0
for all t ≥ 0 with sgn = 1 − 1 . Then, W is a Brownian motion w.r.t.
{x>0} {x≤0}
(F ) . Indeed, the quadratic variation of W satisfies
t t≥0
(cid:90) t √ (cid:90) t
[W] = sgn(B + x )2ds = 1ds = t,
t s 0
0 0
6 HEFTER AND HERZWURM
and thus L´evy’s characterization can be applied. Now, consider the SDE (3) where
the driving Brownian motion W has the particular form given in (11). Due to this
construction of W, we see that the solution of (3) is given by
√
(12) X = (B + x )2,
t t 0
√ (cid:12) √ (cid:12)
since X = (cid:12)B + x (cid:12) and hence
t t 0
(cid:90) t(cid:112) (cid:90) t √ √
2 X dW = 2 |B + x |·sgn(B + x )dB
s s t 0 s 0 s
0 0
(cid:90) t √
= 2 (B + x )dB
s 0 s
0
√
= B2 −t+2 x B .
t 0 t
Moreover, Tanaka’s formula [21, Prop. III.6.8] given by
(cid:90) t
|B −a| = |a|+ sgn(B −a)dB +2LB(a), a ∈ R,
t s s t
0
√
where LB(a) denotes the local time of B in a, yields for a = − x that
0
√ √ √
|B + x | = x +W +2LB(− x )
t 0 0 t t 0
(cid:18) (cid:19)
√ √
= x +W +max 0, sup −( x +W ) ,
0 t 0 s
0≤s≤t
where the second equality follows from Skorokhod’s lemma, see [21, Lem. III.6.14].
Finally, using sup(−A) = −inf(A) for A ⊆ R leads to the solution of (3) given by
(cid:32) √ (cid:18) √ (cid:19)−(cid:33)2
(13) X = (W + x )+ inf W + x ,
t t 0 s 0
0≤s≤t
where x− = −min(0,x) denotes the negative part of x. We stress that the explicit
solution (13) of the SDE (3) holds for any Brownian motion W and does not depend
on the particular construction given in (11), see [21, Cor. V.3.23], since pathwise
uniqueness and strong existence holds for the SDE (3). Hence the unique strong
solution of the SDE (3) is given by (13).
Remark 1. It is well-known that the solution of the SDE (3) can be expressed in
terms of B in (12), see, e.g., [26, Ex. IX.3.16]. However, we are not aware of a result
regarding the explicit form of the strong solution given by (13).
In the context of SDEs, the somehow explicit solution of X by means of B in (12)
is rather useless for strong approximation. The concept of strong solutions entails
a functional dependence of the solution process and the input Brownian motion
appearing in the SDE, which is W in our case. We thus seek to “construct” the
solution X out of W.
Remark 2. Equation (12) clearly yields
√
(X ) =d (cid:0)(W + x )2(cid:1) ,
t t≥0 t 0 t≥0
cf. [21, Thm. III.6.17]. However, this equation is only valid in the distributional
sense and does not hold pathwise. Note that the Brownian path attains its running
minimum whenever the solution hits 0. More precisely, we have
√
(14) X = 0 ⇔ W ≤ − x ∧ W = inf W ,
t t 0 t s
0≤s≤t
OPTIMAL STRONG APPROXIMATION OF THE 1-D SQUARED BESSEL PROCESS 7
cf. Figure 2.
2
Brownian path W
solution X
1.5
1
0.5
0
−0.5
−1
−1.5
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
t
Figure 2. Brownian path and corresponding solution (13) with ini-
tial condition x = 0.5. The solution hits the boundary 0 according
0
to (14).
Remark 3. The expression
√ (cid:18) √ (cid:19)−
(W + x )+ inf W + x
t 0 s 0
0≤s≤t
appearing in (13) is known as a reflected Brownian motion (with reflecting barrier
√
at 0) starting in x ≥ 0, see [29].
0
We now turn to the more general case of SDE (2) with arbitrary b ∈ R.
Proposition 1. The unique strong solution of the SDE (2) is given by
(cid:32) (cid:18) (cid:19)−(cid:33)2
X = u +e−bt inf ebsu ,
t t 2 2 s
0≤s≤t
where u denotes the unique strong solution of the SDE
b √
(15) du = − u dt+dW , u = x , t ≥ 0.
t t t 0 0
2
In particular, the unique strong solution of the SDE (3) is given by (13).
8 HEFTER AND HERZWURM
Remark 4. NotethatthesolutiontothelinearSDE(15)iscalledOrnstein-Uhlenbeck
process and can be solved explicitly by
(cid:18)√ (cid:90) t (cid:19)
(16) u = e−bt x + ebsdW , t ≥ 0,
t 2 0 2 s
0
see, e.g., [21, Ex. V.6.8)].
Proof of Proposition 1. Analogoustotheabovederivation,westartwithaBrownian
motion B and consider the SDE
b √
duB = − uBdt+dB , uB = x , t ≥ 0.
t 2 t t 0 0
Moreover, we assume that the Brownian motion appearing in (2) has the particular
form
(cid:90) t
W = sgn(uB)dB , t ≥ 0.
t s s
0
Then, Itˆo’s formula shows
(cid:16) (cid:17) (cid:16) (cid:17)
d (cid:0)uB(cid:1)2 = 1−b(cid:0)uB(cid:1)2 dt+2uBdB
t t t t
(cid:16) (cid:17)
= 1−b(cid:0)uB(cid:1)2 dt+2(cid:12)(cid:12)uB(cid:12)(cid:12)dW .
t t t
Hence the solution of SDE (2) is given by
(17) X = (cid:0)uB(cid:1)2, t ≥ 0.
t t
On the other hand, the Tanaka-Meyer formula [21, Thm. III.7.1(v)] applied to
u˜B = uBebt
2
t t
combined with the explicit expression (16) yields
(cid:12)(cid:12)u˜Bt (cid:12)(cid:12) = √x0 +(cid:90) te2bssgn(cid:0)u˜Bs (cid:1)dBs +2Λut˜B(0)
0
√ (cid:90) t
= x0 + e2bsdWs +2Λut˜B(0),
0
where Λu˜B(0) denotes the semimartingale local time of u˜B at 0. Thus, Skorokhod’s
lemma [21, Lem. III.6.14] shows
(cid:12)(cid:12)u˜Bt (cid:12)(cid:12) = √x0 +(cid:90) te2bsdWs +(cid:18) inf √x0 +(cid:90) se2budWu(cid:19)−,
0≤s≤t
0 0
and consequently
(cid:18) (cid:19)−
(cid:12)(cid:12)uBt (cid:12)(cid:12) = ut +e−2bt inf e2bsus .
0≤s≤t
It remains to apply [21, Cor. V.3.23]. (cid:3)
Remark 5. Consider the situation of the proof of Proposition 1. The Tanaka-Meyer
formula applied to uB yields
t
(cid:12)(cid:12)uB(cid:12)(cid:12) = √x − b (cid:90) tuBsgn(cid:0)uB(cid:1)ds+(cid:90) tsgn(cid:0)uB(cid:1)dB +2ΛuB(0)
t 0 2 s s s s t
0 0
√ b (cid:90) t(cid:12) (cid:12) (cid:90) t
= x − (cid:12)uB(cid:12)ds+ dW +2ΛuB(0).
0 2 s s t
0 0
OPTIMAL STRONG APPROXIMATION OF THE 1-D SQUARED BESSEL PROCESS 9
(cid:12) (cid:12)
Hence Z = (cid:12)uB(cid:12) is the solution of the reflected SDE on the domain D = ]0,∞[
t t
given by
b √
(18) dZ = − Z dt+dW +dK , Z = x , t ≥ 0,
t t t t 0 0
2
where K is a process of bounded variation with variation increasing only when
Z ∈ ∂D = {0}. For details on reflected SDEs we refer to [29]. In view of (17), we
t
can express the solution X to the SDE (2) by
X = (Z )2, t ≥ 0.
t t
3. Equidistant Methods for SDE (3)
As shown in the previous section, the SDE (3) admits the explicit solution (13).
ThisimmediatelylinkstheproblemofapproximatingSDE(3)toglobaloptimization
under the Wiener measure. We refer to [27] for results on global optimization under
the Wiener measure.
In this section we show sharp (up to constants) upper and lower bounds for eeq
p
for all p ≥ 1. Recall that eeq denotes the n-th minimal error corresponding to
p
the approximation of SDE (3) and algorithms that are based on equidistant grids,
see (5). We show that the convergence rate of eeq is equal to 1/2. Moreover, we
p
show that the drift-implicit Euler scheme converges with the optimal rate 1/2.
For the proofs we exploit results from [3] on the asymptotic error distribution of
the infimum of a Brownian motion approximated by equidistant points on the unit
interval. For n ∈ N we define
√ (cid:18) (cid:19)
(19) δ = n· min W − inf W ,
n k s
0≤k≤n n 0≤s≤1
where W denotes a standard Brownian motion. The following result is due to [3,
Thm. 1 and Lem. 6].
Theorem 1 ([3]).
(i) The sequence (δ ) converges in distribution.
n n∈N
(ii) For all p ∈ [1,∞[ the sequence (δp) is uniformly integrable. In particular,
n n∈N
for all p ∈ [1,∞[ we have
(cid:18) (cid:18)(cid:12) (cid:12)p(cid:19)(cid:19)1/p
(20) E (cid:12)(cid:12) min Wk − inf Ws(cid:12)(cid:12) (cid:52) n−1/2.
(cid:12)0≤k≤n n 0≤s≤1 (cid:12)
Remark 6. In [3], the limiting distribution of the sequence (δ ) is given explicitly
n n∈N
by means of three-dimensional Bessel processes.
Recall that the solution of (3) is given by X = (Y )2 with
1 1
√ (cid:18) √ (cid:19)−
(21) Y = (W + x )+ inf W + x .
1 1 0 s 0
0≤s≤1
Moreover, for n ∈ N we define the approximation X(cid:98)(n) of X by
1 1
(cid:16) (cid:17)2
X(cid:98)(n) = Y(cid:98)(n) ,
1 1
where
√ (cid:18) √ (cid:19)−
Y(cid:98)(n) = (W + x )+ min W + x
1 1 0 k 0
0≤k≤n n
10 HEFTER AND HERZWURM
serves as an approximation of Y . Here, the global infimum is simply replaced by
1
the discrete minimum over n equidistant knots.
The following upper bound is a consequence of Theorem 1(ii).
Theorem 2. For every 1 ≤ p < ∞ we have
(cid:16) (cid:17)
e X(cid:98)(n) (cid:52) n−1/2.
p 1
In particular, the n-th minimal error satisfies
eeq(n) (cid:52) n−1/2.
p
Proof. At first, note that 0 ≤ Y(cid:98)(n) ≤ Y and
1 1
(cid:12) (cid:12) (cid:18) √ (cid:19)− (cid:18) √ (cid:19)−
(cid:12)Y −Y(cid:98)(n)(cid:12) = inf W + x − min W + x ≤ min W − inf W .
(cid:12) 1 1 (cid:12) s 0 k 0 k s
0≤s≤1 0≤k≤n n 0≤k≤n n 0≤s≤1
Hence we get
(cid:18) (cid:19)
(cid:12) (cid:12) (cid:16) (cid:17) (cid:12) (cid:12)
(cid:12)X −X(cid:98)(n)(cid:12) = Y +Y(cid:98)(n) ·(cid:12)Y −Y(cid:98)(n)(cid:12) ≤ 2Y · min W − inf W .
(cid:12) 1 1 (cid:12) 1 1 (cid:12) 1 1 (cid:12) 1 k s
0≤k≤n n 0≤s≤1
Finally, Cauchy-Schwarz inequality yields
E(cid:16)(cid:12)(cid:12)(cid:12)X1 −X(cid:98)1(n)(cid:12)(cid:12)(cid:12)p(cid:17) ≤ 2p ·(cid:32)E(cid:0)Y12p(cid:1)·E(cid:32)(cid:12)(cid:12)(cid:12)(cid:12)0m≤ki≤nnWnk −0≤ins≤f1Ws(cid:12)(cid:12)(cid:12)(cid:12)2p(cid:33)(cid:33)1/2
(cid:52) n−p/2
due to (20) and
E(Yr) < ∞, 1 ≤ r < ∞,
1
since Y =d (cid:12)(cid:12)W +√x (cid:12)(cid:12), see Remark 2. (cid:3)
1 1 0
The natural extension of X(cid:98)(n) to an approximation on the whole equidistant grid
1
with mesh size 1/n is given by
(cid:16) (cid:17)2
X(cid:98)(n) = Y(cid:98)(n) , k = 0,...,n,
k k
n n
where
(cid:16) √ (cid:17) (cid:18) √ (cid:19)−
Y(cid:98)(n) = W + x + min W + x , k = 0,...,n.
k k 0 l 0
n n 0≤l≤k n
Let us stress that this scheme can be expressed by the following Euler-type scheme
(cid:16) (cid:16) (cid:17) (cid:17)
(22) Y(cid:98)(n) = max Y(cid:98)(n) + W −W ,0 , k = 0,...,n−1,
k+1 k k+1 k
n n n n
√
and Y(cid:98)(n) = x . In the context of reflected SDEs, scheme (22) is simply the
0 0
projected Euler scheme for a reflected Brownian motion. Due to
√
x+|x| x+ x2
max(x,0) = = , x ∈ R,
2 2
the drift-implicit Euler scheme (8) for b = 0 reads
(cid:18) (cid:18)(cid:114) (cid:16) (cid:17) (cid:19)(cid:19)2
(23) X(cid:98)n,imp = max X(cid:98)n,imp + W −W ,0 , k = 0,...,n−1.
k+1 k k+1 k
n n n n
