# Hybrid, adaptive, and positivity preserving numerical methods for the Cox-Ingersoll-Ross model

We introduce an adaptive Euler method for the approximate solution of the Cox-Ingersoll-Ross short rate model. An explicit discretisation is applied over an adaptive mesh to the stochastic differential equation (SDE) governing the square root of the solution, relying upon a class of path-bounded timestepping strategies which work by reducing the stepsize as solutions approach a neighbourhood of zero. The method is hybrid in the sense that a backstop method is invoked if the timestep becomes too small, or to prevent solutions from overshooting zero and becoming negative. Under parameter constraints that imply Feller's condition, we prove that such a scheme is strongly convergent, of order at least 1/2. Under Feller's condition we also prove that the probability of ever needing the backstop method to prevent a negative value can be made arbitrarily small. Numerically, we compare this adaptive method to fixed step schemes extant in the literature, both implicit and explicit, and a novel semi-implicit adaptive variant. We observe that the adaptive approach leads to methods that are competitive over the entire domain of Feller's condition.

• 3 publications
• 4 publications
• 1 publication
08/31/2019

### Strong convergence of an adaptive time-stepping Milstein method for SDEs with one-sided Lipschitz drift

We introduce explicit adaptive Milstein methods for stochastic different...
12/31/2021

### A Strongly Monotonic Polygonal Euler Scheme

Rate of convergence results are presented for a new class of explicit Eu...
05/20/2021

### A flexible split-step scheme for MV-SDEs

We present an implicit Split-Step explicit Euler type Method (dubbed SSM...
09/08/2019

### Semi-explicit discretization schemes for weakly-coupled elliptic-parabolic problems

We prove first-order convergence of the semi-explicit Euler scheme combi...
09/22/2019

### Contractivity of Runge-Kutta methods for convex gradient systems

We consider the application of Runge-Kutta (RK) methods to gradient syst...
12/25/2019

### Enforcing strong stability of explicit Runge–Kutta methods with superviscosity

A time discretization method is called strongly stable, if the norm of i...
05/08/2019

### Modern theory of hydraulic fracture modeling with using explicit and implicit schemes

The paper presents novel results, obtained on the basis of the modified ...

## 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) dX(t)=κ(λ−X(t))dt+σ√X(t)dW(t), t∈[0,T]; X(0)=X0>0,

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) Vn+1=g0(Vn+Δt(κ(λ−g1(Vn))+σ√g2(Vn))ΔWn

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 .

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,

 α=4κλ−σ28, β=−4κ8, γ=σ2,

we can write,

 (3) dY(t)=(αY(t)+βY(t))dt+γdWt, t∈[0,T];Y(0)=√X0∈R+,

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

 [f(x)−f(y)](x−y)≤β(x−y)2, for all x,y∈R+,

which can be seen by noting that

 f(x)−f(y)=(x−y)[β−αxy].

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

 (4) Y(t)=Y(0)+∫t0(αY(s)+βY(s))ds+∫t0γdWt,t≥0.

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) κλ>2σ2.

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

 P[Y(t)>0, t≥0]=1.

Assumption 1 provides inverse moment bounds as follows:

###### Lemma 2.

Let be a solution of (3), where Assumption 1 holds, and let . For any , and for , there exists such that

 (6) E[1Y(s)p∣∣∣Ft]≤C(p,T)Y(t)p,a.s.
###### Proof.

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

 (7) E[1X(t)]≤eκtX0% andE[1X(t)p]≤C(p,T)Xp0,

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.

Let be a solution of (3), where Assumption 1 holds, and let . For any and any , there exist constants , such that

 (8) E[supu∈[0,T]Y(u)p∣∣∣Ft]≤M1,p(1+Y(t)p),a.s.,

and

 (9) E[supu∈[0,T]Y(u)p]≤M2,p.
###### Proof.

The proof of (8) is an application of [4, Lemma 2.1] to conditional expectations requiring an invocation of Jensen’s inequality when . Eq. (9) is provided by [8, Lemma 3.2]. ∎

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

 (10) √|a1+a2| ≤ √|a1|+√|a2|; (11) |a1a2| ≤ 12(a21+a22); (12) (a1+…+an)p ≤ np(|a1|p+⋯+|an|p).

