Strong order 1/2 convergence of full truncation Eulerapproximations to the Cox–Ingersoll–Ross process

# Strong order 1/2 convergence of full truncation Euler approximations to the Cox–Ingersoll–Ross process

## Abstract

We study convergence properties of the full truncation Euler scheme for the Cox–Ingersoll–Ross process in the regime where the boundary point zero is inaccessible. Under some conditions on the model parameters (precisely, when the Feller ratio is greater than three), we establish the strong order 1/2 convergence in of the scheme to the exact solution. This is consistent with the optimal rate of strong convergence for Euler approximations of stochastic differential equations with globally Lipschitz coefficients, despite the fact that the diffusion coefficient in the Cox–Ingersoll–Ross model is not Lipschitz.

Keywords: Cox–Ingersoll–Ross process, strong convergence order, explicit full truncation Euler scheme

Mathematics Subject Classification (2010): 60H35, 65C05, 65C30

\WarningFilter

latexText page

## 1 Introduction

Let be a filtered probability space and let be a 1-dimensional -adapted Brownian motion. A Cox–Ingersoll–Ross (CIR) process is defined by the stochastic differential equation (SDE):

 dvt=k(θ−vt)dt+ξ√vtdWvt, (1.1)

where , , and are strictly positive real numbers. The SDE admits a unique strong solution, which is strictly positive when by the Feller test [24]. The CIR process was originally introduced in finance to model the short-term interest rate [8]. Nowadays, due to its desirable properties, like non-negativity, mean-reversion and analytical tractability, it plays a key role in the field of option pricing, for instance when modeling squared volatilities in the Heston model [18]. For a given , the CIR process has a noncentral chi-squared conditional distribution and its increments can be simulated exactly [6]. However, when pricing financial derivatives written on an underlying process modeled by a -dimensional SDE, with CIR dynamics in one or more components, we need to evaluate

 E[f(S)], (1.2)

where is the discounted payoff. This expectation is rarely available in closed form, nor can be sampled exactly. In this case, we employ Monte Carlo simulation methods [14] and approximate the solution to the SDE using a suitable discretization scheme. Due to the non-zero probability of the approximation process becoming negative, the standard Euler–Maruyama scheme applied to (1.1) is not well-defined. Possible remedies are to set the process equal to zero when it turns negative (absorption fix), e.g., the full truncation Euler (FTE) scheme studied in this paper, or reflect it in the origin (reflection fix). An overview of the explicit Euler schemes with different fixes at the boundary can be found in [28]. Alternatively, we can use an implicit Euler or a Milstein scheme to discretize the CIR process.

Weak convergence is important when estimating expectations of payoffs such as the one in (1.2). However, strong convergence plays a crucial role in multilevel Monte Carlo methods [12, 13, 25] and may be required for some complex path-dependent derivatives. Furthermore, pathwise convergence follows automatically [26].

The classical convergence theory [20, 27] does not apply to the CIR process because the square-root diffusion coefficient is not Lipschitz. Consequently, a considerable amount of research has been devoted to the numerical approximation of (1.1) and alternative approaches have been employed to prove the strong convergence of various discretizations for the CIR process.

Results in the literature concerned with positive strong convergence rates for numerical approximations of the CIR process were restricted to the regime where the boundary point is inaccessible until just recently. Strong convergence, at best with a logarithmic rate, of different Euler schemes including the partial truncation, the full truncation, the reflection and the symmetrized Euler schemes was established in [1, 10, 15, 19, 28]. The first non-logarithmic rate was obtained in [3], where it was shown that the symmetrized Euler scheme converges strongly with the standard order 1/2 to the exact solution, although in a very restricted parameter regime. Strong convergence with order 1/2 of the backward (drift-implicit) Euler–Maruyama (BEM) scheme was later established in [11], and this rate was improved to 1 in [2, 30]. Recently, [5] proved the strong convergence with order 1 of the symmetrized Milstein scheme under some restrictive conditions on the parameters, whereas [7] proved the strong convergence of a modified Euler–Maruyama scheme with an order between 1/6 and 1 that depends on the parameters.

In the last few years, there has been some development in the accessible boundary case and polynomial rates of strong convergence for an order of up to 1/2 were established in [16] and [23] for the truncated Milstein scheme and the BEM scheme, respectively.

