Hybrid, adaptive, and positivity preserving numerical methods for the Cox-Ingersoll-Ross model
aa r X i v : . [ q -f i n . C P ] F e b HYBRID, ADAPTIVE, AND POSITIVITY PRESERVINGNUMERICAL METHODS FOR THECOX-INGERSOLL-ROSS MODEL
C ´ONALL KELLY, GABRIEL LORD, AND HERU MAULANA
Abstract.
We introduce an adaptive Euler method for the approxi-mate solution of the Cox-Ingersoll-Ross short rate model. An explicitdiscretisation is applied over an adaptive mesh to the stochastic differ-ential equation (SDE) governing the square root of the solution, relyingupon a class of path-bounded timestepping strategies which work by re-ducing the stepsize as solutions approach a neighbourhood of zero. Themethod is hybrid in the sense that a backstop method is invoked if thetimestep becomes too small, or to prevent solutions from overshootingzero and becoming negative. Under parameter constraints that implyFeller’s condition, we prove that such a scheme is strongly convergent, oforder at least 1 /
2. Under Feller’s condition we also prove that the prob-ability of ever needing the backstop method to prevent a negative valuecan be made arbitrarily small. Numerically, we compare this adaptivemethod to fixed step schemes extant in the literature, both implicit andexplicit, and a novel semi-implicit adaptive variant. We observe thatthe adaptive approach leads to methods that are competitive over theentire domain of Feller’s condition. Introduction
The Cox-Ingersoll-Ross (CIR) process is a short rate model used in thepricing of interest rate derivatives and is given by the following Itˆo-typestochastic differential equation (SDE),(1) dX ( t ) = κ ( λ − X ( t )) dt + σ p X ( t ) dW ( t ) , t ∈ [0 , T ]; X (0) = X > , where W ( t ) is a Wiener Process and κ, λ, and σ are positive parameters,and T ∈ [0 , ¯ T ] for some fixed ¯ T < ∞ . Solutions of (1) are almost surely(a.s.) non-negative: in general they can achieve a value of zero but will bereflected back into the positive half of the real line immediately. Moreover, if2 κλ ≥ σ , 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 X ( t ) has,conditional upon X ( s ) for 0 ≤ s < t , a non-central chi-square distributionwith lim t →∞ E [ X ( t )] = λ and lim t →∞ Var[ X ( t )] = λσ / κ .For Monte Carlo estimation, exact sampling from the known conditionaldistribution of X ( t ) is possible but computationally inefficient and poten-tially restrictive if the innovating Wiener process of (1) is correlated withthat of another process: see [1,5,9,19]. Consequently a substantial literature Date : February 25, 2020.1991
Mathematics Subject Classification.
Key words and phrases.
Cox-Ingersoll-Ross model; Adaptive timestepping; ExplicitEuler-Maruyama method; Strong convergence; Positivity.
Method g ( x ) g ( x ) g ( x )Explicit Euler x x x Absorption x + x + x + Deelstra & Delbaen [7] x x x + Lord et al [19] x x + x + Reflection | x | | x | | x | Higham & Mao [10] x x | x | Diop [4] | x | x x Table 1.
Explicit Euler-Maruyama variants, adapted fromLord et al [19].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 ofthe explicit Euler-Maruyama method leads to schemes of the form(2) V n +1 = g ( V n + ∆ t ( κ ( λ − g ( V n )) + σ p g ( V n ))∆ W n for given functions g , g and g . These functions are selected to ensure thatthe diffusion coefficient remains real-valued (so that (2) is well defined), andin some cases to preserve the positivity of solutions. This approach seeksto accommodate the non-Lipschitz (square-root) diffusion, which facilitatesovershoot when the solutions are close to zero, but it introduces additionalbias to the approximation. A survey of choices common in practice maybe found in [19], and we summarise them in Table 1 using the convention x + := max { , x } .An alternative approach is to transform (1) before discretisation to makethe diffusion coefficient globally Lipschitz. For example, applying the Lam-perti transform Y = √ X yields an auxiliary SDE in Y with a state in-dependent and therefore globally Lipschitz diffusion, but a drift coefficientthat is unbounded when solutions are in a neighbourhood of zero. Thisapproach is effective: a fully implicit Euler discretisation over a uniformmesh that preserves positivity of solutions was proposed in [1] and shownto have uniformly bounded moments. A continuous time extension interpo-lating linearly between mesh points was shown to have strong L p order ofconvergence 1 / p | log( h ) | ) in [8] when 2 κλ > pσ , and acontinuous-time variant based on the same implicit discretisation was shownto have strong L p convergence of order one when 4 κλ > pσ in [2].In this article we will show that a strongly convergent numerical schemecan be constructed by an application of the Lamperti transform to (1) fol-lowed by an explicit Euler-Maruyama discretisation over a procedurally gen-erated adaptive mesh. The purpose of the adaptivity is to manage the non-linear drift response of the discrete tranformed system, a framework forwhich was introduced in [13] for SDEs with one-sided Lipschitz drift andglobally Lipschitz diffusion and extended to allow for monotone coefficientsand a Milstein-type discretisation in [14, 15] respectively. This frameworkimposes maximum and minimum timesteps h max and h min in a fixed ratio ρ and requires the use of a backstop numerical method in the event that the DAPTIVE TIMESTEPPING FOR COX-INGERSOLL-ROSS 3 timestepping strategy attempts to select a stepsize below h min . As in [15],we will use here path-bounded strategies, this time designed to increasethe density of timesteps when solutions approach zero, and we additionallyrequire the backstop method to retake a step when the adaptive strategyovershoots zero. This latter is carried out without discarding samples fromthe Brownian path (preserving the trajectory), and without bridging (pre-serving efficiency).We prove, under parameter constraints that imply the Feller condition(specifically, κλ > σ ), that the order of strong convergence in L is atleast 1 /
2. This parameter constraint is technical, and ensures sufficientlymany finite conditional inverse moments of solutions of (1) (as describedby [4]). We separately prove that, under exactly the Feller condition, theprobability of invoking the backstop method to avoid negative values canbe made arbitrarily small by choosing h max sufficiently small, and provide apractical method for doing so given a user defined tolerance level. The proofrelies upon a finite partitioning of the sample space of trajectories inducedby h max and h min , which allows us to handle the randomness of the numberof timesteps via the Law of Total Probability.Numerically we compare the convergence and efficiency of our hybridadaptive method with a semi-implicit adaptive variant, the fixed-step ex-plicit method due to [10], and the transformed implicit fixed-step methodproposed and analysed in [1, 2, 8], examining the parameter dependence ofthe numerical order of convergence in each case. The numerical conver-gence rates of adaptive methods are seen to outperform those of fixed-stepmethods over the entire domain where Feller’s condition holds.The structure of the article is as follows. In Section 2 we give the formof the SDE governing the Lamperti transform of (1), specify the constraintsplaced upon the parameters for the main strong convergence theorem, andexamine the availability of conditional moment and inverse moment boundsunder these constraints. In Section 3 we set up the framework for our randommesh, characterise the class of path-bounded timestepping strategies anddefine our adaptive numerical method. In Section 4 we present the twomain theorems on strong convergence and positivity, providing illustrativeexamples in the latter case. Finally in Section 5 we numerically compareconvergence and efficiency of several commonly used methods.2. Mathematical Preliminaries
Throughout this article we let ( F t ) t ≥ be the natural filtration of W . Byusing Itˆo’s Formula and applying the transformation Y = √ X we get, dY ( t ) = (cid:18) κ ( λ − X t ) − σ √ X t (cid:19) dt + σ dW t , t ∈ [0 , T ]; Y (0) = p X ∈ R + . By then setting, α = 4 κλ − σ , β = − κ , γ = σ , we can write,(3) dY ( t ) = (cid:18) αY ( t ) + βY ( t ) (cid:19) dt + γdW t , t ∈ [0 , T ]; Y (0) = p X ∈ R + , C ´ONALL KELLY, GABRIEL LORD, AND HERU MAULANA where f ( y ) = α/y + βy is not globally Lipschitz continuous, but when α > β < f ( x ) − f ( y )]( x − y ) ≤ β ( x − y ) , for all x, y ∈ R + , which can be seen by noting that f ( x ) − f ( y ) = ( x − y ) (cid:20) β − αxy (cid:21) . Meanwhile the diffusion coefficient g ( y ) = γ is constant and is thereforeglobally Lipschitz continuous. The SDE (3) has integral form(4) Y ( t ) = Y (0) + Z t (cid:18) αY ( s ) + βY ( s ) (cid:19) ds + Z t γdW t , t ≥ . In order to ensure the a.s. positivity of solutions of (3) and the bound-edness of certain inverse moments of solutions of (3), we will also need thefollowing assumption:
Assumption 1.
Suppose that (5) κλ > σ . Eq. (5) implies the Feller Condition (2 κλ ≥ σ ; see, for example [20]),which ensures that solutions of (1), and therefore (3), remain positive withprobability one: P [ Y ( t ) > , t ≥
0] = 1 . Assumption 1 provides inverse moment bounds as follows:
Lemma 2.
Let ( Y ( t )) t ∈ [0 ,T ] be a solution of (3) , where Assumption 1 holds,and let ≤ t < s ≤ T . For any Y (0) > , and for ≤ p ≤ , there exists C ( p, T ) > such that (6) E (cid:20) Y ( s ) p (cid:12)(cid:12)(cid:12)(cid:12) F t (cid:21) ≤ C ( p, T ) Y ( t ) p , a.s. Proof.
Let ( X ( t )) t ∈ [0 ,T ] be a solution of (1) where Assumption 1 holds. ByLemma A.1 in Bossy & Diop [4],(7) E (cid:20) X ( t ) (cid:21) ≤ e κt X and E (cid:20) X ( t ) p (cid:21) ≤ C ( p, T ) X p , for some C ( p, T ) and any p such that 1 < p < κλσ −
1. Assumption 1 ensuresthat κλσ − >
3, and since Y ( t ) = p X ( t ), (6) follows by Lemma A.1 in [4]as it applies to conditional expectations, the former requiring an additionalapplication of Jensen’s inequality to the first inequality in (7). (cid:3) 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 ( Y ( t )) t ∈ [0 ,T ] be a solution of (3) , where Assumption 1 holds,and let ≤ t ≤ T . For any Y (0) > and any p > , there exist constants M ,p , M ,p < ∞ , such that (8) E " sup u ∈ [0 ,T ] Y ( u ) p (cid:12)(cid:12)(cid:12)(cid:12) F t ≤ M ,p (1 + Y ( t ) p ) , a.s., DAPTIVE TIMESTEPPING FOR COX-INGERSOLL-ROSS 5 and (9) E " sup u ∈ [0 ,T ] Y ( u ) p ≤ M ,p . Proof.
The proof of (8) is an application of [4, Lemma 2.1] to conditionalexpectations requiring an invocation of Jensen’s inequality when 0 < p < (cid:3)
Finally, we will make frequent use of the following elementary inequalities:for n ∈ N and a , . . . , a n ∈ R and p ≥ p | a + a | ≤ p | a | + p | a | ;(10) | a a | ≤
12 ( a + a );(11) ( a + . . . + a n ) p ≤ n p ( | a | p + · · · + | a n | p ) . (12) 3. An Adaptive Numerical Method [13] provided a framework within which to construct timestepping strate-gies for an adaptive explicit Euler-Maruyama numerical scheme applied tononlinear Itˆo-type SDEs of the form(13) dX ( t ) = f ( X ( t )) dt + g ( X ( t )) dB ( t ) , t ∈ [0 , T ] , over a random mesh { t n } n ∈ N on the interval [0 , T ] given by,(14) Y n +1 = Y n + h n +1 f ( Y n ) + g ( Y n )( W ( t n +1 ) − W ( t n )) , where { h n } n ∈ N is a sequence of random timesteps and { t n = P ni =1 h i } Nn =1 with t = 0. The random time step h n +1 is determined by Y n .For completeness, we present here the essential elements of that frame-work, before discussing the dynamical considerations specific to the trans-formed 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 { t n } n ∈ N is an F t -stopping time: i.e. { t n ≤ t } ∈ F t , for all t ≥ , where ( F t ) t ≥ isthe natural filtration of W . We may then define a discrete time filtration {F t n } n ∈ N by F t n = { A ∈ F : A ∩ { t n ≤ t } ∈ F t } , n ∈ N . Assumption 5.
Each h n is F t n − -measurable, and N is a random integersuch that, N = max { n ∈ N : t n − < T } and t N = T, and the length of maximum and minimum stepsizes satisfy h max = ρh min , <ρ < ∞ , and (16) h min ≤ h n +1 ≤ h max ≤ . C ´ONALL KELLY, GABRIEL LORD, AND HERU MAULANA
Remark 6.
The lower bound h min ensures that a simulation over the in-terval [0 , T ] can be completed in a finite number of timesteps, and the upperbound h max prevents stepsizes from becoming too large. The latter is used asa convergence parameter in our examination of the strong convergence of theadaptive method. The random variable N cannot take values outside the fi-nite set { N min , . . . , N max } , where N min := ⌊ T /h max ⌋ and N max := ⌈ T /h min ⌉ . △ W n +1 := W ( t n +1 − W ( t n ) is a Wiener increment over a random in-terval the length of which depends on Y n , through which it depends on { W ( s ) , s ∈ [0 , t n ] } . Therefore △ W n +1 is not independent of F t n ; indeedit is not necessarily normally distributed. Since h n +1 is a bounded F t n -stopping time and F t -measurable, then W ( t n +1 ) − W ( t n ) is F t n -conditionallynormally distributed, by Doob’s optional sampling theorem (see for exam-ple [22]), and for all p ≥ υ p < ∞ such that E [ | W ( t n +1 ) − W ( t n ) ||F t n ] = 0 , a.s. ; E [ | W ( t n +1 ) − W ( t n ) | |F t n ] = h n +1 , a.s. ; E (cid:20)(cid:12)(cid:12)(cid:12)(cid:12)Z st n dW ( r ) (cid:12)(cid:12)(cid:12)(cid:12) p (cid:12)(cid:12)(cid:12)(cid:12) F t n (cid:21) = υ p | s − t n | p , a.s. (17)3.2. Adaptive timestepping strategy.
To ensure strong convergence,our strategy will be to reduce the size of each timestep if discretised so-lutions attempt to enter a neighbourhood of zero. If we wish to control thelikelihood of invoking the backstop to avoid negative values, we will alsoreduce the timestep when solutions grow large.
Definition 7 (A path-bounded time-stepping strategy) . Let { Y n } n ∈ N bea solution of (14) . We say that { h n } n ∈ N is a path-bounded time-steppingstrategy for (14) if the conditions of Assumption 5 are satisfied and thereexist real non-negative constants ≤ Q < R (where R may be infinite if Q = 0 ) such that whenever h min < h n < h max , Q ≤ | Y n | < R, n = 0 , . . . , N − . (18)We now give two examples of path-bounded strategies that are valid for(14), the first with R infinite (which is sufficient to ensure strong conver-gence), and the second with R finite (which is useful if we also wish tominimise the use of the backstop to ensure positivity). Lemma 8.
Let { Y n } n ∈ N be a solution of (14) , and let { h n } n ∈ N be a time-stepping strategy that satisfies Assumption 5. If { h n } n ∈ N satisfies, for some r ≥ , (19) h n +1 := h max × min { , | Y n | r } , n ∈ N , or (20) h n +1 = h max × min( | Y n | r , | Y n | − r ) , then it is path-bounding for (14) .Proof. Suppose that (19) holds. When | Y n | < h n +1 ≤ h max | Y n | r ⇔ | Y n | r ≤ h max h n +1 ≤ h max h min = ρ, n ∈ N , DAPTIVE TIMESTEPPING FOR COX-INGERSOLL-ROSS 7 and when | Y n | ≥ | Y n | ≤
1, so we also have | Y n | r ≤ ρ < ∞ .Hence, when using the strategy defined by (19), | Y n | ≥ ρ /r , | Y n | ≤ ρ /r , n ∈ N , so that (18) holds with Q = 1 /ρ /r and R = ∞ .We can similarly show that (20) is path-bounding for (14), with Q =1 /ρ /r and R = ρ /r . (cid:3) Note that for the strategies defined by (19) and (20), solutions of (14)cannot enter the neigbourhood ( ± /ρ /r ), and therefore terms of the se-quence (1 / | Y n | ) n ∈ N are uniformly bounded from above. This has the effectof controlling inverse moments of the solutions of (14). Moreover for (19),when (15) holds, f ( Y n ) = α | Y n | + β | Y n | ≤ αρ /r + β | Y n | , and therefore that strategy is admissible in the sense of [13, Definition 2.2]with R = αρ /r and R = β . Similarly, (20) is admissible with R = αρ /r + βρ − /r and R = 0.3.3. The adaptive numerical method with backstop.
We consider anadaptive scheme based upon the following explicit Euler-Maruyama discreti-sation of (3) over a random mesh given by,(21) Y n +1 = Y n + h n +1 (cid:18) αY n + βY n (cid:19) + γ ∆ W n +1 . where the timestep sequence is constructed according to (19). For s ∈ [ t n , t n +1 ), the continuous version is given by(22) e Y ( s ) = Y n + Z st n (cid:18) αY n + βY n (cid:19) dr + γ Z st n dW ( r ) , so that e Y ( t n ) = Y n for each n ∈ N .We combine this scheme with a positivity-preserving backstop schemethat is to be applied if the timestepping strategy attempts to select atimestep below h min (in which case we choose h n +1 = h min ) or if the currentselected timestep and subsequently observed Brownian increment △ W n +1 would result in the approximation becoming negative: Definition 9.
Define the map θ : R → R such that θ ( y, z, h ) := y + h (cid:18) αy + βy (cid:19) + γz, so that, if { Y n } n ∈ N is defined by (21) , then Y n +1 = θ ( Y n , ∆ W n +1 , h n +1 ) , n ∈ N . C ´ONALL KELLY, GABRIEL LORD, AND HERU MAULANA
Then we define the continuous form of an adaptive explicit Euler schemeto be the sequence of functions { ( ¯ Y ( s )) s ∈ [ t n ,t n +1 ) } n ∈ N obeying (23) ¯ Y ( s ) = θ (cid:18) ¯ Y n , Z st n dW ( r ) , Z st n dr (cid:19) I { h min
Since the events { Y n +1 < } and { Y n +1 > } are F t n +1 -measurable but not F t n -measurable, if a negative value of Y n +1 is observedfollowing a step of length h n +1 we must retake the step using the backstopmethod, which will ensure positivity over that step by (25) . This introducesan element of backtracking into the algorithm, but as long as the originallycomputed stepsize h n +1 and Brownian increment are retained we can stay onthe same trajectory while avoiding the use of a Brownian bridge. Theorem15 in Section 4.2, illustrated by Example 16, demonstrates that it is alwayspossible to choose h max to ensure that this particular use of the backstop canbe avoided with probability − ε , for arbitrarily small ε ∈ (0 , , on eachtrajectory. 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 timesteppingstrategy. Second, we investigate the likelihood that the adaptive part of themethod generates a negative value (triggering the use of the backstop to en-sure positivity) and show how h max may be chosen to control the probabilityof this occurring.4.1. Strong convergence of the adaptive method with path-boundedtimestepping strategy.Lemma 11.
Let ( Y ( t )) t ∈ [0 ,T ] be a solution of (3) and let { t n } n ∈ N be arandom mesh such that each t n is an F t -stopping time. Fix n ∈ N andsuppose that t n ≤ s ≤ T , where T ∈ [0 , ¯ T ] . Then, for any ≤ p ≤ , wehave E (cid:2) | Y ( s ) − Y ( t n ) | p (cid:12)(cid:12) F t n (cid:3) ≤ p γ p υ p | s − t n | p/ + ¯ L n,p | s − t n | p , a.s., DAPTIVE TIMESTEPPING FOR COX-INGERSOLL-ROSS 9 where ¯ L n,p := 2 p (cid:18) α p C ( p, T ) Y ( t n ) p + | β | p M ,p (1 + Y ( t n ) p ) (cid:19) is an F t n -measurable random variable with finite expectation, and C ( T ) , M ,p are the constants defined by (6) and (8) in the statements of Lemmas2 and 3 respectively .Proof. Solutions of (3) satisfy the integral equation Y ( s ) = Y ( t n ) + Z st n (cid:18) αY ( u ) + βY ( u ) (cid:19) du + Z st n γdW ( u ) , t n ≤ s ≤ T, and therefore Y ( s ) − Y ( t n ) = Z st n (cid:18) αY ( u ) + βY ( u ) (cid:19) du + γ ( W ( s ) − W ( t n )) , t n ≤ s ≤ T. Using the triangle and Cauchy-Schwarz inequalities, and the elementaryinequality (12) with n = 2, | Y ( s ) − Y ( t n ) | p ≤ p (cid:12)(cid:12)(cid:12)(cid:12)Z st n (cid:18) αY ( u ) + βY ( u ) (cid:19) du (cid:12)(cid:12)(cid:12)(cid:12) p + 2 p γ p | W ( s ) − W ( t n ) | p ≤ p | s − t n | p − Z st n (cid:12)(cid:12)(cid:12)(cid:12) αY ( u ) + βY ( u ) (cid:12)(cid:12)(cid:12)(cid:12) p du + 2 p γ p | W ( s ) − W ( t n ) | p ≤ p | s − t n | p − (cid:18)Z st n α p Y ( u ) p du + Z st n | β | p Y ( u ) p du (cid:19) +2 p γ p | W ( s ) − W ( t n ) | p , for s ∈ [ t n , T ]. Now apply conditional expectations on both sides withrespect to F t n and (17) to get, E (cid:2) | Y ( s ) − Y ( t n ) | p (cid:12)(cid:12) F t n (cid:3) ≤ p γ p E (cid:2) | W ( s ) − W ( t n ) | p (cid:12)(cid:12) F t n (cid:3) + 2 p | s − t n | p − (cid:18) E (cid:20)Z st n α p Y ( u ) p du (cid:12)(cid:12)(cid:12)(cid:12) F t n (cid:21) + E (cid:20)Z st n | β | p Y ( u ) p du (cid:12)(cid:12)(cid:12)(cid:12) F t n (cid:21)(cid:19) ≤ p γ p υ p | s − t n | p/ + 2 p | s − t n | p − (cid:18) α p Z st n C ( p, T ) Y ( t n ) p du + | β | p Z st n M ,p (1 + Y ( t n ) p ) du (cid:19) , a.s, where we have used (6) and (8) from the statement of Lemma 2 at the laststep. Therefore E (cid:2) | Y ( s ) − Y ( t n ) | p (cid:12)(cid:12) F t n (cid:3) ≤ p γ p υ p | s − t n | p/ + 2 p (cid:18) α p C ( p, T ) Y ( t n ) p + βM ,p (1 + Y ( t n ) p ) (cid:19) | s − t n | p , a.s, as required. (cid:3) Lemma 12.
Let ( Y ( t )) t ∈ [0 ,T ] be a solution of (3) and let Assumption 1 hold.Let { t n } n ∈ N arise from the adaptive timestepping strategy satisfying (18) in Definition 7 for some < Q < R , and formulate the Taylor expansion of f ( Y ( s )) around Y ( t n ) , where f is as given in (15) , as (26) f ( Y ( s )) = f ( Y ( t n )) + R f ( s, t n , Y ( t n )) , s ∈ [ t n , t n +1 ] , where (27) R f ( s, t n , Y ( t n )) = Z Df ( Y ( t n ) + τ ( Y ( s ) − Y ( t n ))) ( Y ( s ) − Y ( t n )) dτ. For any ≤ p ≤ , the p th conditional moment of R f ( s, t n , Y ( t n )) satisfies (28) E (cid:2) | R f | p (cid:12)(cid:12) F t n (cid:3) ≤ K n,p h p/ n +1 , a.s., where K n,p is an a.s. finite and F t n -measurable random variable given by K n,p = 2 p | β | p (cid:16) p γ p υ p + ¯ L n,p h p/ n +1 (cid:17) +2 p α p p C ( p, T ) Y ( t n ) p (cid:16) p γ p υ / p + ¯ L / n, p h p/ n +1 (cid:17) Moreover, there exists K p independent of n such that (29) K p := E [ K n,p ] < ∞ . Proof.
By direct substitution of f ( y ) from (15) into (27), and taking the p th -moment conditional upon F t n , we get E (cid:2) | R f | p (cid:12)(cid:12) F t n (cid:3) = E (cid:20)(cid:12)(cid:12)(cid:12)(cid:12) ( Y ( s ) − Y ( t n )) (cid:18) β − αY ( s ) Y ( t n ) (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) p (cid:12)(cid:12)(cid:12)(cid:12) F t n (cid:21) . Using the triangle inequality and (12) we get E (cid:2) | R f | p (cid:12)(cid:12) F t n (cid:3) ≤ p | β | p E (cid:20) | Y ( s ) − Y ( t n ) | p (cid:12)(cid:12)(cid:12)(cid:12) F t n (cid:21) + 2 p α p Y ( t n ) p E (cid:20)(cid:12)(cid:12)(cid:12)(cid:12)(cid:18) ( Y ( s ) − Y ( t n )) · Y ( s ) (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) p (cid:12)(cid:12)(cid:12)(cid:12) F t n (cid:21) Next apply Lemma 11 followed by the Cauchy-Schwarz inequality to get E [ | R f | p (cid:12)(cid:12) F t n ] ≤ p | β | p h p/ n +1 (cid:16) p γ p υ p + ¯ L n,p h p/ n +1 (cid:17) +2 p α p Y ( t n ) p p E [ | Y ( s ) − Y ( t n ) | p |F t n ] s E (cid:20) | Y ( s ) | p (cid:12)(cid:12)(cid:12)(cid:12) F t n (cid:21) ≤ p | β | p h p/ n +1 (cid:16) p γ p υ p + ¯ L n,p h p/ n +1 (cid:17) +2 p α p p C ( p, T ) Y ( t n ) p p E [ | Y ( s ) − Y ( t n ) | p |F t n ]Again applying Lemma 11 and the elementary inequality (10) this be-comes E [ | R f | p (cid:12)(cid:12) F t n ] ≤ p | β | p h p/ n +1 (cid:16) p γ p υ p + ¯ L n,p h p/ n +1 (cid:17) +2 p h p/ n +1 α p p C ( p, T ) Y ( t n ) p (cid:16) p γ p υ / p + ¯ L / n, p h p/ n +1 (cid:17) , DAPTIVE TIMESTEPPING FOR COX-INGERSOLL-ROSS 11 from which the statement of the Lemma follows when we observe that thea.s. finiteness of K n,p is ensured by Assumption 1, and (29) is ensured byLemmas 2 & 3. (cid:3) Lemma 13.
Let ( Y ( t n )) t n ∈ [0 ,T ] be the solution of (3) with initial value Y (0) = Y = √ X . Let (cid:16) e Y ( s ) (cid:17) s ∈ [ t n ,t n +1 ] be a solution of (22) over theinterval [ t n , t n +1 ] and { h n } n ∈ N be a sequence of random timesteps defined by (3) and { t n = P ni =1 h i } Nn =1 with t = 0 . Then for n, p ∈ N , there exists an F t n -measurable random variable ¯ K n with finite expectation K n := E [ ¯ K n ] < ∞ such that (30) E (cid:2) E ( t n +1 ) (cid:12)(cid:12) F t n (cid:3) − E ( t n ) ≤ Z t n +1 t n E (cid:2) E ( r ) |F t n (cid:3) + ¯ K n h n +1 , a.s. where the error E ( s ) := Y ( s ) − e Y ( s ) , s ∈ [ t n , t n +1 ] .Proof. For s ≥ t n , we subtract (22) from (4) to get E ( s ) = Y ( s ) − e Y ( s )= (cid:20) Y ( t n ) + Z st n f ( Y ( r )) dr + γ Z st n dW ( r ) (cid:21) − (cid:20) Y n + Z st n f ( Y n ) dr + γ Z st n dW ( r ) (cid:21) = E ( t n ) + Z st n ˜ f ( Y ( r ) , Y n ) dr, (31)where f is defined as in (15) and ˜ f ( Y ( r ) , Y n ) = f ( Y ( r )) − f ( Y n ). Applyingthe Itˆo formula and setting s = t n +1 , we can write, E ( t n +1 ) = E ( t n ) + 2 Z t n +1 t n E ( r ) ˜ f ( Y ( r ) , Y n ) dr. By (26) in the statement of Lemma 12, ˜ f ( Y ( r ) , Y n ) = ˜ f ( Y ( t n ) , Y n )+ R f ( r, t n , Y ( t n ),where R f is defined in (27). This, and an application of (11) gives E ( t n +1 ) − E ( t n ) = 2 Z t n +1 t n E ( r ) R f ( r, t n , Y ( t n )) dr + 2 Z t n +1 t n E ( r ) ˜ f ( Y ( t n ) , Y n ) dr ≤ Z t n +1 t n E ( r ) dr + Z t n +1 t n R f ( r, t n , Y ( t n )) dr +2 Z t n +1 t n E ( r ) ˜ f ( Y ( t n ) , Y n ) dr. (32) Consider the third term on the RHS of (32), and substitute (31) into theintegrand: Z t n +1 t n E ( r ) ˜ f ( Y ( t n ) , Y n ) dr = E ( t n ) ˜ f ( Y ( t n ) , Y n ) Z t n +1 t n dr + ˜ f ( Y ( t n ) , Y n ) Z t n +1 t n Z rt n du dr + ˜ f ( Y ( t n ) , Y n ) Z t n +1 t n Z rt n R f ( u, t n , Y ( t n )) du dr ≤ β Z t n +1 t n E ( t n ) dr + ˜ f ( Y ( t n ) , Y n ) h n +1 + ˜ f ( Y ( t n ) , Y n ) Z t n +1 t n Z rt n R f ( u, t n , Y ( t n )) du dr ≤ ˜ f ( Y ( t n ) , Y n ) h n +1 + ˜ f ( Y ( t n ) , Y n ) Z t n +1 t n Z rt n R f ( u, t n , Y ( t n )) du dr, (33)Now substitute (33) into (32), to get(34) E ( t n +1 ) − E ( t n ) ≤ Z t n +1 t n E ( r ) dr + Z t n +1 t n R f ( r, t n , Y ( t n )) dr + 2 ˜ f ( Y ( t n ) , Y n ) h n +1 + 2 ˜ f ( Y ( t n ) , Y n ) Z t n +1 t n Z rt n R f ( u, t n , Y ( t n )) du dr. Apply expectations to both sides of (34), conditional upon F t n , to get E (cid:2) E ( t n +1 ) |F t n (cid:3) − E ( t n ) ≤ Z t n +1 t n E [ E ( r ) |F t n ] dr + Z t n +1 t n E [ R f ( r, t n , Y ( t n )) |F t n ] dr + 2 ˜ f ( Y ( t n ) , Y n ) h n +1 + 2 ˜ f ( Y ( t n ) , Y n ) Z t n +1 t n Z rt n E [ R f ( u, t n , Y ( t n )) |F t n ] du dr, a.s. Apply the bound (28) in the statement of Lemma 12 with p = 1 , E (cid:2) E ( t n +1 ) |F t n (cid:3) − E ( t n ) ≤ Z t n +1 t n E [ E ( r ) |F t n ] dr + (2 ˜ f ( Y ( t n ) , Y n ) + K n, ) h n +1 + 2 ˜ f ( Y ( t n ) , Y n ) K n, h / n +1 , a.s. Since, by (16) in the statement of Assumption 5, h n +1 ≤ h max ≤
1, thestatement of the Lemma now follows, with ¯ K n := 2 ˜ f ( Y ( t n ) , Y n ) + K n, +2 ˜ f ( Y ( t n ) , Y n ) K n, .To see that E [ ¯ K n ] < ∞ , note that by (18) in Definition 7, and since β < f ( Y n ) ≤ α/Q . The finiteness of E [ K n, ] and E [ K n, ] is given by (6)in the statement of Lemma 2 with p = 6 , (cid:3) DAPTIVE TIMESTEPPING FOR COX-INGERSOLL-ROSS 13
Theorem 14.
Let ( Y ( t )) t ∈ [0 ,T ] be the solution of (3) with initial value Y (0) = Y = √ X . Let ( ¯ Y ( t )) t ∈ [0 ,T ] be a solution of (23) with initial value ¯ Y (0) = Y (0) and path-bounded timestepping strategy { h n } n ∈ N satisfying theconditions of Definition 7 for some < Q < R , with R possibly infinite.There exists C > , independent of h max , such that E [ | Y ( T ) − ¯ Y ( T ) | ] ≤ Ch max . Proof.
From (30) in the statement of Lemma 13, we have that when h n +1 ≥ h min , E (cid:2) E ( t n +1 ) |F t n (cid:3) − E ( t n ) ≤ Z t n +1 t n E (cid:2) E ( r ) |F t n (cid:3) dr + ¯ K n h n +1 , a.s. Suppose that h n +1 < h min and Y n +1 is generated from Y n via an applicationof the backstop method over a single step of length h min . This correspondsto single application of the map ϕ in Definition 9 and therefore the relation(24) is satisfied.Define the positive constant Γ = C ∨ F t n -measurable randomvariables ¯Γ n, = C ∨ ¯ K n , for n ∈ N . Noting again that, by (16) in thestatement of Assumption 5, h n +1 ≤ h max ≤
1, we see that the combinedmethod given by (23) satisfies, on almost all trajectories,(35) E (cid:2) E ( t n +1 ) |F t n (cid:3) − E ( t n ) ≤ Γ Z t n +1 t n E (cid:2) E ( r ) |F t n (cid:3) dr + ¯Γ n, h n +1 . Sum both sides of (35) over n = 0 , . . . , N − E " N − X n =0 (cid:0) E (cid:2) E ( t n +1 ) |F t n (cid:3) − E ( t n ) (cid:1) ≤ Γ E " N − X n =0 Z t n +1 t n E (cid:2) E ( r ) |F t n (cid:3) dr + E " N − X n =0 ¯Γ n, h n +1 . Consider first the LHS of (36). Since N ≤ N max , we can write E " N − X n =0 (cid:0) E (cid:2) E ( t n +1 ) |F t n (cid:3) − E ( t n ) (cid:1) = N max − X n =0 E (cid:2)(cid:0) E (cid:2) E ( t n +1 ) |F t n (cid:3) − E ( t n ) (cid:1) I { N ≥ n +1 } (cid:3) = N max − X n =0 (cid:0) E (cid:2) E (cid:2) E ( t n +1 ) |F t n (cid:3) (cid:12)(cid:12) σ ( N ≥ n + 1) (cid:3) − E (cid:2) E ( t n ) (cid:12)(cid:12) σ ( N ≥ n + 1) (cid:3)(cid:1) P [ N ≥ n + 1] , (37)a.s., where σ ( N ≥ n +1) is the σ -algebra generated by the event { N ≥ n +1 } .Since N is a F t n -stopping time, the event { N ≤ n } ∈ F t n and therefore sois { N ≤ n } c = { N ≥ n + 1 } , and therefore σ ( N ≥ n + 1) ⊆ F t n . By the tower property of conditional expectations we get E " N − X n =0 (cid:0) E (cid:2) E ( t n +1 ) |F t n (cid:3) − E ( t n ) (cid:1) = N max − X n =0 (cid:0) E (cid:2) E ( t n +1 ) (cid:12)(cid:12) N ≥ n + 1 (cid:3) − E (cid:2) E ( t n ) (cid:12)(cid:12) N ≥ n + 1 (cid:3)(cid:1) P [ N ≥ n + 1]= N max − X n =0 (cid:0) E (cid:2) E ( t n +1 ) I { N ≥ n +1 } (cid:3) − E (cid:2) E ( t n ) I { N ≥ n +1 } (cid:3)(cid:1) = E (cid:2) E ( t N ) (cid:3) = E (cid:2) E ( T ) (cid:3) . Similarly, the RHS of (36) can be writtenΓ E " N − X n =0 Z t n +1 t n E (cid:2) E ( r ) |F t n (cid:3) dr + N − X n =0 ¯Γ n, h n +1 = Γ E " N max − X n =0 Z t n +1 t n E (cid:2) E ( r ) |F t n (cid:3) I { N ≥ n +1 } dr ( I ) + E " N max − X n =0 ¯Γ n, h n +1 I { N ≥ n +1 } ( II ) . Consider first (I). Since t n , t n +1 , and I { N ≥ n +1 } are F t n -measurable wecan bring I { N ≥ n +1 } inside the conditional expectation, and exchange theorder of integration and conditional expectation as follows E " N max − X n =0 Z t n +1 t n E (cid:2) E ( r ) |F t n (cid:3) I { N ≥ n +1 } dr = N max − X n =0 E (cid:20)Z t n +1 t n E (cid:2) E ( r ) |F t n (cid:3) I { N ≥ n +1 } dr (cid:21) = N max − X n =0 E (cid:20)Z t n +1 t n E (cid:2) E ( r ) I { N ≥ n +1 } |F t n (cid:3) dr (cid:21) = N max − X n =0 E (cid:20) E (cid:20)Z t n +1 t n E ( r ) I { N ≥ n +1 } dr (cid:12)(cid:12)(cid:12)(cid:12) F t n (cid:21)(cid:21) = N max − X n =0 E (cid:20)Z t n +1 t n E ( r ) I { N ≥ n +1 } dr (cid:21) = E " N max − X n =0 Z t n +1 t n E ( r ) I { N ≥ n +1 } dr = E " N max − X n =0 Z t n +1 t n E ( r ) I { r ≤ T } dr = E (cid:20)Z T E ( r ) dr (cid:21) = Z T E (cid:2) E ( r ) (cid:3) dr. (38) DAPTIVE TIMESTEPPING FOR COX-INGERSOLL-ROSS 15
Finally consider (II). We have, since E (cid:2) ¯Γ n, I { N ≥ n +1 } (cid:3) ≤ E (cid:2) ¯Γ n, (cid:3) ≤ C ∨ K =: Γ for all n = 0 , . . . , N max − N max = ⌈ T /h min ⌉ , and by Assumption16, ρh min = h max ≤ E " N max − X n =0 ¯Γ n, h n +1 I { N ≥ n +1 } ≤ E " h N max − X n =0 ¯Γ n, I { N ≥ n +1 } = h N max − X n =0 E (cid:2) ¯Γ n, I { N ≥ n +1 } (cid:3) ≤ h Γ N max ≤ h Γ (cid:18) Th min + 1 (cid:19) ≤ ( ρT + 1)Γ h max . (39)Substituting (37), (38), and (39) back into (36) we get E [ E ( T ) ] ≤ Γ Z T E [ E ( r ) ] dr + ( ρT + 1)Γ h max . Since this inequality holds if T is varied continuously over [0 , ¯ T ], for any¯ T < ∞ , an application of Gronwall’s inequality gives the result. (cid:3) On the positivity of an adaptive method with path-boundedstrategy.
In this section we assume Feller’s condition (2 κλ ≥ σ ) to en-sure that solutions of (1) remain a.s. positive, but we do not require thatAssumption 1 holds.4.2.1. Probability of positivity over a single step.
Consider the timesteppingstrategy defined by (19) with r = 1, satisfying Definition 7 with R = ∞ .The probability of solutions of (21) becoming negative after a single stepwith this strategy, and hence triggering a use of the backstop method, isgiven by P [ Y k +1 < | Y k = y >
0] = Φ ( a ( y )) , y > h min , where Φ( x ) = √ π R x −∞ e − s / ds , and a ( y ) = − y − (cid:16) αy + βy (cid:17) h max (1 ∧ y r ) γ p h max (1 ∧ y r ) . Figure 1 presents two surface plots of these one-step probabilities against y and h max . Feller’s condition is satisfied in both cases. In Figure 1 (a), whenAssumption 1 holds, we see that the probability of invoking the backstopto avoid a negative value is highest when h max is large and Y n is close toor above 1, in which case the timestepping strategy will tend to select h n +1 to be close to h max . This probability drops off rapidly as h max reduces. InFigure 1 (b), when Assumption 1 does not hold, the highest probabilitiesof invoking the backstop for preserving positivity when Y n is close to h min .However the maximum probability is significantly lower than any seen inFigure 1 (a). κλ > σ κλ < σ (a): σ = 0 . , λ = 0 . , κ = 2 (b): σ = 0 . , λ = 0 . , κ = 1 Figure 1.
Surface plots of probabilities of solutions of (14)with (15) and timestepping strategy (19) becoming negativeover a single step, for h max ∈ [0 . ,
1] and Y = y ∈ [ h min , . ρ = 2 . In (a) Assumption 1 holds. In (b)
Assump-tion 1 does not hold.4.2.2.
Probability of positivity over a full trajectory.
If we require path-bounded strategies where
R < ∞ , it is possible to derive an upper limiton h max that is sufficient, over the entire trajectory, to keep the probabil-ity of needing the backstop scheme to prevent a negative value below somearbitrarily small tolerance. Our analysis reworks and extends the approachtaken in the proof of [16, Theorem 4.3], using adaptive timestepping to han-dle unboundedness in the drift term. By allowing the use of the backstop toensure a minimum timestep, along with an application of the Law of TotalProbability, we can avoid fixing the random number of steps N . Theorem 15.
Let { Y n } Nn =0 be a solution of (23) , with initial value Y > ,evaluated on a random mesh { h n } Nn =1 satisfying the conditions of Definition7 with R < ∞ . Suppose also that Y ∈ (0 , R ) . Then, for each ε ∈ (0 , there exists ¯ h max ( ε ) > such that, for all h max ∈ (0 , ¯ h max ( ε )) P [ R N ] > − ε. where R N := T Nj =0 { Y j > } .Proof. Since, by (25), the backstop method will ensure positivity over asingle step if h n +1 = h min , the event { Y n +1 > } is equivalent to the followingunion: (( △ W n +1 p h n +1 > − γ Y n h n +1 + α p h n +1 Y n + β p h n +1 Y n !) ∩ { h n +1 > h min } ) ∪ { h n +1 = h min } . DAPTIVE TIMESTEPPING FOR COX-INGERSOLL-ROSS 17
Moreover, when Y n > (cid:26) u ∈ − γ (cid:18) Qh max + α √ h min R + β p h max R (cid:19)(cid:27) ⊆ ( u ∈ − γ Y n h n +1 + α p h n +1 Y n + β p h n +1 Y n !) . For each i = N min , . . . N max , define Ω i := { ω ∈ Ω : N ( ω ) = i } , so that { Ω } N max i = N min is a finite partition of the sample space Ω. On each Ω i define thesequence of sub-events R n ( i ) = { Y n > , Y n − > , . . . , Y > , Y > } ∩ Ω i , n = 0 , , . . . i. Recall that the random variable △ W n +1 / p h n +1 is distributed condi-tionally upon F t n like a standard Normal random variable. Moreover, R n ( i ) ∈ F t n for n = 0 , . . . , i and i = N min , . . . , N max . Let Φ denote thedistribution function of a standard Normal random variable, and suppose { ξ n } n ∈ N is a sequence of mutually independent standard normal randomvariables. Then P [ Y n +1 > |R n ( i )]= P [ { Y n +1 > } ∩ {R n ( i ) } ] P [ R n ( i )] = E (cid:2) E (cid:2) I { Y n +1 > }∩{R n ( i ) } |F t n (cid:3)(cid:3) P [ R n ( i )] ≥ E (cid:2) E (cid:2) I { Y n +1 > }∩{R n ( i ) } |F t n (cid:3)(cid:3) = E [ P [ { Y n +1 > } ∩ {R n ( i ) }|F t n ]] ≥ E [ P [ h n +1 = h min |F t n ]]+ E " P " △ W n +1 p h n +1 > − γ (cid:18) Qh max + α √ h min R + β p h max R (cid:19) (cid:12)(cid:12)(cid:12)(cid:12) F t n ≥ E " P " △ W n +1 p h n +1 > − γ (cid:18) Qh max + α √ h min R + β p h max R (cid:19) (cid:12)(cid:12)(cid:12)(cid:12) F t n = P (cid:20) ξ n +1 > − γ (cid:18) Qh max + α √ h max R √ ρ + β p h max R (cid:19)(cid:21) = 1 − Φ (cid:18) − γ (cid:18) Qh max + p h max (cid:18) αR √ ρ + βR (cid:19)(cid:19)(cid:19) = Φ (cid:18) γ (cid:18) Qh max + p h max (cid:18) αR √ ρ + βR (cid:19)(cid:19)(cid:19) , n = 0 , . . . , i − . Since Y > P [ R ( i ) | Ω i ] = 1 and therefore, since R n ( i ) ⊆ Ω i , Φtakes values on [0 , i ≤ N max , P [ R i ( i ) | Ω i ] ≥ P [ R i ( i )] = P " i \ n =0 R n ( i ) = i Y n =1 P [ R n ( i ) |R n − ( i ) , . . . , R ( i )]= i − Y n =0 P [ Y n +1 > |R n ( i )] ≥ Φ (cid:18) γ (cid:18) Qh max + p h max (cid:18) αR √ ρ + βR (cid:19)(cid:19)(cid:19) i ≥ Φ (cid:18) γ (cid:18) Qh max + p h max (cid:18) αR √ ρ + βR (cid:19)(cid:19)(cid:19) N max , for i = N min , . . . , N max . Multiplying through by P [Ω i ] and applying theLaw of Total Probability by summing both sides over i = N min , . . . , N max gives P [ R N ] = N max X i = N min P [ R i ( i ) | Ω i ] P [Ω i ] ≥ N max X i = N min Φ (cid:18) γ (cid:18) Qh max + p h max (cid:18) αR √ ρ + βR (cid:19)(cid:19)(cid:19) N max P [Ω i ]= Φ (cid:18) γ (cid:18) Qh max + p h max (cid:18) αR √ ρ + βR (cid:19)(cid:19)(cid:19) N max . Fix ε ∈ (0 , h max ∈ (0 , ¯ h max ( ε )), we have(41) Φ (cid:18) γ (cid:18) Qh max + p h max (cid:18) αR √ ρ + βR (cid:19)(cid:19)(cid:19) N max ≥ − ε. To (41), apply the following inequality due to [21]1 √ π Z x − x e − s / ds > p − e − x / , x ∈ R + , along with the fact that N max = ρT /h max , leading us to seek h max so that
12 + 12 vuuuut − exp − (cid:16) Qh max + √ h max (cid:16) αR √ ρ + βR (cid:17)(cid:17) γ ρTh max ≥ − ε. DAPTIVE TIMESTEPPING FOR COX-INGERSOLL-ROSS 19
Thus we derive the bound˜ h max ( ε ) := sup ( ¯ h ∈ (0 ,
1) : Qh + √ h (cid:18) αR √ ρ + βR (cid:19) ≥ r ln (cid:16) − (2(1 − ε ) hρT − (cid:17) − γ , h ∈ (0 , ¯ h ) ) . ¯ h max ( ε ) is uniquely defined for each ε ∈ (0 ,
1) because(42) g ( h ) := Qh + √ h (cid:18) αR √ ρ + βR (cid:19) − q − γ ln (cid:0) − (2(1 − ε ) h/ρT − (cid:1) is continuous on R + with lim h → + g ( h ) = ∞ , and therefore there is aneighbourhood of zero corresponding to (0 , ¯ h max ( ε )) within which g is posi-tive. (cid:3) Note that if we extend the interval of simulation [0 , T ] and keep ε ∈ (0 , N max in (41). This will leadto a reduction in the bound ¯ h max ( ε ), in a way that is characterised by (42).More generally g ( h ), as defined by (42) in the proof of Theorem 15, providesa practical guide for choosing h max in order to control the probability ofinvoking the backstop to avoid negative values. Example 16.
Consider two adaptive timestepping strategies based on (20) with r = 1 and ρ = 2 , , each used to simulate a single trajectory of (3) over the interval [0 , using the adaptive method (23) . In each case, we wishto choose h max so that the probability of requiring the backstop in order toavoid negative values on that trajectory is less than ε , and this will hold forany Y ∈ ( h min , R ) .Table 2 shows the value of ¯ h ( ε ) for a range of tolerances ε for param-eter sets where Assumption 1 is satisfied ( κλ > σ ), and where it is not( κλ < σ ). The resulting bounds on h max are determined by substituting allparameters into (42) and solving g ( h ) = 0 for h using the fsolve commandin Maple with 20 digits of precision. We report the first 4 significant digitsin each case, which is sufficient to illustrate the sensitivity of these boundsto the choice of ρ and ε . Numerical simulation
Given (1) and its associated transformation (3), we compare our hybridadaptive method (23), referred to in this section as
Explicit Adaptive(EA) , to a natural semi-implicit variant constructed by replacing the updateequation (21) with(43) Y n +1 = (1 − βh n +1 ) − (cid:20) Y n + h n +1 αY n + γ △ W n +1 (cid:21) , referred to in this section as Semi-Implicit Adaptive (SIA) . In both caseswe will use the adaptive timestepping strategy given by (19) with r = 1 (notethat we see similar results when r = 2). We also compare to two fixed stepmethods: the explicit discretisation of (1) analysed by [10] given by X n +1 = X n + hκ ( λ − X n ) + σ p | X n |△ W n +1 , ρ = 2 , Q = 0 . , R = 64 ε ¯ h max ( ε ) σ = 0 . − . × − λ = 0 .
05 10 − . × − κ = 2 10 − . × − ε ¯ h max ( ε ) σ = 0 . − . × − λ = 0 .
05 10 − . × − κ = 1 10 − . × − ρ = 2 , Q = 0 . , R = 256 ε ¯ h max ( ε ) σ = 0 . − . × − λ = 0 .
05 10 − . × − κ = 2 10 − . × − ε ¯ h max ( ε ) σ = 0 . − . × − λ = 0 .
05 10 − . × − κ = 1 10 − . × − Table 2.
Bounds on h max ensuring positivity (without theuse of the backstop) of trajectories of (24) with probabilityat least 1 − ε where the path-bounded timestepping strategysatisfies Definition 7. Assumption 1 is satisfied ( κλ > σ )for tables in the left column, and violated ( κλ < σ ) fortables in the right column.referred to in this section as Explicit Fixed (EF) , and the drift implicitsquare root discretisation of (3) proposed and analysed in [1, 2, 8], given by Y n +1 = Y n + γ △ W n +1 − βh ) + s ( Y n + γ △ W n +1 ) − βh ) + αh − βh , and referred to in this section as Implicit Fixed (IF) .In the first part, we will compare the strong convergence of these methodsin the mean stepsize, and the corresponding numerical efficiency. In the sec-ond part we look at strong convergence in h max , and explore the dependenceon model parameters.5.1. Strong convergence and efficiency.
Throughout the section, wetake ρ = 2 . We solve using EA and SIA with values of h max = 2 − i , i = 4 , . . . , M sample trajectories to estimate p E [ | X ( T ) − X N | ],the root mean square error (RMSE), at a final time T = 1. To computeerror estimates we first generate a reference solution using IF over a meshwith stepsize h = 2 − , using a Brownian bridge to ensure values for theadaptive approximations are on the reference trajectory. To ensure thatwe are comparing adaptive and fixed step schemes of similar average cost,when solving using IF and EF we take as the fixed step h mean the averageof all timesteps h ( m ) n taken by EA over each path and each realisation ω m , m = 1 , . . . , M so that h mean = 1 M M X m =1 N ( ω m ) N ( ωm ) X n =1 h ( ω m ) n . In Figure 2 we examine strong convergence for these methods by plottingRMSE against h mean with M = 1000 on a log-log scale, and efficiency by DAPTIVE TIMESTEPPING FOR COX-INGERSOLL-ROSS 21 -3 -2 hmean10 -5 -4 -3 R M SE EA EF IF SIA -2 -1 CPU Time (Second)10 -4 -3 R M SE EA EF IF SIA (a) (b) -3 -2 hmean10 -5 -4 -3 R M SE EA EF IF SIA -2 -1 CPU Time (Second)10 -4 -3 R M SE EA EF IF SIA (c) (d)
Figure 2.
Convergence and efficiency of methods applied to(1) with λ = 0 . σ = 0 . Y = 0 . κ = 2 in (a) and (b) ,and κ = 1 in (c) and (d) .plotting RMSE against average compute time (cputime) again with M =1000.In Figure 2 (a), Assumption 1 holds, and the estimated error at each valueof h mean is comparable for all methods, and the numerical order appearsto be close to one. In Figure 2 (b) we also see comparable efficiencies asmeasured by CPU time for this example. In Figure 2 (c), Assumption 1 doesnot hold, and we see first that the numerical order of EF has reduced, andthe estimated error at each value of h mean is lowest for EA and SIA , whichalso demonstrate the fastest CPU times in Figure 2 (d) for lower RMSEvalues.5.2.
Parameter dependence of the strong convergence rate.
Finally,we investigate numerically the dependence of the rate of strong convergenceon the value of the parameter a := σ / (2 κλ ). Note that Assumption 1 cor-responds to a < .
25, and Feller’s condition corresponds to a ≤
1. For 30uniformly spaced values of a in the interval [0 . , . L as the slope of the corresponding a R a t e EA EF IF SIA a R a t e EA EF IF SIA (a) (b)
Figure 3.
Estimation of strong convergence order for meth-ods applied to (1) with λ = 0 . Y = 0 . κ = 2 in (a) ,and κ = 0 . (b) . Here, a = σ / (2 κλ ). In both cases, As-sumption 1 holds to the left of the vertical line at a = 0 . a = 1.error over a range of values of h mean computed by generated strong conver-gence plots as in Figure 2 and using the polyfit command in MATLAB toestimate the order of strong convergence for each method.In Figure 3 we observe that EF maintains the highest numerical order ofconvergence: at or close to one while Feller’s condition holds. The reductionin order outside of this region, which is visible for all methods, occurs moresharply in Figure 3 (b) when κ is small. This is followed by SIA , whichmaintains a numerical order of convergence close to EA when Assumption 1holds, but reduces more quickly outside this region. The difference is morepronounced in Figure 3 (a), and this may be because updates using the SIA method are subject to a damping factor (1 − βh n +1 ) − , as can be seen in(43), which has greater effect for larger values of κ . Finally, we observe for EF an uptick in convergence rate as a approaches zero, and this is consistentwith the notion that it should display order one convergence in the absenceof noise.Figure 4 demonstrates how frequently the backstop was invoked in theproduction of Figure 3, where we separately track usage to ensure positivityand usage to bound below the stepsize at h min . We see that when Assump-tion 1 holds, we do not require the backstop to avoid negative values, thoughas we move to the boundary of that region we start to use it to bound thestepsize at h min for a small proportion of steps. Note also that the usage toavoid negative values increases with a when κ = 2 only, whereas usage tobound the timestep increases with a for both κ = 2 , .
2, but more rapidlyin the latter case.
DAPTIVE TIMESTEPPING FOR COX-INGERSOLL-ROSS 23 % B a cks t op EA =2 SIA =2 EA =0.2 SIA =0.2 % B a cks t op EA =2 SIA =2 EA =0.2 SIA =0.2 (a) : Avoid Y n < (b) : Ensure h n = h min Figure 4.
Percentage of times the backstop was invoked inthe production of Figure 3 in order to (a) avoid a negativevalue or (b) bound the timestep from below by h min . Acknowledgement(s)
The authors are grateful to Professor Alexandra Rodkina, of the Univer-sity of the West Indies at Mona, Jamaica, and Ms Fandi Sun, of Heriot-WattUniversity, Edinburgh, UK, for useful discussion in preparing this work.
Funding
The first and third authors were supported by a grant from
LembagaPengelola Dana Pendidikan (LPDP) Republik Indonesia . References [1] Alfonsi, A.
On the discretization schemes for the CIR (and Bessel squared) processes ,Monte Carlo Methods and Applications. (2005), pp. 355-384.[2] Alfonsi, A.[3] Berkaoui, A., Bossy, M. and Diop, A. BEuler Scheme for SDEs with Non-LipschitzDiffusion Coefficient: Strong Convergence , ESAIM Probability and Statistics, :1(2008), pp. 1-11.[4] Bossy, M. and Diop, A. An Efficient Discretization Scheme for One DimensionalSDEs with a Diffusion Coefficient Function of the Form | x | α , α ∈ [ , Exact Simulation of Stochastic Volatility and Other AffineJump Diffusion Processes , Oper. Res. 54, 217-231. 2006.[6] Cox, J. C., Ingersoll, J. E., and Ross, S. A.
A Theory of the Term Structure of InterestRates , Econometrica, :2 (1985), pp 385-407.[7] Deelstra, G. and Delbaen, F. Convergence of Discretized Stochastic (Interest Rate)Process with Stochastic Drift Term , Applied Stochastic Models and Data Analysis, :1 (1998), pp. 77-84.[8] Dereich, S., Neuenkirch, A., and Szpruch, L.. An Euler-Type Method for The StrongApproximation of The Cox-Ingersoll-Ross Process , Proceedings of The Royal SocietyA Mathematical Physical and Engineering Sciences (2012), pp. 1105-1115.[9] Glassermann, P.
Monte Carlo Methods in Financial Engineering , Springer, New York,2008. [10] Higham, D. J., and Mao, X.
Convergence of the Monte Carlo Simulations Involvingthe Mean Reverting Square Root Process , Journal of Computational Finance, :3(2005), pp. 35-62.[11] Hutzenthaler, M., Jentzen, A., and Kloeden, P. E. Strong Convergence of An ExplicitNumerical Method for SDEs with Non-Globally Lipschitz Continuous Coefficients ,Annals of Applied Probability, :4 (2012), pp. 1611-1641.[12] Jafari, M. A. and Abbasian, S. The Moments for Solution of the Cox-Ingersoll-RossInterest Rate Model , Journal of Finance and Economics, :1 (2017), pp. 34-37.[13] Kelly, C. and Lord, G. J. Adaptive Timestepping Strategies for Nonlinear StochasticSystems , IMA Journal of Numerical Analysis, :3 (2018), pp. 1523-1549.[14] Kelly, C. and Lord, G. J. Adaptive Euler Methods for Stochastic Systems with Non-Globally Lipschitz Coefficients , arXiv:1805.11137 (2018), 21 pages.[15] Kelly, C., Lord, G. J., and Sun, F.
Strong convergence of an adaptive time-steppingMilstein method for SDEs with one-sided Lipschitz drift , arXiv:1909.00099 (2019), 20pages.[16] Kelly, C., Rodkina, A., and Rapoo, E. M.
Adaptive timestepping for pathwise sta-bility and positivity of strongly discretised nonlinear stochastic differential equations ,Journal of Computational and Applied Mathematics, (2018), pp. 39-57.[17] Kloeden, P. E. and Platen, E.
Numerical Solution of Stochastic Differential Equations ,Stochastic Modeling and Applied Probability. Springer Berlin Heidelberg, 2011.[18] Liu, W. and Mao, X.
Almost sure stability of the Euler-Maruyama method with ran-dom variable stepsize for stochastic differential equations , Numer. Algorithms, :2(2017), pp. 573–592.[19] Lord, R., and Koekkoek R., and Van Dijk, D. A comparison of biased simulationschemes for stochastic volatility models , Quantitative Finance, :2 (2010), pp. 177-194.[20] Mao, X. Stochastic Differential Equations and Applications(2nd Edition) , WoodheadPublishing UK, 2007.[21] Sasvari, Z. and Chen, H. Tight Bounds for the Normal Distribution.
The AmericanMathematical Monthly , :1 (1999), pp. 76.[22] Shiryaev, A. N. Probability (2nd edition), Springer, Berlin, 1996.[23] Stuart, A. M. and Humphries, A. R. Dynamical Systems and Numerical Analysis .Cambridge University Press, 1996.[24] Zeytun, S. and Gupta, A.
A Comparative Study of the Vasicek and the CIR Model ofthe Short Rate , Berichte des Fraunhofer ITWM, Kaiserslautern, Germany, 2007.
E-mail address : [email protected] School of Mathematical Sciences, University College Cork, Republic ofIreland
E-mail address : [email protected] Department of Mathematics, Radboud University, Netherlands
E-mail address : [email protected]@fmipa.unp.ac.id