## Abstract

We analyse the strong approximation of the Cox–Ingersoll–Ross (CIR) process in the regime where the process does not hit zero by a positivity preserving drift-implicit Euler-type method. As an error criterion, we use the *p*th mean of the maximum distance between the CIR process and its approximation on a finite time interval. We show that under mild assumptions on the parameters of the CIR process, the proposed method attains, up to a logarithmic term, the convergence of order 1/2. This agrees with the standard rate of the strong convergence for global approximations of stochastic differential equations with Lipschitz coefficients, despite the fact that the CIR process has a non-Lipschitz diffusion coefficient.

## 1. Introduction and main result

In computational finance, a lot of effort has been given to the so-called Cox–Ingersoll–Ross (CIR) process recently. The CIR process under consideration has the following form:
1.1
Here, *W*=(*W*_{t})_{t≥0} is a one-dimensional Brownian motion, *κ*,*λ*≥0, *θ*>0 and *x*_{0}>0. It is well known that equation (1.1) admits a unique strong solution that is non-negative (Karatzas & Shreve 1991, ch. 5). The CIR process was originally proposed by Cox *et al.* (1985) as a model for short-term interest rates. Nowadays, this model is widely used in financial modelling, e.g. as the volatility process in the Heston model (Heston 1993).

One of the main objectives in mathematical finance is the pricing of (path-dependent) derivatives. If the asset prices or the interest rates dynamics are modelled by a *d*-dimensional stochastic differential equation (SDE) with solution (*S*_{t})_{t∈[0,T]}, then this corresponds to the quadrature problem
where is the discounted payoff of the derivative. Typically, explicit formulae for such quantities are unknown and have to be approximated by Monte Carlo methods that are based on approximate solutions of the SDE on [0,*T*]. Here, the knowledge of the global strong convergence rate of the approximation is important, in particular, if the efficient multi-level Monte Carlo method (Giles 2007, 2008) is to be used.

The strong global approximation of equation (1.1) has been studied in several articles. Strong convergence (without a rate or with a logarithmic rate) of several discretization schemes has been shown by Deelstra & Delbaen (1998), Alfonsi (2005), Higham & Mao (2005), Bossy & Diop (2007) and Gyöngy & Rásonyi (2011). In Alfonsi (2005), a general framework for the analysis of strong approximation of the CIR process was presented, along with extensive simulation studies. Moreover, the strong approximation of more general Ait-Sahalia-type interest-rate models was analysed in Higham *et al.* (2011). However, only in Berkaoui *et al.* (2008) were non-logarithmic convergence rates obtained. There, it was shown that a symmetrized Euler method had strong convergence order 1/2 under restrictive assumptions on the parameters of the equation (see §2). The difficulties in obtaining strong convergence rates for the approximations of the CIR process are due to its square-root coefficient. Thus, the standard theory that relies on the global Lipschitz assumption does not apply (Milstein 1995; Kloeden & Platen 1999).

In this article, we focus on the regime where the CIR process does not hit zero, i.e. where
an assumption that is often fulfilled in interest-rate models. By the Feller test, it is true if and only if 2*κλ*≥*θ*^{2} (Karatzas & Shreve 1991, ch. 5). We will not directly do a numerical analysis for the CIR process, but for a coordinate transformation thereof. We consider the process , which, by Itô’s formula, satisfies
1.2
with
This transformation, known as the Lamperti transformation, allows us to shift the nonlinearity from the diffusion coefficient into the drift coefficient. Note that the drift
satisfies for *α*>0, the one-sided Lipschitz condition
This property is crucial to control the error propagation of drift-implicit Euler schemes (Higham *et al.* 2002, 2011).

The drift-implicit Euler method with stepsize *Δ*>0 for equation (1.2) leads to the numerical scheme
1.3
with and
Recalling that *α*,*γ*>0 and *β*<0, equation (1.3) has the unique positive solution
which we call the *drift-implicit square-root Euler method*. Transforming back, i.e.
gives a strictly positive approximation of the original CIR process. This scheme has already been suggested by Alfonsi (2005), but convergence results have not been established. Using piecewise linear interpolation, i.e.
we obtain a global approximation of the CIR process on [0,*T*]. Exploiting the structure of SDE (1.2), we establish the following theorem.

