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 pth 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=(Wt)t≥0 is a one-dimensional Brownian motion, κ,λ≥0, θ>0 and x0>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 (St)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.
Let 2κλ>θ2, x0>0 and T>0. Then, for all there exists a constant Kp>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 pth 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 Cp>0 depends only p,κ,λ,θ,x0 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 pth 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 (Bt)t≥0 is a Brownian motion independent of (Wt)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.
For our error analysis, we need to control the inverse moments of the CIR process. As Xt follows a non-central χ2 distribution, we have for p>−(2κλ/θ2), where 1F1 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 xq≤ceεx for x≥0. Hence, we arrive at the following.
Let 2 κλ>θ2, T>0 and q≥0. Therefore,
Furthermore, we need a smoothness result for equation (1.2).
Let 2 κλ>θ2 and T>0. Then, for all q≥1, we have and
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 qth 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
Let 2 κλ>θ2. For T>0 and 1≤p<2κλ/θ2, there exists c>0 such that for Δ∈(0,1/2].
Without loss of generality, we assume that Δ<T. For the error at time point kΔ, we have the recursion with Multiplying both sides by ek+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 yk+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.
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.
Let 2 κλ>θ2, T>0 and Δ∈(0,1/2]. Then, for all q≥1, we have
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.
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