## 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) dX(t)=f(X(t))dt+g(X(t))dB(t),t∈[0,T],

over a random mesh on the interval given by,

 (14) Yn+1=Yn+hn+1f(Yn)+g(Yn)(W(tn+1)−W(tn)),

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) f(y)=αy+βy,g(y)=γ,

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

 Ftn={A∈F:A∩{tn≤t}∈Ft},n∈N.
###### Assumption 5.

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

 N=max{n∈N:tn−1

and the length of maximum and minimum stepsizes satisfy , and

 (16) hmin≤hn+1≤hmax≤1.
###### 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

 E[|W(tn+1)−W(tn)||Ftn] = 0,a.s.; E[|W(tn+1)−W(tn)|2|Ftn] = hn+1,a.s.; (17) E[∣∣∣∫stndW(r)∣∣∣p∣∣∣Ftn] = υp|s−tn|p2,a.s.

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).

Let be a solution of (14). We say that is a path-bounded time-stepping strategy for (14) if the conditions of Assumption 5 are satisfied and there exist real non-negative constants (where may be infinite if ) such that whenever ,

 (18) Q≤|Yn|

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.

Let be a solution of (14), and let be a time-stepping strategy that satisfies Assumption 5. If satisfies, for some ,

 (19) hn+1:=hmax×min{1,|Yn|r},n∈N,

or

 (20) hn+1=hmax×min(|Yn|r,|Yn|−r),

then it is path-bounding for (14).

###### Proof.

Suppose that (19) holds. When ,

 hn+1≤hmax|Yn|r⇔1|Yn|r≤hmaxhn+1≤hmaxhmin=ρ,n∈N,

and when it is obvious that , so we also have . Hence, when using the strategy defined by (19),

 |Yn|≥1ρ1/r,1|Yn|≤ρ1/r,n∈N,

so that (18) holds with and .

We can similarly show that (20) is path-bounding for (14), with and . ∎

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,

 f(Y2n)=α|Yn|2+β|Yn|2≤αρ1/r+β|Yn|2,

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) Yn+1=Yn+hn+1(αYn+βYn)+γΔWn+1.

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

 (22) ˜Y(s)=Yn+∫stn(αYn+βYn)dr+γ∫stndW(r),

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.

Define the map such that

 θ(y,z,h):=y+h(αy+βy)+γz,

so that, if is defined by (21), then

 Yn+1=θ(Yn,ΔWn+1,hn+1), n∈N.

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

 (23) ¯Y(s)=θ(¯Yn,∫stndW(r),∫stndr)I{hmin0}+φ(¯Yn,∫stndW(r),∫stndr)I{hn+1=hmin}∩{Yn+1>0}+φ(¯Yn,∫stndW(r),∫stndr)I{hmin

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

 (24) E[∣∣φ(¯Yn,ΔWn+1,hmin)−Y(tn+hmin)∣∣p∣∣Ftn]−∣∣¯Yn−Y(tn)∣∣p≤C1∫stnE[|¯Y(r)−Y(r)|p|Ftn]dr+C2hp/4+1min, n∈N, a.s.,

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

 (25) φ(y,ΔWn+1,h)>0a.s.,y∈R+,h∈[hmin,hmax].
###### 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

 E[|Y(s)−Y(tn)|p∣∣Ftn]≤2pγpυp|s−tn|p/2+¯Ln,p|s−tn|p, a.s.,

where

 ¯Ln,p:=22p(αpC(p,T)Y(tn)p+|β|pM1,p(1+Y(tn)p))

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

 Y(s)=Y(tn)+∫stn(αY(u)+βY(u))du+∫stnγdW(u),tn≤s≤T,

and therefore

 Y(s)−Y(tn)=∫stn(αY(u)+βY(u))du+γ(W(s)−W(tn)),tn≤s≤T.

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

 |Y(s)−Y(tn)|p ≤ 2p∣∣ ∣∣∫stn(αY(u)+βY(u))du∣∣ ∣∣p+2pγp|W(s)−W(tn)|p ≤ 2p|s−tn|p−1∫stn∣∣∣αY(u)+βY(u)∣∣∣pdu+2pγp|W(s)−W(tn)|p ≤ 22p|s−tn|p−1(∫stnαpY(u)pdu+∫stn|β|pY(u)pdu) +2pγp|W(s)−W(tn)|p,

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

 E[|Y(s)−Y(tn)|p∣∣Ftn]≤2pγpυp|s−tn|p/2+22p(αpC(p,T)Y(tn)p+βM1,p(1+Y(tn)p))|s−tn|p, a.s,

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) f(Y(s))=f(Y(tn))+Rf(s,tn,Y(tn)),s∈[tn,tn+1],