In this paper, we study the full truncation Euler scheme proposed in [28]. This scheme preserves the positivity of the original process, is easy to implement and hence arguably the most widely used scheme in practice. Perhaps most importantly, it is found empirically to produce the smallest bias of all explicit Euler schemes with different fixes at the boundary [28]. Consider a uniform grid

 T=Nδt,tn=nδt,∀n∈{0,1,...,N}. (1.3)

We introduce the discrete-time auxiliary process

 ~vtn+1=~vtn+k(θ−~v+tn)δt+ξ√~v+tnδWvtn,~v0=v0, (1.4)

where and , its continuous-time interpolation

 ~vt=~vtn+k(θ−~v+tn)(t−tn)+ξ√~v+tn(Wvt−Wvtn) (1.5)

as well as the non-negative process

 ¯vt=~v+tn, (1.6)

whenever . The convergence in of this scheme was proved in [28]. The convergence rate, however, remained an open question and our Theorem 1.1 is the first result to address it, to the best of our knowledge. For convenience, define the Feller ratio

 ν=2kθξ2. (1.7)

We establish the strong convergence in with order 1/2 of the scheme in the inaccessible boundary case, specifically, for a Feller ratio above three. Hence, we obtain the optimal strong convergence rate for the numerical approximation of SDEs with globally Lipschitz coefficients [21, 29]. The main and novel idea of the proof is to weight the difference between the process and its approximation by the former raised to a suitably chosen negative power and prove the strong convergence with a rate of the weighted error. This, in turn, allows us to derive an upper bound for the actual error.

###### Theorem 1.1.

Suppose that and let . Then the FTE scheme converges strongly in with order 1/2, i.e., there exist and a constant such that, for all ,

 supt∈[0,T]E[|vt−¯vt|p]1p≤CN−0.5. (1.8)

We mention that the assumption on the Feller ratio from Theorem 1.1, i.e., , appears in the literature as a sufficient condition for the strong convergence with a rate of several discretization schemes for the CIR process. For example, this condition ensures in Corollary 4.1 in [7] the strong convergence with order 1/2 (and order 1 if ) of the modified Euler–Maruyama scheme, and in Proposition 3.1 in [30] the strong convergence with order 1 of the BEM scheme.

In [17], a lower error bound was recently derived for all discretization schemes for the CIR process based on equidistant evaluations of the Brownian driver in the regime where the boundary point is accessible. In light of this result, we demonstrate numerically that the FTE scheme achieves an optimal performance – in the sense – in half of the regime where the boundary point is accessible, where by optimal we mean that the empirical convergence rate is the best possible for equidistant discretization schemes for the CIR process.

The remainder of this paper is structured as follows. In Section 2, we prove the convergence of the scheme. In Section 3, we conduct numerical tests for the rate of convergence that validate and complement our theoretical findings. Finally, Section 4 contains a short discussion.

## 2 Convergence analysis

We need to control the polynomial moments of the CIR process and its FTE discretization.

###### Lemma 2.1.

The CIR process has bounded moments, i.e.,

 Missing dimension or its units for \hskip (2.1)
###### Proof.

Follows from [11] or Theorem 3.1 in [22]. ∎

###### Lemma 2.2.

The FTE scheme has uniformly bounded moments, i.e.,

 supN≥1E[supt∈[0,T]|~vt|p]<∞,∀p≥1. (2.2)
###### Proof.

Follows from a simple application of the Burkholder–Davis–Gundy (BDG) inequality and Proposition 3.7 in [9]. ∎

By construction, the FTE approximation is non-negative. However, an important step in the convergence analysis lies in analyzing the behaviour of the auxiliary process at the boundary. The next result derives a polynomial upper bound in the time step size on the probability of becoming negative. Similar results were established for the symmetrized Euler scheme, in Lemma 3.7 in [4], and for the symmetrized Milstein scheme, in Lemma 2.2 in [5]. However, the full truncation Euler scheme has led to different technical challenges and the arguments employed in the proofs of the aforementioned results do not apply here.

###### Proposition 2.3.

Suppose that and let

 ¯ν=inf{x>0:(4x∨ν)(ν−x)νx(ν−x−1)<1.99√νπeν2−1}. (2.3)

Then there exist and a constant such that, for all ,

 sup0≤n≤NP(~vtn≤0)≤CN−ν+¯ν+1. (2.4)
###### Proof.

