## 1. Introduction

The Cox-Ingersoll-Ross (CIR) process is a short rate model used in the pricing of interest rate derivatives and is given by the following Itô-type stochastic differential equation (SDE),

(1) |

where is a Wiener Process and and are positive parameters, and for some fixed . Solutions of (1) are almost surely (a.s.) non-negative: in general they can achieve a value of zero but will be reflected back into the positive half of the real line immediately. Moreover, if , referred to as the Feller condition, solutions will be a.s. positive. No closed form solution of (1) is available, though it is known that has, conditional upon for

, a non-central chi-square distribution with

and Var.For Monte Carlo estimation, exact sampling from the known conditional distribution of

is possible but computationally inefficient and potentially restrictive if the innovating Wiener process of (1) is correlated with that of another process: see [1, 5, 9, 19]. Consequently a substantial literature has developed on the efficient numerical approximation of solutions of (1); we now highlight the parts which are relevant to our analysis.An approach that seeks to directly discretise (1) using some variant of the explicit Euler-Maruyama method leads to schemes of the form

(2) |

for given functions , and . These functions are selected to ensure that the diffusion coefficient remains real-valued (so that (2) is well defined), and in some cases to preserve the positivity of solutions. This approach seeks to accommodate the non-Lipschitz (square-root) diffusion, which facilitates overshoot when the solutions are close to zero, but it introduces additional bias to the approximation. A survey of choices common in practice may be found in [19], and we summarise them in Table 1 using the convention .

Method | |||
---|---|---|---|

Explicit Euler | |||

Absorption | |||

Deelstra & Delbaen [7] | |||

Lord et al [19] | |||

Reflection | |||

Higham & Mao [10] | |||

Diop [4] |

An alternative approach is to transform (1) before discretisation to make the diffusion coefficient globally Lipschitz. For example, applying the Lamperti transform yields an auxiliary SDE in with a state independent and therefore globally Lipschitz diffusion, but a drift coefficient that is unbounded when solutions are in a neighbourhood of zero. This approach is effective: a fully implicit Euler discretisation over a uniform mesh that preserves positivity of solutions was proposed in [1]

and shown to have uniformly bounded moments. A continuous time extension interpolating linearly between mesh points was shown to have strong

order of convergence (up to a factor of ) in [8] when , and a continuous-time variant based on the same implicit discretisation was shown to have strong convergence of order one when in [2].In this article we will show that a strongly convergent numerical scheme can be constructed by an application of the Lamperti transform to (1) followed by an explicit Euler-Maruyama discretisation over a procedurally generated adaptive mesh. The purpose of the adaptivity is to manage the nonlinear drift response of the discrete tranformed system, a framework for which was introduced in [13] for SDEs with one-sided Lipschitz drift and globally Lipschitz diffusion and extended to allow for monotone coefficients and a Milstein-type discretisation in [14, 15] respectively. This framework imposes maximum and minimum timesteps and in a fixed ratio and requires the use of a backstop numerical method in the event that the timestepping strategy attempts to select a stepsize below . As in [15], we will use here path-bounded strategies, this time designed to increase the density of timesteps when solutions approach zero, and we additionally require the backstop method to retake a step when the adaptive strategy overshoots zero. This latter is carried out without discarding samples from the Brownian path (preserving the trajectory), and without bridging (preserving efficiency).

We prove, under parameter constraints that imply the Feller condition (specifically, ), that the order of strong convergence in is at least . This parameter constraint is technical, and ensures sufficiently many finite conditional inverse moments of solutions of (1) (as described by [4]). We separately prove that, under exactly the Feller condition, the probability of invoking the backstop method to avoid negative values can be made arbitrarily small by choosing sufficiently small, and provide a practical method for doing so given a user defined tolerance level. The proof relies upon a finite partitioning of the sample space of trajectories induced by and

, which allows us to handle the randomness of the number of timesteps via the Law of Total Probability.