### Theorem 1.1

*Let 2κλ>θ*^{2}*, x*_{0}*>0 and T>0. Then, for all
*
*there exists a constant K*_{p}*>0 such that
*
*for all Δ∈(0,1/2].*

Hence, we obtain the optimal strong convergence rate for the approximation of SDEs with Lipschitz coefficients (Müller-Gronbach 2002). The only price to be paid is the restriction on *p* that arises from the need to control the inverse *p*th moments of the CIR process, which are infinite for *p*≥2*κλ*/*θ*^{2}.

The remainder of this article is structured as follows: in the next section, we give a short overview on discretization schemes for the CIR process, while the proof of theorem 1.1 is given in §3.

## 2. Numerical methods for the Cox–Ingersoll–Ross process

Discretization schemes for the CIR process have been proposed in numerous articles, including Deelstra & Delbaen (1998), Alfonsi (2005, 2010), Higham & Mao (2005), Jäckel & Kahl (2006), Bossy & Diop (2007), Moro & Schurz (2007), Berkaoui *et al.* (2008), Günther *et al.* (2008) and Gyöngy & Rásonyi (2011). In this short summary, we mainly focus on Euler-type methods.

In Berkaoui *et al.* (2008), the authors studied a symmetrized Euler method for the approximation of equation (1.1), i.e.
2.1
for *k*=0,1,…. Under the assumption
they found that
where the constant *C*_{p}>0 depends only *p*,*κ*,*λ*,*θ*,*x*_{0} and *T*. As a consequence of lemma 3.5 in §3, the piecewise linear interpolation of this scheme satisfies the same error estimate with respect to the *p*th mean maximum distance as our drift-implicit square-root Euler scheme. However, for the symmetrized Euler scheme, the assumptions on the parameters of the CIR process are clearly more restrictive. Moreover, it was shown in Bossy & Diop (2007) that the symmetrized Euler scheme converged weakly with rate 1 if 2*κλ*≥2*θ*^{2}.

The truncated Euler scheme
2.2
was analysed in Deelstra & Delbaen (1998), while the scheme
2.3
was studied in Higham & Mao (2005). Both schemes do not preserve positivity, but satisfy
for without further restrictions on the parameters of the equation. In particular, an application of Gyöngy & Rásonyi (2011, theorem 3.1) yielded a logarithmic convergence rate for the Euler schemes (2.2) and (2.3). Moreover, if 2*κλ*≥*θ*^{2}, it follows from Gyöngy (1998) that the Euler schemes (2.1)–(2.3) have a pathwise convergence rate of 1/2−*ε* for all *ε*>0, i.e. we have
for . The asymptotic error distribution of these schemes can be deduced from Neuenkirch & Zähle (2009): for 2*κλ*≥*θ*^{2}, it holds that
for , where (*B*_{t})_{t≥0} is a Brownian motion independent of (*W*_{t})_{t≥0} and (*Φ*_{t})_{t≥0} is given by

While no explicit solution of equation (1.1) is known, the finite-dimensional distributions can be characterized in terms of a non-central *χ*^{2} distribution (Cox *et al.* 1985). Thus, equation (1.1) can be simulated exactly at a finite number of time points using the Markov property (Glassermann 2004; Broadie & Kaya 2006).

However, the algorithms for the exact simulation of the CIR process are strongly problem dependent: the number of degrees of freedom of the non-central *χ*^{2} random variable, which has to be simulated in each step, is 4*κλ*/*θ*^{2}. Thus, the computational cost of the algorithms depends strongly on *κ*,*λ* and *θ*. The same problem, i.e. strong dependence of the computational cost of the algorithm on the parameters of the equation, arises also for the exact sampling algorithm introduced in Beskos & Roberts (2005), which can be also applied to equation (1.1) (Beskos *et al.* 2006).

While for the simulation of the CIR process at a single point, the exact simulation methods are useful, discretization schemes remain superior if a full sample path of the CIR process has to be simulated or if the CIR process is part of a system of SDEs (Higham & Mao 2005). Moreover, results on strong convergence rates for approximations of equation (1.1) have an interest of their own because they are one of the most prominent examples for an SDE, whose coefficients do not satisfy the standard global Lipschitz assumption.

## 3. Proof of theorem 1.1