We first show that exists and derive bounds which imply that the exponent in (2.4) is negative. Note that as the value of increases, the left-hand side term of the inequality in (2.3) decreases and the right-hand side term increases. Hence, decreases as increases. In particular, .

Suppose that and fix . Define, for brevity,

 αN=12(1−kTN)=1−kδt2. (2.5)

First, consider the sequence given by

 Missing dimension or its units for \hskip (2.6)

As , one can clearly see that for all . We will show by induction that

 Missing dimension or its units for \hskip (2.7)

where

 φν=1−¯νν and ην=(ν−¯ν)(4¯ν∨ν)ν¯ν. (2.8)

Since

 φνην=¯ν4¯ν∨ν≤14≤(1−αN)2, (2.9)

(2.7) holds when . Suppose that (2.7) holds for some , then

 cj+1≤αN−α2N+(1−αN−φνj−1+ην)2 (2.10)

and some simple computations lead to the following sufficient condition for the induction step,

 (j−1+ην)2(1−2αN)+(j−1+ην)(1−2αN)+(j−1)(1−φν)+ην(1−φν)−φν≥0, (2.11)

which clearly holds. For convenience, define another sequence given by

 aj=2(αN−cj)ξ2δt,∀0≤j≤N, (2.12)

such that

 Missing dimension or its units for \hskip (2.13)

A sequence similar to was analyzed in [4]. However, a sharper lower bound than the one obtained in Lemma 3.6 in [4] is needed for our purposes. We now show that, for large enough,

 S1=N−2∑n=0n∏j=0exp{−ajkθδt}≤2√ν(1−kδt)πeν2(1−kδt), (2.14)

a bound which will be of use later in the proof. Using (2.7), we get

 S1 =N−2∑n=0exp{νn∑j=1(cj−αN)} ≤N−2∑n=0exp{νn∑j=1(1−2αN−φνj−1+ην)} ≤N−2∑n=0exp{νkTnN−νφν∫n0(x+ην)−1dx} ≤ηνφννN−1∑n=0eνkTnN(ην+n)−νφν. (2.15)

Let for the rest of the proof. Since , the Hurwitz (generalized Riemann) zeta function

 ζ(νφν,ην)=∞∑n=0(ην+n)−νφν (2.16)

converges and hence, for large enough, we have

 N−1∑n=1(eνkTnN−1−ϵ)(ην+n)−νφν≤∑n>log(1+ϵ)νkTN0(eνkT−1)(ην+n)−νφν≤ϵη−νφνν, (2.17)

which implies that

 N−1∑n=0eνkTnN(ην+n)−νφν≤(1+ϵ)N−1∑n=0(ην+n)−νφν≤(1+ϵ)ζ(νφν,ην). (2.18)

However,

 ζ(νφν,ην)=η−νφνν+∞∑n=1(ην+n)−νφν≤η−νφνν+∫∞0(ην+x)−νφνdx, (2.19)

and hence

 ζ(νφν,ην)≤η−νφνν+η−νφν+1ννφν−1. (2.20)

Combining (2), (2.18) and (2.20) and using (2.3) and (2.8), we deduce that

 S1≤1.99(1+ϵ)√νπeν2. (2.21)

However, for large enough, we have

 √νπeν2≤(1+ϵ)√ν(1−kδt)πeν2(1−kδt) (2.22)

and the upper bound on in (2.14) follows.

Second, recall from (1.4) that, for all ,

 ~vtN−j=~vtN−j−1+kθδt−k~v+tN−j−1δt+ξ√~v+tN−j−1(WvtN−j−WvtN−j−1), (2.23)

and note that

 P(~vtN−j≤0)=P(~vtN−j−1≤−kθδt)+P(~vtN−j≤0,~vtN−j−1>0). (2.24)

Let be the natural filtration generated by and consider the shorthand notations and for the conditional expectation and probability. Conditioning on , we get

 P(~vtN−j≤0,~vtN−j−1>0) =E[\mathds1~vtN−j−1>0EtN−j−1[\mathds1~vtN−j≤0]] =E[\mathds1w>0PtN−j−1(Z≤−kθδt+w+(1−kδt)ξ√w+δt)], (2.25)

where and independent of . Using a standard inequality for the lower tail of the normal distribution, namely

 Missing dimension or its units for \hskip (2.26)