Numerically we compare the convergence and efficiency of our hybrid adaptive method with a semi-implicit adaptive variant, the fixed-step explicit method due to [10], and the transformed implicit fixed-step method proposed and analysed in [1, 8, 2], examining the parameter dependence of the numerical order of convergence in each case. The numerical convergence rates of adaptive methods are seen to outperform those of fixed-step methods over the entire domain where Feller’s condition holds.

The structure of the article is as follows. In Section 2 we give the form of the SDE governing the Lamperti transform of (1), specify the constraints placed upon the parameters for the main strong convergence theorem, and examine the availability of conditional moment and inverse moment bounds under these constraints. In Section 3 we set up the framework for our random mesh, characterise the class of path-bounded timestepping strategies and define our adaptive numerical method. In Section 4 we present the two main theorems on strong convergence and positivity, providing illustrative examples in the latter case. Finally in Section 5 we numerically compare convergence and efficiency of several commonly used methods.

## 2. Mathematical Preliminaries

Throughout this article we let be the natural filtration of . By using Itô’s Formula and applying the transformation we get,

By then setting,

we can write,

(3) |

where is not globally Lipschitz continuous, but when it satisfies a one-sided Lipschitz condition with constant :

which can be seen by noting that

Meanwhile the diffusion coefficient is constant and is therefore globally Lipschitz continuous. The SDE (3) has integral form

(4) |

In order to ensure the a.s. positivity of solutions of (3) and the boundedness of certain inverse moments of solutions of (3), we will also need the following assumption:

###### Assumption 1.

Suppose that

(5) |

Eq. (5) implies the Feller Condition (; see, for example [20]), which ensures that solutions of (1), and therefore (3), remain positive with probability one:

Assumption 1 provides inverse moment bounds as follows:

###### Lemma 2.

###### Proof.

Let be a solution of (1) where Assumption 1 holds. By Lemma A.1 in Bossy & Diop [4],

(7) |

for some and any such that . Assumption 1 ensures that , and since , (6) follows by Lemma A.1 in [4] as it applies to conditional expectations, the former requiring an additional application of Jensen’s inequality to the first inequality in (7). ∎

We also need the following bounds on positive moments of solutions of (3), which apply under Feller’s condition and in particular under Assumption 1.

###### Lemma 3.

###### Proof.

Finally, we will make frequent use of the following elementary inequalities: for and and ,

(10) | |||||

(11) | |||||

(12) |

## 3. An Adaptive Numerical Method

[13] provided a framework within which to construct timestepping strategies for an adaptive explicit Euler-Maruyama numerical scheme applied to nonlinear Itô-type SDEs of the form

(13) |

over a random mesh on the interval given by,

(14) |

where is a sequence of random timesteps and with . The random time step is determined by .

For completeness, we present here the essential elements of that framework, before discussing the dynamical considerations specific to the transformed CIR model (3) corresponding to (13) with

(15) |

which lead to our proposed timestepping strategy.

### 3.1. Framework for a random mesh

###### Definition 4 ([18]).

Suppose that each member of the sequence is an -stopping time: i.e. , where is the natural filtration of . We may then define a discrete time filtration by

###### Assumption 5.

Each is -measurable, and is a random integer such that,

and the length of maximum and minimum stepsizes satisfy , and

(16) |

###### Remark 6.

The lower bound ensures that a simulation over the interval can be completed in a finite number of timesteps, and the upper bound

prevents stepsizes from becoming too large. The latter is used as a convergence parameter in our examination of the strong convergence of the adaptive method. The random variable

cannot take values outside the finite set , where and .is a Wiener increment over a random interval the length of which depends on , through which it depends on . Therefore is not independent of

; indeed it is not necessarily normally distributed. Since

is a bounded -stopping time and -measurable, then is -conditionally normally distributed, by Doob’s optional sampling theorem (see for example [22]), and for all there exists such that(17) |

### 3.2. Adaptive timestepping strategy

To ensure strong convergence, our strategy will be to reduce the size of each timestep if discretised solutions attempt to enter a neighbourhood of zero. If we wish to control the likelihood of invoking the backstop to avoid negative values, we will also reduce the timestep when solutions grow large.