In the following, we will denote constants by *c*, regardless of their value.

### (a) Preliminaries

For our error analysis, we need to control the inverse moments of the CIR process. As *X*_{t} follows a non-central *χ*^{2} distribution, we have
for *p*>−(2*κλ*/*θ*^{2}), where _{1}*F*_{1} denotes the confluent hypergeometric function, and
else (Hurd & Kuznetsov 2008, theorem 3.1). As
see formula 13.1.5 in Abramowitz & Stegun (1964, p. 504), it follows for *p*>−(2*κλ*/*θ*^{2}) that is bounded on [0,*T*], i.e.
3.1

Moreover, from Hurd & Kuznetsov (2008, theorem 3.1) or Bossy & Diop (2007, lemma A.2), we also have
if and only if . In general, one can estimate polynomial moments against exponential moments because for *q*,*ε*>0, there exists *c*>0 such that *x*^{q}≤*ce*^{εx} for *x*≥0. Hence, we arrive at the following.

### Lemma 3.1

*Let 2 κλ*>*θ*^{2}, *T*>0 *and q*≥0. *Therefore*,

Furthermore, we need a smoothness result for equation (1.2).

### Lemma 3.2

*Let 2 κλ*>*θ*^{2} *and T*>0. *Then, for all q*≥1, *we have*
*and*

### Proof.

For 0≤*s*≤*t*, we have
By the Cauchy–Schwarz inequality, one has
As , the first assertion follows now from equation (3.1), lemma 3.1, Minkowski’s inequality and the smoothness of Brownian motion in the *q*th mean. For the second assertion, in addition, we use the well-known fact that the modulus of continuity of Brownian motion satisfies
for *Δ*∈(0,1/2] (e.g. Fischer & Nappo 2010). The third assertion follows from
and
■

### (b) Error bound for the implicit Euler scheme for *Y*

### Proposition 3.3

*Let 2 κλ*>*θ*^{2}. *For T*>0 *and* 1≤*p*<2*κλ*/*θ*^{2}, *there exists c*>0 *such that*
*for Δ*∈(0,1/2].

### Proof.

Without loss of generality, we assume that *Δ*<*T*. For the error
at time point *kΔ*, we have the recursion
with
Multiplying both sides by *e*_{k+1}, we obtain
As
we have
and it follows that
3.2

Now, it remains to analyse the local error
Note that
for *a*,*b*>0. Thus,
and we obtain
An application of Hölder’s and Minkowski’s inequality yields
for *q*,*q*′>1 with 1/*q*+1/*q*′=1. Now, lemma 3.2 gives
3.3
for all *q*>1. As
applying (3.1) with *q*>1 such that *pq*<2*κλ*/*θ*^{2} yields
and using (3.2) completes the proof of the proposition. ■

### (c) Moment bounds for the implicit Euler scheme for *Y*

Now we show that all moments of the approximation scheme are uniformly bounded. Multiplying (1.3) by *y*_{k+1} yields
It follows that
3.4
with
3.5
and we obtain by induction that
3.6
for all *k*=0,1,…,⌈*T*/*Δ*⌉. This allows us to show the following.

### Lemma 3.4

*Let 2 κλ*≥*θ*^{2} *and Δ*>0, *T*>0. *Then, for all p*≥1, *we have*

### (d) Error bound for the drift-implicit square-root Euler scheme

Now denote by the piecewise linear interpolation of the CIR process with stepsize *Δ*>0, i.e.

### Lemma 3.5

*Let 2 κλ*>*θ*^{2}, *T*>0 *and Δ*∈(0,1/2]. *Then, for all q*≥1, *we have*

### Proof.

Combining the equality with the Cauchy–Schwarz inequality and lemma 3.2, we get Now, the assertion follows from ■

Owing to
and lemma 3.5, it only remains to control
As
an application of Hölder’s inequality with *ε*>0 such that (1+*ε*)*p*<2*κλ*/*θ*^{2} and proposition 3.3 give
It remains to apply lemmas 3.2 and 3.4 to finish the proof of theorem 1.1.

## Acknowledgements

The authors thank Martin Altmayer for valuable comments on an earlier version of the manuscript.

- Received August 22, 2011.
- Accepted November 14, 2011.

- This journal is © 2011 The Royal Society