and the arithmetic mean-geometric mean (AM-GM) inequality, we deduce that

 P(~vtN−j≤0,~vtN−j−1>0) ≤12√ν(1−kδt)πE[\mathds1w>0exp{−(kθδt+w+(1−kδt))22ξ2w+δt}] ≤e−ν2(1−kδt)2√ν(1−kδt)πE[exp{−a1abs∞(~vtN−j−1)}], (2.27)

where if and otherwise. Let , then

 E[exp{−aiabs∞(~vtN−j)}] =E[exp{−aiabs∞(~vtN−j)}(\mathds1~vtN−j−1>0+\mathds1~vtN−j−1≤0)] ≤E[\mathds1~vtN−j−1>0EtN−j−1[exp{−aiabs∞(~vtN−j)}]] +P(−kθδt<~vtN−j−1≤0). (2.28)

Denote as before and let be the inner expectation on the right-hand side of (2), i.e.,

 I=EtN−j−1[exp{−aiabs∞(w+kθδt−kw+δt+ξ√w+δtZ)}], (2.29)

where independent of . There are two possible outcomes, namely , in which case , and , which is treated now:

 I Missing dimension or its units for \hskip =exp{−aikθδt−wai(1−kδt)+12wa2iξ2δt}Φ(−z0−aiξ√wδt), (2.30)

where

 z0=−kθδt+w(1−kδt)ξ√wδt (2.31)

and is the standard normal CDF. We deduce from (2.13) that

 I≤exp{−aikθδt−wai+1}, (2.32)

and hence that

 E[exp{−aiabs∞(~vtN−j)}] Extra open brace or missing close brace +P(−kθδt<~vtN−j−1≤0). (2.33)

For large enough, putting together (2.14), (2.24), (2) and (2), we can prove by induction that, for all ,

 sup0≤n≤lP(~vtN−n≤0) ≤e−ν2(1−kδt)2√ν(1−kδt)πl∑n=0E[exp{−an+1abs∞(~vtN−l−1)}]n∏j=0exp{−ajkθδt} Missing or unrecognized delimiter for \big (2.34)

Taking in (2) and since , we obtain

 sup0≤n≤NP(~vtn≤0)≤e−ν2(1−kδt)2√ν(1−kδt)πN−1∑n=0exp{−an+1v0}n∏j=0exp{−ajkθδt}. (2.35)

Using (2.7) yields

 S2 =N−1∑n=0exp{−an+1v0}n∏j=0exp{−ajkθδt} Missing or unrecognized delimiter for \Bigg ≤ηνφννexp{ν(v0θ+kT)}N−1∑n=0(ην+n)−νφνexp{−2v0φνNξ2T(n+ην)}. (2.36)

However, for all , we have that such that , and hence

 S2 ≤ηνφννexp{ν(v0θ+kT)}N−1∑n=0e−νφν(kθTv0)νφνN−νφν Missing or unrecognized delimiter for \right (2.37)

Combining (2.22), (2.35) and (2), we conclude that there exists such that

 sup0≤n≤NP(~vtn≤0)≤(1+ϵ)e−ν22√νπ(kθTηνv0)νφνexp{ν(v0θ+kT−φν)}N−ν+¯ν+1 (2.38)

for all . ∎

Next, we bound the difference between the two continuous-time approximations and .

###### Proposition 2.4.

Suppose that and let . Then there exist and a constant such that, for all ,

 supt∈[0,T]E[|~vt−¯vt|p]1p≤CN−0.5. (2.39)
###### Proof.

Let . For convenience of notation, define

 Δ~vt=~vt−~v¯t, (2.40)

where for all . From the triangle inequality, we have

 |~vt−¯vt|≤|~v¯t−~v+¯t|+|Δ~vt|=|~v¯t|\mathds1~v¯t≤0+|Δ~vt| (2.41)

and hence, using Hölder’s inequality, we get

 supt∈[0,T]E[|~vt−¯vt|p] ≤2p−1supN≥1sup0≤n≤NE[|~vtn|p(p+ϵ)ϵ]ϵp+ϵsup0≤n≤NP(~vtn≤0)pp+ϵ +2p−1supt∈[0,T]E[|Δ~vt|p]. (2.42)

We can bound the last term on the right-hand side from above as follows,

 |Δ~vt|p ≤(kθδt+k¯vtδt+ξ√¯vt∣∣Wvt−Wv¯t∣∣)p ≤3p−1kpθp(δt)p+3p−1kp|~v¯t|p(δt)p+3p−1ξp|~v¯t|p2∣∣Wvt−Wv¯t∣∣p, (2.43)