###### Definition 7 (A path-bounded time-stepping strategy).

We now give two examples of path-bounded strategies that are valid for (14), the first with infinite (which is sufficient to ensure strong convergence), and the second with finite (which is useful if we also wish to minimise the use of the backstop to ensure positivity).

###### Lemma 8.

###### Proof.

Note that for the strategies defined by (19) and (20), solutions of (14) cannot enter the neigbourhood , and therefore terms of the sequence are uniformly bounded from above. This has the effect of controlling inverse moments of the solutions of (14). Moreover for (19), when (15) holds,

and therefore that strategy is admissible in the sense of [13, Definition 2.2] with and . Similarly, (20) is admissible with and .

### 3.3. The adaptive numerical method with backstop

We consider an adaptive scheme based upon the following explicit Euler-Maruyama discretisation of (3) over a random mesh given by,

(21) |

where the timestep sequence is constructed according to (19). For , the continuous version is given by

(22) |

so that for each .

We combine this scheme with a positivity-preserving backstop scheme that is to be applied if the timestepping strategy attempts to select a timestep below (in which case we choose ) or if the current selected timestep and subsequently observed Brownian increment would result in the approximation becoming negative:

###### Definition 9.

Then we define the continuous form of an adaptive explicit Euler scheme to be the sequence of functions obeying

(23) |

for , , where , satisfies the conditions of Assumption 5, and the map satisfies

(24) |

for some non-negative constants and , independent of N, and

(25) |

###### Remark 10.

Since the events and are -measurable but not -measurable, if a negative value of is observed following a step of length we must retake the step using the backstop method, which will ensure positivity over that step by (25). This introduces an element of backtracking into the algorithm, but as long as the originally computed stepsize and Brownian increment are retained we can stay on the same trajectory while avoiding the use of a Brownian bridge. Theorem 15 in Section 4.2, illustrated by Example 16, demonstrates that it is always possible to choose to ensure that this particular use of the backstop can be avoided with probability , for arbitrarily small , on each trajectory.

## 4. Main Results

In this section, we first demonstrate strong convergence of solutions of (23) to those of (3) under Assumption 1 and a path-bounded timestepping strategy. Second, we investigate the likelihood that the adaptive part of the method generates a negative value (triggering the use of the backstop to ensure positivity) and show how may be chosen to control the probability of this occurring.

### 4.1. Strong convergence of the adaptive method with path-bounded timestepping strategy

###### Lemma 11.

Let be a solution of (3) and let be a random mesh such that each is an -stopping time. Fix and suppose that , where . Then, for any , we have

where

is an -measurable random variable with finite expectation, and , are the constants defined by (6) and (8) in the statements of Lemmas 2 and 3 respectively .

###### Proof.

Solutions of (3) satisfy the integral equation

and therefore

Using the triangle and Cauchy-Schwarz inequalities, and the elementary inequality (12) with ,

for . Now apply conditional expectations on both sides with respect to and (17) to get,

where we have used (6) and (8) from the statement of Lemma 2 at the last step. Therefore

as required.

∎

###### Lemma 12.

Let be a solution of (3) and let Assumption 1 hold. Let arise from the adaptive timestepping strategy satisfying (18) in Definition 7 for some , and formulate the Taylor expansion of around , where is as given in (15), as

(26) |

where

(27) |

For any , the conditional moment of

satisfies

(28) |

where is an a.s. finite and -measurable random variable given by

Moreover, there exists independent of such that

(29) |

###### Proof.

By direct substitution of from (15) into (27), and taking the -moment conditional upon , we get

Using the triangle inequality and (12) we get

Next apply Lemma 11 followed by the Cauchy-Schwarz inequality to get

###### Lemma 13.

###### Proof.

For , we subtract (22) from (4) to get

(31) | |||||

where is defined as in (15) and . Applying the Itô formula and setting , we can write,

By (26) in the statement of Lemma 12, , where is defined in (27). This, and an application of (11) gives

(32) | |||||

Consider the third term on the RHS of (32), and substitute (31) into the integrand:

(33) | |||||