where

 (27)

For any , the conditional moment of

satisfies

 (28) E[|Rf|p∣∣Ftn]≤Kn,php/2n+1,a.s.,

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

Moreover, there exists independent of such that

 (29) Kp:=E[Kn,p]<∞.
###### Proof.

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

 E[|Rf|p∣∣Ftn]=E[∣∣ ∣∣(Y(s)−Y(tn))(β−αY(s)Y(tn))∣∣ ∣∣p∣∣∣Ftn].

Using the triangle inequality and (12) we get

 E[|Rf|p∣∣Ftn]≤2p|β|pE[|Y(s)−Y(tn)|p∣∣∣Ftn]+2pαpY(tn)pE[∣∣ ∣∣((Y(s)−Y(tn))⋅1Y(s))∣∣ ∣∣p∣∣∣Ftn]

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

 E[|Rf|p∣∣Ftn] ≤ 2p|β|php/2n+1(2pγpυp+¯Ln,php/2n+1) +2pαpY(tn)p√E[|Y(s)−Y(tn)|2p|Ftn] ⎷E[1|Y(s)|2p∣∣∣Ftn] ≤ 2p|β|php/2n+1(2pγpυp+¯Ln,php/2n+1) +2pαp√C(p,T)Y(tn)2p√E[|Y(s)−Y(tn)|2p|Ftn]

Again applying Lemma 11 and the elementary inequality (10) this becomes

 E[|Rf|p∣∣Ftn]≤2p|β|php/2n+1(2pγpυp+¯Ln,php/2n+1) +2php/2n+1αp√C(p,T)Y(tn)2p(2pγpυ1/22p+¯L1/2n,2php/2n+1),

from which the statement of the Lemma follows when we observe that the a.s. finiteness of is ensured by Assumption 1, and (29) is ensured by Lemmas 2 & 3. ∎

###### Lemma 13.

Let be the solution of (3) with initial value . Let be a solution of (22) over the interval and be a sequence of random timesteps defined by (3) and with . Then for , there exists an -measurable random variable with finite expectation such that

 (30) E[E(tn+1)2∣∣Ftn]−E(tn)2≤∫tn+1tnE[E(r)2|Ftn]+¯Knh2n+1,a.s.

where the error , .

###### Proof.

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

 (31) E(s) = Y(s)−˜Y(s) = [Y(tn)+∫stnf(Y(r))dr+γ∫stndW(r)] −[Yn+∫stnf(Yn)dr+γ∫stndW(r)] = E(tn)+∫stn~f(Y(r),Yn)dr,

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

 E(tn+1)2=E(tn)2+2∫tn+1tnE(r)~f(Y(r),Yn)dr.

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

 (32) E(tn+1)2−E(tn)2 = 2∫tn+1tnE(r)Rf(r,tn,Y(tn))dr+2∫tn+1tnE(r)~f(Y(tn),Yn)dr ≤ ∫tn+1tnE(r)2dr+∫tn+1tnRf(r,tn,Y(tn))2dr +2∫tn+1tnE(r)~f(Y(tn),Yn)dr.

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

 (33) ∫tn+1tnE(r)~f(Y(tn),Yn)dr = E(tn)~f(Y(tn),Yn)∫tn+1tndr+~f(Y(tn),Yn)2∫tn+1tn∫rtndudr +~f(