and hence

 supt∈[0,T]E[|Δ~vt|p] ≤3p−1(kθT)pN−p+3p−1(kT)psupN≥1sup0≤n≤NE[|~vtn|p]N−p Extra open brace or missing close brace (2.44)

where . Substituting back into (2) with (2) and using Lemma 2.2 and Proposition 2.3 concludes the proof. ∎

Before we prove the main result of this paper, we need the following technical auxiliary result.

###### Lemma 2.5.

Suppose that . Then we can find such that

 ν>2β+1>ν+2q−√(ν+2q−1)2−4q(q−1) (2.45)

and

 2(ν−¯ν−1)(ν−β−1)>νq. (2.46)
###### Proof.

Note that (2.45) and (2.46) are equivalent to

 (νq)∨(ν−1)(ν−¯ν−1)<(ν−¯ν−1)(ν−2q−1+√(ν+2q−1)2−4q(q−1)) (2.47)

since spans the interval of values associated with (2.47) for

 2(ν−¯ν−1)(ν−β−1). (2.48)

Furthermore, one can easily check that

 2q<√(ν+2q−1)2−4q(q−1), (2.49)

which implies that

 1+¯ν<ν−νqν−2q−1+√(ν+2q−1)2−4q(q−1) (2.50)

is equivalent to (2.47). However, we can rewrite (2.50) as

 1+¯ν<ν(1−2q+1−ν+√(ν+2q−1)2−4q(q−1)4(2ν−q−1)). (2.51)

Let and define the right-hand side function

 fq(x) =(q+1+x)(1−4q−(3q+x)+√(3q+x)2−4q(q−1)4(q+1+2x)) =1+(2q+1+2x)xq+1+2x+(2q+2+2x)q(q−1)2(q+1+2x)(3q+x+√(3q+x)2−4q(q−1)). (2.52)

Hence, we deduce that

 fq(x)≥1+x+q(q−1)4(3q+x)≥1+q−112≥1.083>1.057=1+¯ν|ν=3≥1+¯ν, (2.53)

which concludes the proof. ∎

With these results at our disposal, we are now ready to prove the main theorem.

###### Proof.

(Theorem 1.1) Let and fix any which satisfies (2.45) and (2.46). For convenience of notation, define

 evt=vt−~vt,ev0=0, (2.54)

such that

 devt=−k(vt−¯vt)dt+ξ(√vt−√¯vt)dWvt. (2.55)

Since , the CIR process has almost surely strictly positive paths. Applying Itô’s formula to the function yields

 v−βt|evt|q =−β∫t0v−(β+1)u|evu|qdvu+q∫t0v−βu|evu|q−1sgn(evu)devu+12β(β+1)∫t0v−(β+2)u|evu|qd⟨v⟩u +12q(q−1)∫t0v−βu|evu|q−2d⟨ev⟩u−βq∫t0v−(β+1)u|evu|q−1sgn(evu)d⟨v,ev⟩u, (2.56)

where if and otherwise, and hence

 v−βt|evt|q =−βkθ∫t0v−(β+1)u|evu|qdu+βk∫t0v−βu|evu|qdu−βξ∫t0v−(β+0.5)u|evu|qdWvu −qk∫t0v−βu|evu|q−1sgn(evu)(vu−¯vu)du+qξ∫t0v−βu|evu|q−1sgn(evu)(√vu−√¯vu)dWvu Missing dimension or its units for \hskip −βqξ2∫t0v−(β+0.5)u|evu|q−1sgn(evu)(√vu−√¯vu)du. (2.57)

We can show that the two stochastic integrals in (2) are true martingales by a simple application of Hölder’s inequality and Lemmas 2.1, 2.2 and 2.5. Taking expectations on both sides, since

 vu−¯vu=evu+~vu−¯vu and sgn(evu)evu=|evu|, (2.58)

we deduce that

 Missing or unrecognized delimiter for \Big =(β−q)kE[∫t0v−βu|evu|qdu]−qkE[∫t0v−βu|evu|q−1sgn(evu)(~vu−¯vu)du] −βqξ2E[∫t0v−(β+0.5)u|evu|q−2evu(√vu−√¯vu)du] +12q(q−1)ξ2E[∫t0v−βu|evu|q−2∣∣√vu−√¯vu∣∣2du] −12βξ2(ν−β−1)E[∫