# Approximating Explicitly the Mean-Reverting CEV Process

Authors: N. Halidias; I. S. Stamatiou
Journal: Journal of Probability and Statistics
 (2015)

---

## Abstract

We are interested in the numerical solution of mean-reverting CEV processes that appear in financial mathematics models and are described as nonnegative solutions of certain stochastic differential equations with sublinear diffusion coefficients of the form(xt)q, where 1/2<q<1. Our goal is to construct explicit numerical schemes that preserve positivity. We prove convergence of the proposed SD scheme with rate depending on the parameter q. Furthermore, we verify our findings through numerical experiments and compare with other positivity preserving schemes. Finally, we show how to treat the two-dimensional stochastic volatility model with instantaneous variance process given by the above mean-reverting CEV process.


---

## Body

## 1. Introduction

Consider the following stochastic models in Itô form:(1)St=S0+∫0tμ·Sudu+∫0tVup·SudWu,t∈0,T,Vt=V0+∫0tk1-k2Vsds+∫0tk3VsqdW~st∈0,T,where St represents the underlying financially observable variable, Vt is the instantaneous volatility when p=1 or the instantaneous variance when p=1/2, and the Wiener processes (Wt), (W~t) have correlation ρ.We assume that(Vt) is a mean-reverting CEV process of the above form, with the coefficients ki>0 for i=1,2,3 and q>1/2, since the process (Vt) has to be nonnegative. To be more precise the above restriction on q implies that (Vt) is positive; that is, 0 is unattainable, as well as nonexplosive; that is, ∞ is unattainable, as can be verified by Feller’s classification of boundaries [17, Proposition 5.22] (see also Appendix A). The steady-state level of Vt is k1/k2 and the rate of mean-reversion is k2.System (1) for p=q=1/2 is the Heston model. When q=1 we get the Brennan-Schwartz model [2, Section II] that despite its simple form cannot provide analytical expressions for St.Process(Vt) for q=1/2, also known as the CIR process [6, (13)], by the initials of the authors that proposed it for the term structure of interest rates, has received a lot of attention and we just mention the latest two contributions to the study of such processes (see [4, 5] and references therein).Process(Vt) for 1/2≤q≤1 has been also considered for the dynamics of the short-term interest rate [3, (1)]. The stationary distribution of the process has also been derived in [7, Proposition 2.2].We aim for a positivity preserving scheme for the process(Vt). The scheme that we propose and denote as semidiscrete (SD) preserves the analytical property of (Vt) staying positive. The idea of the semidiscrete method is that we discretize a part of the original SDE and then apply Itô’s formula (cf. [8] where the method originally appeared and [5, 9, 23]). The explicit Euler scheme fails to preserve positivity, as well as the standard Milstein scheme. We intend to apply the semidiscrete method for the numerical approximation of (Vt) in model (1) with 1/2<q<1 and compare with other positivity preserving methods such as the balanced implicit method (BIM) (introduced by [11, (3.2)] with the positivity preserving property [25, Section 5] and its stability properties [13]; see also [14] for an extended balanced method with better stability behavior) and the balanced Milstein method (BMM) [25, Theorem 5.9] (we give in Appendix B the form of all the above schemes for the approximation of (Vt)). Finally, we approximate the stochastic volatility model (1) with p=1/2. In [15] a thorough treatment can be found, where also another stochastic volatility model is suggested.Section2 provides the setting and the main results, Theorems 1 and 2, concerning the L2-convergence of the proposed semidiscrete method to the true solution of mean-reverting CEV processes of the form of the stochastic volatility in (1). The rate of mean-square convergence in Theorem 1 is logarithmic and in Theorem 2 is polynomial with magnitude 1/2(q-1/2). The main ingredient of the approach we adopt, inspired by [16], is a change of the initial Brownian motion (Wt) to another Brownian motion (W^t) justified by Lévy’s martingale characterization of Brownian motion.Section3 is devoted to the logarithmic rate of convergence of the proposed semidiscrete scheme, while Section 4 concerns the proof of the polynomial rate of convergence. In Section 5 we briefly discuss the case where we do not alter the initial Brownian motion (Wt). This approach produces reduced convergence rate. Finally, Section 6 presents illustrative figures where the behavior of the proposed scheme, regarding the order of convergence, is shown and a comparison with BIM and BMM schemes is given. In Section 7 we treat the full model (1) for a special case. Concluding remarks are in Section 8 and in Appendix B we briefly present numerical schemes for the integration of the variance-volatility process (Vt).


## 2. The Setting and the Main Results

We consider the following SDE:(2)xt=x0+∫0tk1-k2xsds+∫0tk3xsqdWs,t∈0,T,where k1,k2,k3 are positive and 1/2<q<1. Then, Feller’s test implies that there is a unique strong solution such that xt>0 a.s. when x0>0 a.s. Let (3)fθx,y=k1-k21-θx-k32/41+k2θΔx2q-1-k2θy︸f1x,y+k32/41+k2θΔx2q-1︸f2x,(4)gx,y=k3xq-1/2y, where f(x,x)=a(x)=k1-k2x and g(x,x)=b(x)=k3xq.Let the partition0=t0<t1<⋯<tN=T with Δ=T/N and consider the following process: (5)ytSDq=ytn+f1ytn,yt·Δ+∫tntf2ytnds+∫tntsgn⁡zsgytn,ysdWs,with y0=x0 a.s. or more explicitly (6)ytSDq=ytn+k1-k21-θytn-k3241+k2θΔytn2q-1-k2θyt·Δ+∫tntk3241+k2θΔytn2q-1ds+k3ytnq-1/2∫tntsgn⁡zsysdWs,for t∈(tn,tn+1], where θ∈[0,1] represents the level of implicitness and (7)zt=yn+k321+k2θΔytnq-1/2Wt-Wtn,with (8)yn≔ytn1-k2Δ1+k2θΔ+k1Δ1+k2θΔ-k3241+k2θΔ2ytn2q-1Δ.Process (6) is well defined when yn≥0 and this is true when 1/(1+k2θΔ)(k3)2≤4(k2∧k1) and Δ(2-θ)≤1/k2. Furthermore, (6) has jumps at nodes tn. Solving for yt, we end up with the following explicit scheme: (9)ytSDq=yn+∫tntk3241+k2θΔ2ytn2q-1ds+k31+k2θΔytnq-1/2∫tntsgn⁡zsysdWs, with solution in each step given by [1, (4.39), page 123](10)ytSDq=zt2, which has the pleasant feature ytSD(q)≥0.Inspired by [16] we remove the term sgn⁡(zs) from (6) by considering the process (11)W~t=∫0tsgn⁡zsdWs, which is a martingale with quadratic variation W~t,W~t=t and thus a standard Brownian motion with respect to its own filtration, justified by Lévy’s theorem [17, Theorem 3.16, page 157]. Therefore, the compact form of (6) becomes (12)ytSD=x0+∫0tk1-k21-θys^-k2θys~ds+∫ttn+1k1-k21-θytn-k3241+k2θΔytn2q-1-k2θytds+k3∫0tys^q-1/2ysdW~s, for t∈(tn,tn+1] where (13)s^=tj,s∈tj,tj+1,j=0,…,n,s~=tj+1,fors∈tj,tj+1,t,fors∈tn,tj=0,…,n-1.Consider also the process (14)x~t=x0+∫0tk1-k2x~sds+∫0tk3x~sqdW~s,t∈0,T. The process (xt) of (2) and the process (x~t) of (14) have the same distribution. We show in the following that Esup0≤t≤T|ytSD(q)-x~t|2→0 as Δ↓0; thus the same holds for the unique solution of (2); that is, Esup0≤t≤T|ytSD(q)-xt|2→0 as Δ↓0. To simplify notation we write W~,(x~t) as W,(xt). We end up with the following explicit scheme: (15)ytSDq=yn+∫tntk3241+k2θΔ2ytn2q-1ds+k31+k2θΔytnq-1/2∫tntysdWs, where yn is as in (8).Assumption A.
Let the parametersk1,k2,k3 be positive such that 1/1+k2θΔ(k3)2≤4(k2∧k1) and consider Δ>0 such that Δ2-θ<1/k2, for θ∈[0,1]. Moreover assume x0≥0 a.s. and E(x0)p<A for some p≥4.Theorem 1 (logarithmic rate of convergence).
Let Assumption A hold. The semidiscrete scheme (15) converges to the true solution of (2) in the mean-square sense with rate given by (16)Esup0≤t≤TytSDq-xt2≤Cln⁡Δ-1⁡, where C is independent of Δ and given by (17)C≔326ϵk34T2e6T2k22+k2T,where 0<ϵ<q-1/2.Assumption B.
Let Assumption A hold where nowx0∈R and x0>0.Theorem 2 (polynomial rate of convergence).
Let Assumption B hold. Then the semidiscrete scheme (15) converges to the true solution of (2) in the mean-square sense with rate given by (18)Esup0≤t≤TytSDq-xt2≤CΔq-1/2,where (19)C≔16k32TA2x0+k1T2A^4q-22e6k22T2+CHKϵ-1x01-qνλ, and CHK is the constant described in (83) and λ is an appropriately chosen positive parameter which satisfies (84) and always exists, ν(λ)≔λ/2(1-q)2(k3)2-1, and ϵ>1.In the following sections we write for simplicityytSD or yt for ytSD(q).


## 3. Logarithmic Rate of Convergence

### 3.1. Moment Bounds

Lemma 3 (moment bound for SD approximation).
It holds that(20)Esup0≤t≤Tytp≤ApEx0+k1Tp, for any p>2, where Ap≔expp(p-1)/2(k3)2p-1/2p+2p-1/pT.Proof of Lemma3.
We first observe that(yt) is bounded in the following way: (21)0≤yt≤x0+∫0tk1ds+k3∫0tys^q-1/2ysdWs≤x0+k1T+k3∫0tys^q-1/2ysdWs≔ut, a.s., where the lower bound comes from the construction of (yt) and the upper bound follows from a comparison theorem. We will bound (ut) and therefore (yt), since 0≤yt≤ut a.s. Set the stopping time τR≔inf⁡{t∈[0,T]:ut>R}, for R>0 with the convention inf⁡∅=∞. Application of Itô’s formula on (ut∧τR)p implies (22)ut∧τRp=x0+k1Tp+pp-12k32∫0t∧τRusp-2ys^2q-1ysds+pk3∫0t∧τRusp-1ys^q-1/2ysdWs≤x0+k1Tp+pp-12k32∫0t∧τRusp-1ys^2q-1ds+Mt≤x0+k1Tp+pp-12k32∫0t∧τRp-12pusp+2p-1pys^2q-1pds+Mt≤x0+k1Tp+pp-12k32p-12p+2p-1p∫0t∧τRuspds+Mt, where in the second step we have used the fact that 0≤yt≤ut, in the third step the inequality xp-1y≤ϵp-1/pxp+1/pϵp-1yp, valid for x∧y≥0 and p>1 with ϵ=1/2, and in the final step the fact that 1/2<q<1 and Mt≔pk3∫0t∧τR‍(us)p-1(ys^)q-1/2ysdWs. Taking expectations in the above inequality and using that Mt is a local martingale vanishing at 0, we get (23)Eut∧τRp≤Ex0+k1Tp+pp-12k32p-12p+2p-1p∫0tEus∧τRpds≤Ex0+k1Tpexp⁡pp-12k32p-12p+2p-1pT≤ApEx0+k1Tp, where we have applied the Gronwall inequality [18, (7)]. We have that (24)yt∧τRp=yτRpIτR≤t+ytpIt<τR≥ytpIt<τR. Thus taking expectations in the above inequality and using the estimated upper bound for E(ut∧τR)p we arrive at (25)EytpIt<τR≤Eyt∧τRp≤Eut∧τRp≤ApEx0+k1Tp,and, taking the limit as R→∞, we get (26)limR→∞⁡EytpIt<τR≤ApEx0+k1Tp. Let us fix t. The sequence of stopping times τR is increasing in R and t∧τR→t as R→∞, and thus the sequence (yt)pI{t<τR} is nondecreasing in R and (yt)pI{t<τR}→(yt)p as R→∞. Application of the monotone convergence theorem implies (27)Eytp≤ApEx0+k1Tp, for any p>2. Using again Itô’s formula on (ut)p, taking the supremum, and then using Doob’s martingale inequality on the diffusion term we bound Esup0≤t≤T(ut)p and thus Esup0≤t≤T(yt)p.Lemma 4 (error bound for SD scheme).
Letns be an integer such that s∈[tns,tns+1]. Then (28)Eys-ys^p≤A^pΔp/2,Eys-ys~p<A~pΔp/2, for any p>0, where the positive quantities A^p,A~p do not depend on Δ.Proof of Lemma4.
First we take ap>2. We get that(29)ys-ys^p=∫tnssk1-k21-θyu^-k2θyu~du+∫tnstns+1k2θys^du-∫stns+1k2θysdu+∫stnsk1-k21-θytns-k3241+k2θΔytns2q-1du+k3∫tnssyu^q-1/2yudWup≤5p-1∫tnssk1-k21-θyu^-k2θyu~dup+k2pθpys^ptns+1-tnsp+k2pθpysptns+1-sp+∫stnsk1-k21-θytns-k3241+k2θΔytns2q-1dup+k3p∫tnssyu^q-1/2yudWup≤5p-1tns-sp-1∫tnssk1-k21-θyu^-k2θyu~pdu+k2pθpys^p+yspΔp+k1-k21-θytns-k3241+k2θΔytns2q-1pΔp+k3p∫tnssyu^q-1/2yudWup, where we have used the Cauchy-Schwarz inequality. Taking expectations in the above inequality and using Lemma 3 and Doob’s martingale inequality on the diffusion term we conclude (30)Eys-ys^p≤A^pΔp/2, where the positive quantity A^p, except on p, depends also on the parameters k1,k2,k3,θ,q, but not on Δ. Now, for 0<p<2, we get (31)Eys-ys^p≤Eys-ys^2p/2≤A^pΔp/2, where we have used Jensen’s inequality for the concave function ϕ(x)=xp/2. Following the same lines, we can show that (32)Eys-ys~p≤A~pΔp/2, for any 0<p, where the positive quantity A~p, except on p, depends also on the parameters k1,k2,k3,θ,q, but not on Δ.For the rest of this section we rewrite again the compact form of (12) in the following way: (33)ytSD=x0+∫0tfθys^,ys~ds+∫0tgys^,ysdWs︸ht+∫ttn+1f1ytn,ytds,where fθ(·,·) is given by (3) and the auxiliary process (ht) is close to (yt) as shown in the next result.Lemma 5 (moment bounds involving the auxiliary process).
For anys∈[0,T] it holds that (34)Ehs-ysp≤CpΔp,Ehsp≤Ch,and for s∈[tn,tn+1] one has that (35)Ehs-ys^p≤C^pΔp/2,Ehs-ys~p≤C~pΔp/2,for any p>0, where the positive quantities Cp,C^p,C~p,Ch do not depend on Δ.Proof of Lemma5.
We have that(36)hs-ysp=∫stn+1f1ytn,ysdup≤tn+1-spf1ytn,ysp, for any p>0, where we have used (33). Using Lemma 3 we get the left part of (34). Now for p>2 and noting that (37)Ehsp≤2p-1Ehs-ysp+2p-1Eysp≤2p-1CpΔp+2p-1ApEx0+k1Tp≤Ch,we get the right part of (34), where we have used Lemma 3. The case 0<p<2 follows by Jensen’s inequality as in Lemma 4.
Furthermore, fors∈[tn,tn+1] and p>2, we derive that (38)Ehs-ys^p≤2p-1Ehs-ysp+2p-1Eys-ys^p≤2p-1CpΔp+2p-1A^pΔp/2≤C^pΔp/2, where we have used (30) and in the same manner (39)Ehs-ys~p≤2p-1CpΔp+2p-1A~pΔp/2≤C~pΔp/2. The case 0<p<2 follows by Jensen’s inequality.


### 3.2. Convergence of the Auxiliary Process(ht) to (xt) in L1

We will use the representation (33) and write (40)ht-xt=∫0tfθys^,ys~-fθxs,xsds+∫0tgys^,ys-gxs,xsdWs.Proposition 6.
Let Assumption A hold. Then one has(41)sup0≤t≤TEht-xt≤J3Δq-1/2mem+2k32T1mek2T, for any m>1, where em=e-m(m+1)/2 and (42)J3≔2k32TA2Ex0+k1T2A^4q-2.Proof of Proposition6.
Let the nonincreasing sequence{em}m∈N with em=e-m(m+1)/2 and e0=1. We introduce the following sequence of smooth approximations of |x| (method of Yamada and Watanabe [19]): (43)ϕmx=∫0xdy∫0yψmudu, where the existence of the continuous function ψm(u) with 0≤ψm(u)≤2/(mu) and support in (em,em-1) is justified by ∫emem-1‍(du/u)=m. The following relations hold for ϕm∈C2(R,R) with ϕm(0)=0,(44)x-em-1≤ϕmx≤x,ϕm′x≤1,x∈R,ϕm′′x≤2mx,whenem<x<em-1,ϕm′′x=0otherwise.We have that(45)Eht-xt≤em-1+Eϕmht-xt.Moreover we find that (46)fθys^,ys~-fθxs,xs=k1-k21-θys^-k2θys~-k1-k2xs=-k21-θys^-xs-k2θys~-xs=k21-θhs-ys^+k2θhs-ys~-k2hs-xs,(47)gys^,ys-gxs,xs2=k3ys^q-1/2ys-k3xsq2≤k32ysys^q-1/2-ysq-1/2+ysq-xsq2≤2k32ysys^q-1/2-ysq-1/22+ysq-xsq2≤2k32ysys^-ys2q-1+ys-xs2≤2k32ysys^-ys2q-1+hs-ys+hs-xs,where we have used properties of Hölder continuous functions and, namely, the fact that xq is q-Hölder continuous for q≤1, that is, |xq-yq|≤|x-y|q, and that xq is 1/2-Hölder continuous since q>1/2. Application of Itô’s formula to the sequence {ϕm}m∈N implies (48)ϕmht-xt=∫0tϕm′hs-xsfθys^,ys~-fθxs,xsds+Mt+12∫0tϕm′′hs-xsgys^,ys~-gxs,xs2ds≤∫0tk21-θhs-ys^+k2θhs-ys~+k2hs-xsds+Mt+2k32∫0t1mhs-xsysys^-ys2q-1+hs-ys+hs-xsds≤k21-θ∫0ths-ys^ds+k2θ∫0ths-ys~ds+2k32mem∫0ths-ysds+k2∫0ths-xsds+Mt+2k32mem∫0tysys^-ys2q-1ds+2k32Tm, where in the second step we have used (46) and (47) and the properties of ϕm and (49)Mt≔∫0tϕm′hu-xugyu^,yu~-gxu,xudWu.Taking expectations in the above inequality yields (50)Eϕmht-xt≤k21-θ∫0tEhs-ys^ds+k2θ∫0tEhs-ys~ds+2k32mem∫0tEhs-ysds+2k32mem∫0tEysys^-ys2q-1ds+2k32Tm+k2∫0tEhs-xsds≤k21-θTC^1Δ+k2θTC~1Δ+2k32TC1memΔ+k2∫0tEhs-xsds+2k32mem∫0tEys2Eys^-ys4q-2ds+2k32Tm≤k2T1-θC^1+θC~1Δ+2k32TC1memΔ+k2∫0tEhs-xsds+2k32TmemA2Ex0+k1T2A^4q-2Δq-1/2+2k32Tm, where we have used Lemma 5 in the second step and the Hölder inequality and Lemmas 3 and 4 in the third step and the fact that EMt=0 (the function d(u)=ϕm′(hu-xu)(g(yu^,yu~)-g(xu,xu)) belongs to the space M2([0,t];R) of real-valued measurable Ft-adapted processes such that E∫0t‍|d(u)|2du<∞; thus [20, Theorem 1.5.8] implies EMt=0). Thus (45) becomes (51)Eht-xt≤em-1+J1Δ+2k32TC1Δmem+J3Δq-1/2mem+2k32T1m+k2∫0tEhs-xsds≤J3Δq-1/2mem+2k32T1m+k2∫0tEhs-xsds≤J3Δq-1/2mem+2k32T1mek2t, where in the second step we have used the asymptotic relations Δκ=o(Δq-1/2) as Δ↓0 for any κ≥1/2,em-1=o(1/m) as m→∞,Δ=o(Δκ/mem) for any κ≤1 as m→∞; in the last step we have used the Gronwall inequality and J3 is as defined in Proposition 6 while(52)J1≔k2T1-θC^1+θC~1. Taking the supremum over all 0≤t≤T gives (41).


### 3.3. Convergence of the Auxiliary Process(ht) to (xt) in L2

Proposition 7.
Let Assumption A hold. Then one has(53)Esup0≤t≤Tht-xt2≤Cϵln⁡Δ-1, where Cϵ is independent of Δ and given by Cϵ≔323/2ϵ(k3)4T2e6T2(k2)2+k2T, where 0<ϵ<q-1/2.Proof of Proposition7.
We estimate the difference|Et|2≔|ht-xt|2. It holds that (54)Et2=∫0tfθys^,ys~-fθxs,xsds+∫0tgys^,ys-gxs,xsdWs2≤2T∫0tk21-θhs-ys^+k2θhs-ys~+k2Es2ds+2Mt2≤6Tk221-θ2∫0ths-ys^2ds+6Tk22θ2∫0ths-ys~2ds+6Tk22∫0tEs2ds+2Mt2,where in the second step we have used the Cauchy-Schwarz inequality and (46) and (55)Mt≔∫0tgyu^,yu~-gxu,xudWu. Taking the supremum over all t∈[0,T] and then expectations we have (56)Esup0≤t≤T⁡Et2≤6Tk221-θ2∫0TEhs-ys^2ds+6Tk22θ2∫0TEhs-ys~2ds+6Tk22∫0TEsup0≤l≤s⁡El2ds+2Esup0≤t≤T⁡Mt2≤6T2k221-θ2A^2Δ+6T2k22θ2A~2Δ+6Tk22∫0TEsup0≤l≤s⁡El2ds+8EMT2, where in the second step we have used Lemma 4 and Doob’s martingale inequality with p=2, since Mt is an R-valued martingale that belongs to L2. We find that (57)EMT2≔E∫0Tgys^,ys-gxs,xsdWs2=E∫0Tgys^,ys-gxs,xs2ds≤2k32E∫0Tysys^-ys2q-1+hs-ys+hs-xsds≤2k32∫0TEysys^-ys2q-1ds+2k32∫0TEhs-ysds+2k32∫0TEEsds, where we have used (47). Now, Lemmas 3, 4, and 5 imply (58)EMT2≤J6Δ2q-1+2k32TC1Δ+2k32∫0TEEsds≤J6Δq-1/2+2k32∫0TEEsds, where we have used the asymptotic relations Δl=o(Δq-1/2) for all l≥1/2 as Δ↓0 and the quantity J6 is given by J6≔2(k3)2TA2E(x0+k1T)2A^4q-2.
Relation (56) becomes (59)Esup0≤t≤T⁡Et2≤8J6Δq-1/2+J5Δ+6Tk22∫0TEsup0≤l≤s⁡El2ds+16k32∫0TEsup0≤l≤s⁡Elds≤8J6Δq-1/2+16k32TJ3Δq-1/2mem+3k32T1mek2T+6Tk22∫0TEsup0≤l≤s⁡El2ds≤16k32TJ3ek2T+6T2k22Δq-1/2mem+32k34T2ek2T+6T2k221m, where we have used Proposition 6 in the second step with the sequence em as defined there and Gronwall’s inequality in the last step and the asymptotic relation Δκ=o(Δκ/mem) as m→∞, for any κ>0, and J5 is independent of Δ and given by J5≔6T2(k2)2((1-θ)2A^2+θ2A~2).
We takem=ln⁡Δ-λ, with λ>0, to be specified soon and note that eln⁡Δ-λ=o(Δ-λ) as Δ↓0, since eln⁡n=o(n) as n→∞. Moreover we have that (60)Δq-1/2em=Δq-1/2e-m2/2em/2=Δq-1/2e-ln⁡Δ-λ/2e1/2ln⁡Δ-λ=Δq-1/2-3λ/2e1/2ln⁡Δ-λΔ-λ. Now, since q>1/2 there is an ϵ>0 small enough such that q-1/2-ϵ>0. We take λ=2ϵ/3 and conclude that (61)Δq-1/2em=Δq-1/2-ϵe1/2ln⁡Δ2ϵ/3Δ-2ϵ/3⟶0, as Δ→0 which in turn implies the asymptotic relation Δq-1/2/mem=o(1/m) as m→∞, with the logarithmic rate stated before. In the same way we can show Δ1/4/mem=o(1/m) as m→∞, by taking ϵ<1/4. We finally arrive at (62)Esup0≤t≤T⁡Et2≤32k34T2ek2T+6T2k221ln⁡Δ-2ϵ/3, by taking 0<ϵ<q-1/2, which implies (53).


### 3.4. Proof of Theorem1

In order to finish the proof of Theorem1 we just use the triangle inequality, Lemma 5, and Proposition 7 to get(63)Esup0≤t≤T⁡yt-xt2≤2Esup0≤t≤T⁡ht-yt2+2Esup0≤t≤T⁡Et2≤2C2Δ2+2Cϵln⁡Δ-1≤Cln⁡Δ-1,where C=C(k2,k3,ϵ,T), given in the statement of Theorem 1.


## 3.1. Moment Bounds

Lemma 3 (moment bound for SD approximation).
It holds that(20)Esup0≤t≤Tytp≤ApEx0+k1Tp, for any p>2, where Ap≔expp(p-1)/2(k3)2p-1/2p+2p-1/pT.Proof of Lemma3.
We first observe that(yt) is bounded in the following way: (21)0≤yt≤x0+∫0tk1ds+k3∫0tys^q-1/2ysdWs≤x0+k1T+k3∫0tys^q-1/2ysdWs≔ut, a.s., where the lower bound comes from the construction of (yt) and the upper bound follows from a comparison theorem. We will bound (ut) and therefore (yt), since 0≤yt≤ut a.s. Set the stopping time τR≔inf⁡{t∈[0,T]:ut>R}, for R>0 with the convention inf⁡∅=∞. Application of Itô’s formula on (ut∧τR)p implies (22)ut∧τRp=x0+k1Tp+pp-12k32∫0t∧τRusp-2ys^2q-1ysds+pk3∫0t∧τRusp-1ys^q-1/2ysdWs≤x0+k1Tp+pp-12k32∫0t∧τRusp-1ys^2q-1ds+Mt≤x0+k1Tp+pp-12k32∫0t∧τRp-12pusp+2p-1pys^2q-1pds+Mt≤x0+k1Tp+pp-12k32p-12p+2p-1p∫0t∧τRuspds+Mt, where in the second step we have used the fact that 0≤yt≤ut, in the third step the inequality xp-1y≤ϵp-1/pxp+1/pϵp-1yp, valid for x∧y≥0 and p>1 with ϵ=1/2, and in the final step the fact that 1/2<q<1 and Mt≔pk3∫0t∧τR‍(us)p-1(ys^)q-1/2ysdWs. Taking expectations in the above inequality and using that Mt is a local martingale vanishing at 0, we get (23)Eut∧τRp≤Ex0+k1Tp+pp-12k32p-12p+2p-1p∫0tEus∧τRpds≤Ex0+k1Tpexp⁡pp-12k32p-12p+2p-1pT≤ApEx0+k1Tp, where we have applied the Gronwall inequality [18, (7)]. We have that (24)yt∧τRp=yτRpIτR≤t+ytpIt<τR≥ytpIt<τR. Thus taking expectations in the above inequality and using the estimated upper bound for E(ut∧τR)p we arrive at (25)EytpIt<τR≤Eyt∧τRp≤Eut∧τRp≤ApEx0+k1Tp,and, taking the limit as R→∞, we get (26)limR→∞⁡EytpIt<τR≤ApEx0+k1Tp. Let us fix t. The sequence of stopping times τR is increasing in R and t∧τR→t as R→∞, and thus the sequence (yt)pI{t<τR} is nondecreasing in R and (yt)pI{t<τR}→(yt)p as R→∞. Application of the monotone convergence theorem implies (27)Eytp≤ApEx0+k1Tp, for any p>2. Using again Itô’s formula on (ut)p, taking the supremum, and then using Doob’s martingale inequality on the diffusion term we bound Esup0≤t≤T(ut)p and thus Esup0≤t≤T(yt)p.Lemma 4 (error bound for SD scheme).
Letns be an integer such that s∈[tns,tns+1]. Then (28)Eys-ys^p≤A^pΔp/2,Eys-ys~p<A~pΔp/2, for any p>0, where the positive quantities A^p,A~p do not depend on Δ.Proof of Lemma4.
First we take ap>2. We get that(29)ys-ys^p=∫tnssk1-k21-θyu^-k2θyu~du+∫tnstns+1k2θys^du-∫stns+1k2θysdu+∫stnsk1-k21-θytns-k3241+k2θΔytns2q-1du+k3∫tnssyu^q-1/2yudWup≤5p-1∫tnssk1-k21-θyu^-k2θyu~dup+k2pθpys^ptns+1-tnsp+k2pθpysptns+1-sp+∫stnsk1-k21-θytns-k3241+k2θΔytns2q-1dup+k3p∫tnssyu^q-1/2yudWup≤5p-1tns-sp-1∫tnssk1-k21-θyu^-k2θyu~pdu+k2pθpys^p+yspΔp+k1-k21-θytns-k3241+k2θΔytns2q-1pΔp+k3p∫tnssyu^q-1/2yudWup, where we have used the Cauchy-Schwarz inequality. Taking expectations in the above inequality and using Lemma 3 and Doob’s martingale inequality on the diffusion term we conclude (30)Eys-ys^p≤A^pΔp/2, where the positive quantity A^p, except on p, depends also on the parameters k1,k2,k3,θ,q, but not on Δ. Now, for 0<p<2, we get (31)Eys-ys^p≤Eys-ys^2p/2≤A^pΔp/2, where we have used Jensen’s inequality for the concave function ϕ(x)=xp/2. Following the same lines, we can show that (32)Eys-ys~p≤A~pΔp/2, for any 0<p, where the positive quantity A~p, except on p, depends also on the parameters k1,k2,k3,θ,q, but not on Δ.For the rest of this section we rewrite again the compact form of (12) in the following way: (33)ytSD=x0+∫0tfθys^,ys~ds+∫0tgys^,ysdWs︸ht+∫ttn+1f1ytn,ytds,where fθ(·,·) is given by (3) and the auxiliary process (ht) is close to (yt) as shown in the next result.Lemma 5 (moment bounds involving the auxiliary process).
For anys∈[0,T] it holds that (34)Ehs-ysp≤CpΔp,Ehsp≤Ch,and for s∈[tn,tn+1] one has that (35)Ehs-ys^p≤C^pΔp/2,Ehs-ys~p≤C~pΔp/2,for any p>0, where the positive quantities Cp,C^p,C~p,Ch do not depend on Δ.Proof of Lemma5.
We have that(36)hs-ysp=∫stn+1f1ytn,ysdup≤tn+1-spf1ytn,ysp, for any p>0, where we have used (33). Using Lemma 3 we get the left part of (34). Now for p>2 and noting that (37)Ehsp≤2p-1Ehs-ysp+2p-1Eysp≤2p-1CpΔp+2p-1ApEx0+k1Tp≤Ch,we get the right part of (34), where we have used Lemma 3. The case 0<p<2 follows by Jensen’s inequality as in Lemma 4.
Furthermore, fors∈[tn,tn+1] and p>2, we derive that (38)Ehs-ys^p≤2p-1Ehs-ysp+2p-1Eys-ys^p≤2p-1CpΔp+2p-1A^pΔp/2≤C^pΔp/2, where we have used (30) and in the same manner (39)Ehs-ys~p≤2p-1CpΔp+2p-1A~pΔp/2≤C~pΔp/2. The case 0<p<2 follows by Jensen’s inequality.


## 3.2. Convergence of the Auxiliary Process(ht) to (xt) in L1

We will use the representation (33) and write (40)ht-xt=∫0tfθys^,ys~-fθxs,xsds+∫0tgys^,ys-gxs,xsdWs.Proposition 6.
Let Assumption A hold. Then one has(41)sup0≤t≤TEht-xt≤J3Δq-1/2mem+2k32T1mek2T, for any m>1, where em=e-m(m+1)/2 and (42)J3≔2k32TA2Ex0+k1T2A^4q-2.Proof of Proposition6.
Let the nonincreasing sequence{em}m∈N with em=e-m(m+1)/2 and e0=1. We introduce the following sequence of smooth approximations of |x| (method of Yamada and Watanabe [19]): (43)ϕmx=∫0xdy∫0yψmudu, where the existence of the continuous function ψm(u) with 0≤ψm(u)≤2/(mu) and support in (em,em-1) is justified by ∫emem-1‍(du/u)=m. The following relations hold for ϕm∈C2(R,R) with ϕm(0)=0,(44)x-em-1≤ϕmx≤x,ϕm′x≤1,x∈R,ϕm′′x≤2mx,whenem<x<em-1,ϕm′′x=0otherwise.We have that(45)Eht-xt≤em-1+Eϕmht-xt.Moreover we find that (46)fθys^,ys~-fθxs,xs=k1-k21-θys^-k2θys~-k1-k2xs=-k21-θys^-xs-k2θys~-xs=k21-θhs-ys^+k2θhs-ys~-k2hs-xs,(47)gys^,ys-gxs,xs2=k3ys^q-1/2ys-k3xsq2≤k32ysys^q-1/2-ysq-1/2+ysq-xsq2≤2k32ysys^q-1/2-ysq-1/22+ysq-xsq2≤2k32ysys^-ys2q-1+ys-xs2≤2k32ysys^-ys2q-1+hs-ys+hs-xs,where we have used properties of Hölder continuous functions and, namely, the fact that xq is q-Hölder continuous for q≤1, that is, |xq-yq|≤|x-y|q, and that xq is 1/2-Hölder continuous since q>1/2. Application of Itô’s formula to the sequence {ϕm}m∈N implies (48)ϕmht-xt=∫0tϕm′hs-xsfθys^,ys~-fθxs,xsds+Mt+12∫0tϕm′′hs-xsgys^,ys~-gxs,xs2ds≤∫0tk21-θhs-ys^+k2θhs-ys~+k2hs-xsds+Mt+2k32∫0t1mhs-xsysys^-ys2q-1+hs-ys+hs-xsds≤k21-θ∫0ths-ys^ds+k2θ∫0ths-ys~ds+2k32mem∫0ths-ysds+k2∫0ths-xsds+Mt+2k32mem∫0tysys^-ys2q-1ds+2k32Tm, where in the second step we have used (46) and (47) and the properties of ϕm and (49)Mt≔∫0tϕm′hu-xugyu^,yu~-gxu,xudWu.Taking expectations in the above inequality yields (50)Eϕmht-xt≤k21-θ∫0tEhs-ys^ds+k2θ∫0tEhs-ys~ds+2k32mem∫0tEhs-ysds+2k32mem∫0tEysys^-ys2q-1ds+2k32Tm+k2∫0tEhs-xsds≤k21-θTC^1Δ+k2θTC~1Δ+2k32TC1memΔ+k2∫0tEhs-xsds+2k32mem∫0tEys2Eys^-ys4q-2ds+2k32Tm≤k2T1-θC^1+θC~1Δ+2k32TC1memΔ+k2∫0tEhs-xsds+2k32TmemA2Ex0+k1T2A^4q-2Δq-1/2+2k32Tm, where we have used Lemma 5 in the second step and the Hölder inequality and Lemmas 3 and 4 in the third step and the fact that EMt=0 (the function d(u)=ϕm′(hu-xu)(g(yu^,yu~)-g(xu,xu)) belongs to the space M2([0,t];R) of real-valued measurable Ft-adapted processes such that E∫0t‍|d(u)|2du<∞; thus [20, Theorem 1.5.8] implies EMt=0). Thus (45) becomes (51)Eht-xt≤em-1+J1Δ+2k32TC1Δmem+J3Δq-1/2mem+2k32T1m+k2∫0tEhs-xsds≤J3Δq-1/2mem+2k32T1m+k2∫0tEhs-xsds≤J3Δq-1/2mem+2k32T1mek2t, where in the second step we have used the asymptotic relations Δκ=o(Δq-1/2) as Δ↓0 for any κ≥1/2,em-1=o(1/m) as m→∞,Δ=o(Δκ/mem) for any κ≤1 as m→∞; in the last step we have used the Gronwall inequality and J3 is as defined in Proposition 6 while(52)J1≔k2T1-θC^1+θC~1. Taking the supremum over all 0≤t≤T gives (41).


## 3.3. Convergence of the Auxiliary Process(ht) to (xt) in L2

Proposition 7.
Let Assumption A hold. Then one has(53)Esup0≤t≤Tht-xt2≤Cϵln⁡Δ-1, where Cϵ is independent of Δ and given by Cϵ≔323/2ϵ(k3)4T2e6T2(k2)2+k2T, where 0<ϵ<q-1/2.Proof of Proposition7.
We estimate the difference|Et|2≔|ht-xt|2. It holds that (54)Et2=∫0tfθys^,ys~-fθxs,xsds+∫0tgys^,ys-gxs,xsdWs2≤2T∫0tk21-θhs-ys^+k2θhs-ys~+k2Es2ds+2Mt2≤6Tk221-θ2∫0ths-ys^2ds+6Tk22θ2∫0ths-ys~2ds+6Tk22∫0tEs2ds+2Mt2,where in the second step we have used the Cauchy-Schwarz inequality and (46) and (55)Mt≔∫0tgyu^,yu~-gxu,xudWu. Taking the supremum over all t∈[0,T] and then expectations we have (56)Esup0≤t≤T⁡Et2≤6Tk221-θ2∫0TEhs-ys^2ds+6Tk22θ2∫0TEhs-ys~2ds+6Tk22∫0TEsup0≤l≤s⁡El2ds+2Esup0≤t≤T⁡Mt2≤6T2k221-θ2A^2Δ+6T2k22θ2A~2Δ+6Tk22∫0TEsup0≤l≤s⁡El2ds+8EMT2, where in the second step we have used Lemma 4 and Doob’s martingale inequality with p=2, since Mt is an R-valued martingale that belongs to L2. We find that (57)EMT2≔E∫0Tgys^,ys-gxs,xsdWs2=E∫0Tgys^,ys-gxs,xs2ds≤2k32E∫0Tysys^-ys2q-1+hs-ys+hs-xsds≤2k32∫0TEysys^-ys2q-1ds+2k32∫0TEhs-ysds+2k32∫0TEEsds, where we have used (47). Now, Lemmas 3, 4, and 5 imply (58)EMT2≤J6Δ2q-1+2k32TC1Δ+2k32∫0TEEsds≤J6Δq-1/2+2k32∫0TEEsds, where we have used the asymptotic relations Δl=o(Δq-1/2) for all l≥1/2 as Δ↓0 and the quantity J6 is given by J6≔2(k3)2TA2E(x0+k1T)2A^4q-2.
Relation (56) becomes (59)Esup0≤t≤T⁡Et2≤8J6Δq-1/2+J5Δ+6Tk22∫0TEsup0≤l≤s⁡El2ds+16k32∫0TEsup0≤l≤s⁡Elds≤8J6Δq-1/2+16k32TJ3Δq-1/2mem+3k32T1mek2T+6Tk22∫0TEsup0≤l≤s⁡El2ds≤16k32TJ3ek2T+6T2k22Δq-1/2mem+32k34T2ek2T+6T2k221m, where we have used Proposition 6 in the second step with the sequence em as defined there and Gronwall’s inequality in the last step and the asymptotic relation Δκ=o(Δκ/mem) as m→∞, for any κ>0, and J5 is independent of Δ and given by J5≔6T2(k2)2((1-θ)2A^2+θ2A~2).
We takem=ln⁡Δ-λ, with λ>0, to be specified soon and note that eln⁡Δ-λ=o(Δ-λ) as Δ↓0, since eln⁡n=o(n) as n→∞. Moreover we have that (60)Δq-1/2em=Δq-1/2e-m2/2em/2=Δq-1/2e-ln⁡Δ-λ/2e1/2ln⁡Δ-λ=Δq-1/2-3λ/2e1/2ln⁡Δ-λΔ-λ. Now, since q>1/2 there is an ϵ>0 small enough such that q-1/2-ϵ>0. We take λ=2ϵ/3 and conclude that (61)Δq-1/2em=Δq-1/2-ϵe1/2ln⁡Δ2ϵ/3Δ-2ϵ/3⟶0, as Δ→0 which in turn implies the asymptotic relation Δq-1/2/mem=o(1/m) as m→∞, with the logarithmic rate stated before. In the same way we can show Δ1/4/mem=o(1/m) as m→∞, by taking ϵ<1/4. We finally arrive at (62)Esup0≤t≤T⁡Et2≤32k34T2ek2T+6T2k221ln⁡Δ-2ϵ/3, by taking 0<ϵ<q-1/2, which implies (53).


## 3.4. Proof of Theorem1

In order to finish the proof of Theorem1 we just use the triangle inequality, Lemma 5, and Proposition 7 to get(63)Esup0≤t≤T⁡yt-xt2≤2Esup0≤t≤T⁡ht-yt2+2Esup0≤t≤T⁡Et2≤2C2Δ2+2Cϵln⁡Δ-1≤Cln⁡Δ-1,where C=C(k2,k3,ϵ,T), given in the statement of Theorem 1.


## 4. Polynomial Rate of Convergence

We work with the stochastic time change inspired by [21]. We define the process (64)γt≔∫0t128k32q2ys1-q+xs1-q2dsand the stopping time (65)τl≔inf⁡s∈0,T:6Tk22s+γs≥l. The process γ(t) is well defined since xt>0 a.s. and yt≥0 (see Section 2).The difference|Et|2≔|ht-xt|2 is estimated as in Section 3 and we get, as in (56), that (66)Esup0≤t≤τ⁡Et2≤J5Δ+6Tk22∫0τEsup0≤l≤s⁡El2ds+8EMτ2, where τ is a stopping time and J5 independent of Δ is as in the proof of Proposition 7. The main difference here will be the estimation of the last term in (66). The approach in Section 3 resulted in the L1 estimation E|Et| where we used the Yamada-Watanabe approach. Now, we use the Berkaoui approach. Relation (47) becomes (67)gys^,ys-gxs,xs2≤2k32ysys^-ys2q-1+ysq-xsq2≤2k32ysys^-ys2q-1+ysq-xsq2ys1-q+xs1-q2γs′64q2≤2k32ysys^-ys2q-1+18hs-ys2+Es2γs′, where we have used the inequality(68)aq-bqa1-q+b1-q≤2qa-b, valid for all a≥0,b≥0, and 1/2≤q≤1. Consequently, we get the upper bound (69)EMτ2≔E∫0τgys^,ys-gxs,xsdWs2≤J6Δq-1/2+18∫0τEhs-ys2γs′ds+18∫0τEEs2γs′ds≤J6Δq-1/2+18∫0τEhs-ys4Eγs′2ds+18∫0τEEs2γs′ds, where we used the Hölder inequality; J6 independent of Δ is as in proof of Proposition 7. Relation (66) becomes (70)Esup0≤t≤τ⁡Et2≤8J6Δq-1/2+6Tk22∫0τEsup0≤l≤s⁡El2ds+∫0τEhs-ys4Eγs′2ds+∫0τEEs2γs′ds≤8J6Δq-1/2+C4Δ2∫0τE128k32q2ys1-q+xs1-q22ds+∫0τEsup0≤l≤s⁡El26Tk22s+γs′ds≤8J6Δq-1/2+C4128k32q2Δ2∫0τE1xs2-2qds+∫0τEsup0≤l≤s⁡El26Tk22s+γs′ds, where we have used Lemma 5 in the second step. At this point we want to estimate the inverse moments of (xt) and to do so we consider the transformation v=x2-2q and apply Itô’s formula to get (71)vt=v0+∫0t1-2q1-qk32︸K0+21-qk1︸K1vs1-2q/2-2q-21-qk2︸K2vsds+∫0t2k31-q︸K3vsdWs, for t∈[0,T], where v0=(x0)2-2q>0. Denote the drift coefficient of the process (vt) by a(vt) and consider the function (72)αv≔av-λ+K2v+2q-1λ+K01/2q-1k12-2q/2q-1︸ηλv, where λ≥0. Some elementary calculations show that this function attains its minimum at v∗≔k1(2q-1)/η(λ)2-2q and α(v∗)=0, thus (73)av≥λ-K2+ηλv.Consider the process (ζt(λ)) defined through (74)ζtλ=ζ0+∫0tλ-K2+ηλζsds+∫0tK3ζsdWs, for t∈[0,T] with ζ0(λ)=v0. Process (74) is a square root diffusion process and when 2λ/(K3)2-1≥0 or (75)λ≥21-q2k32the process is a CIR process which remains positive if ζ0(λ)>0. By a comparison theorem [17, Proposition 5.2.18] we obtain that vt≥ζt(λ)>0 a.s. or (vt)-1≤(ζt(λ))-1 a.s. or equivalently (xt)2q-2≤(ζt(λ))-1 a.s. The inverse moment bounds of (ζt(λ)) follow [22, (3.1)]:(76)supt∈0,TEζtλp<∞,forp>-2λK32by choosing big enough λ and particularly such that (75) holds strictly. Therefore, (77)Esup0≤t≤τ⁡Et2≤8J6Δq-1/2+∫0τEsup0≤l≤s⁡El26Tk22s+γs′ds. Relation (77) for τ=τl implies (78)Esup0≤t≤τl⁡Et2≤8J6Δq-1/2+∫0τlEsup0≤l≤s⁡El26Tk22s+γs′ds≤8J6Δq-1/2+∫0lEsup0≤j≤u⁡Eτj2du≤8J6elΔq-1/2, where in the last step we have used Gronwall’s inequality. Using again relation (77) for τ=T and under the change of variables u=6T(k2)2s+γs we get (79)Esup0≤t≤T⁡Et2≤8J6Δq-1/2+∫06k22T2+γTEsup0≤j≤u⁡Eτj2du≤8J6Δq-1/2+∫0∞Esup0≤j≤u⁡I6k22T2+γT≥uEτj2du≤∫06k22T2Esup0≤j≤u⁡Eτj2du+∫6k22T2∞P6k22T2+γT≥uEsup0≤j≤u⁡Eτj2∣6k22T2+γT≥udu+8J6Δq-1/2≤8J6e6k22T2Δq-1/2+∫0∞PγT≥uEsup0≤j≤u⁡Eτj2du+8J6Δq-1/2≤16J6e6k22T2Δq-1/2+8J6Δq-1/2∫0∞PγT≥ueudu, where in the last steps we have used (78). We proceed by showing that u→P(γT≥u)eu∈L1(R+). Markov’s inequality implies (80)PγT≥u≤e-ϵuEeϵγT, for any ϵ>0. The following bound holds: (81)γT=∫0T128k32q2ys1-q+xs1-q2ds≤128k32q2∫0Txs2q-2ds,and thus(82)EeϵγT≤Eeϵ128k32q2∫0Txs2q-2ds, where -1<2q-2<0. It remains to bound the exponential inverse moments of (xt) defined through the stochastic integral equation (2). Exponential inverse moments for the CIR process are known [10, Theorem 3.1] and are given by (83)Eeδ∫0tζsλ-1ds≤CHKζ0-1/2νλ-νλ2+8δ/K32, for 0≤δ≤2λ/(K3)2-12(K3)2/8=:ν(λ)2(K3)2/8 where the positive constant CHK is explicitly given in [10, (10)] and depends on the parameters k2,k3,T,q, but is independent of ζ0. Thus the other condition that we require for parameter λ is (84)λ≥21-q2δk3+21-q2k32. When (84) is satisfied then (75) is satisfied too; thus there is actually no restriction on the coefficient δ in (83) since we can always choose appropriately a λ such that (84) holds. Relation (82) becomes(85)EeϵγT≤Eeϵ128k32q2∫0Tvs-1ds≤Eeϵ128k32q2∫0Tζsλ-1ds.We therefore require that (86)128k32q2ϵ≤νλ2K328 and can always find a ϵ>1, such that the above relation holds by choosing appropriately λ as discussed before.Relation (85) becomes (87)EeϵγT≤CHKζ0-νλ/2, and therefore (88)PγT≥u≤CHKx01-qνλe-ϵu, where λ is chosen such that (86) holds with ϵ>1. We conclude (89)Esup0≤t≤T⁡Et2≤16J6e6k22T2Δq-1/2+8J6CHKx01-qνλΔq-1/2∫0∞e1-ϵudu≤C·Δq-1/2,by choosing ϵ>1 where C=C(k1,k2,k3,T,q,ϵ)≔8J62e6(k2)2T2+CHK/ϵ-1(x0)(1-q)ν(λ), is as given in statement of Theorem 2.


## 5. Alternative Approach with Reduced Rate of Convergence

In this section we briefly discuss the case where instead of (12) we use directly (6). Then, Lemmas 3, 4, and 5 still hold; that is, the moment bounds and error bounds of (ytSD), as well as the moment bounds involving the auxiliary process (ht), are true. The proof of the convergence results follows the same lines as in Sections 3 and 4. The main difference is in the estimation (47) that now becomes (90)sgn⁡zsgys^,ys-gxs,xs2≤3k32ys^2q-1yssgn⁡zs-12+ysys^-ys2q-1+hs-ys+hs-xs. The first term on the right-hand side of the above inequality containing the sgn⁡(zs) will contribute in a negative way to the rate of convergence. We do not give all the details but just mention that in order to bound the expectation of that term, which can be done in the following way, (91)Eys^2q-1yssgn⁡zs-12=E4ytn2q-1ysIzs≤0≤4Eytn2q-1ys-ytn2q+4Eytn2qIzs≤0Iytn≤Δ1-2ξ+4Eytn2qIzs≤0Iytn>Δ1-2ξ≤4Eytn2q-1ys-ytn+4Δ2q-4qξ+4Eytn4qPzs≤0∩ytn>Δ1-2ξ, we need to estimate the probability of zt being negative when at the same time ytn>Δ1-2ξ, for 0<ξ<1/2.Lemma 8.
For everyt∈[tn,tn+1] it holds that(92)Pzt≤0∩ytn>Δ1-2ξ≤Ck2,k3,θ,ΔΔ, where Ck2,k3,θ,Δ≔k3/1-k2(2-θ)Δ and Δ(2-θ)<1/k2 and (k3)2/(1+k2θΔ)≤4k2. Relation (92) implies that P{zt≤0}∩{ytn>Δ1-2ξ}=O(Δ) as Δ↓0.Proof of Lemma8.
By the definition (7) of (zt) for t∈[tn,tn+1] and for 0<ξ<1/2, we have that (93)A≔zt≤0∩ytn>Δ1-2ξ=ytnq-1/2Wt-Wtn≤-21+k2θΔk3yn∩ytn>Δ1-2ξ⊆A1∪A2, where (94)A1≔Wt-Wtn≤-21+k2θΔk3ynytn-q+1/2∩ytn≥1,A2≔Wt-Wtn≤-21+k2θΔk3ynytn-q+1/2∩1>ytn>Δ1-2ξ. The following inclusion relations hold for the event A1:(95)A1⊆ΔWn≤-21+k2θΔk3ytn1-k2Δ1+k2θΔ-k32Δ41+k2θΔ2ytn2q-1ytn-q+1/2∩ytn≥1⊆ΔWn≤-21+k2θΔk31-k22-θΔ1+k2θΔ-k32Δ41+k2θΔ2⊆ΔWnt-tn≤-2k31-k22-θΔ1+k2θΔt-tn when Δ(2-θ)<1/k2 and (k3)2/(1+k2θΔ)≤4k2, where ΔWn≔Wt-Wtn. We obtain (96)PG≤-β=∫-∞-β12πe-u2/2du≤∫-∞-βe-u2/2du=∫β∞e-u2/2du≤1βe-β2/2, for every standard normal random variable G, where in the last step we have used [17, (9.20), page 112] valid for β>0. Using the fact that ΔWn/t-tn is a standard normal r.v. and ignoring the exponential term in (96), since its exponent is negative, we get that (97)PA1≤k321-k22-θΔt-tn≤Ck2,k3,θ,ΔΔ. The following inclusion relations hold for the event A2:(98)A2⊆ΔWn≤-21+k2θΔk3ytn1-k21-θΔ1+k2θΔ+k1Δ1+k2θΔ-k32Δ41+k2θΔ2ytn2q-1ytn-q+1/2∩1>ytn>Δ1-2ξ⊆ΔWn≤-21+k2θΔk3Δ1-2ξ1-k21-θΔ1+k2θΔ+k1-k3241+k2θΔΔ1+k2θΔ⊆ΔWnt-tn≤-2k31-k21-θΔ1+k2θΔt-tnΔ1/2-ξ when Δ(1-θ)<1/k2 and (k3)2/(1+k2θΔ)≤4k1. Using again (96) we have that (99)PA2≤k321-k21-θΔt-tnΔξ-1/2e-2/k321-k21-θΔ1+k2θΔ/t-tnΔ1-2ξ≤k321-k21-θΔΔξe-2/k321-k21-θΔ1+k2θΔΔ-2ξ. Taking probabilities in the inclusion relation (93) and using (97) and (99) we get (100)PA≤PA1+PA2≤Ck2,k3,θ,ΔΔ+k321-k21-θΔΔξe-2/k321-k21-θΔ1+k2θΔΔ-2ξ≤Ck2,k3,θ,ΔΔ, since Δξe-Δ-2ξ=o(Δ) as Δ↓0. Finally, note that Ck2,k3,θ,Δ=k3/1-k2(2-θ)Δ→k3 as Δ↓0 which justifies the O(Δ) notation (see, e.g., [24]).


## 6. Numerical Experiments

We discretize the interval[0,T] with a number of steps in power of 2. The semidiscrete (SD) scheme is given by (101)ytn+1SD=ytn1-k2Δ1+k2θΔ+k1Δ1+k2θΔ-k32Δ41+k2θΔ2ytn2q-1+k321+k2θΔytnq-1/2ΔWn2, for n=0,…,N-1, where ΔWn≔Wtn+1-Wtn are the increments of the Brownian motion which are Gaussian random variables with ΔWn~N(0,Δ).The ALF (Alfonsi) scheme [4, Section 3] is an implicit scheme which requires solving the nonlinear equation (102)Yn+1=ytn+1-qk1Yn+1-q/1-q-k2Yn+1-qk322Yn+1-1Δ+k31-qΔWn, and then computing ytn+1ALF=(Yn+1)1/1-q. The estimation of Yn+1 in (102) can be done, for example, with Newton’s method but requires a small enough Δ (in the CIR case, that is, when q=1/2, (102) simplifies to a solution of a quadratic equation). We also consider a scheme recently proposed in [16] using again the SD method, but in a different way, (103)ytn+1HALq=ytn1-k2Δ+k1Δ-qk32Δ2ytn2q-11-q+k31-qΔWn1/1-q,for n=0,…,N-1. Note the similarity in the expressions of (103) and the SD scheme (101) proposed here. This is not strange, because they both rely on the same way of splitting the drift coefficient. In particular, in the explicit HAL scheme, the following process is considered: (104)ytHALq=ytn+f~1ytn·Δ+∫tntf~2ysds+∫tntsgn⁡zsg~ysdWs, for t∈(tn,tn+1] with y0=x0 a.s. where now (105)fx=k1-k2x-qk32/2x2q-1︸f~1x+qk32/2x2q-1︸f~2x,g~x=k3xq,zt=ytn1-k2Δ+k1Δ-qk32Δ2ytn2q-11-q+k31-qWt-Wtn. A comparison with (3) and (4) shows that f~2(x)=2qf2(x) and g~(x)=g(x,x), for θ=0. We write (104) again as (106)ytHALq=ytn+k1-k2ytn-qk322ytn2q-1·Δ+∫tntqk322ys2q-1ds+k3∫tntsgn⁡zsysqdWs, and the process (106) is well defined when (107)k32≤2qk1,Δ≤22k2+qk32. The reader can compare again with (6) for θ=0. Solving for yt, we end up with ytHAL(q)=|zt|1/1-q. The main result in [16] is (108)EytHAL-xt2≤C·Δ2qq-1/2,when (107) holds, implying a rate of convergence at least q(q-1/2) which is bigger than the rate of convergence of the SD scheme proposed here which is at least 1/2(q-1/2) (see Theorem 2).We also consider two more linear-implicit schemes that were stated in the Introduction and discussed in AppendixB. Namely, we compare with the balanced implicit method (BIM) with appropriate weight functions to guarantee positivity ([25, Theorem 5.9]), which reads (109)ytn+1BIMq=ytn+k1Δ+k3ytnqΔWn+ΔWn1+k2Δ+k3ytnq-1ΔWn, and the balanced Milstein method (BMM) with the suggested weight functions [25, Theorem 5.9] that is given by (110)ytn+1BMMq=ytn+k1+Θ-1k2ytnΔ+k3ytnqΔWn+qk32/2ytn2q-1ΔWn21+Θk2Δ+qk32/2ytn2q-2Δ. We take the relaxation parameter Θ to be 1/2 as recommended in [25, (5.10)].We aim to show experimentally the order of convergence for the above positivity preserving methods for the estimation of the true solution of the CEV model (2), that is, the semidiscrete methods SD (101) and the HAL scheme (103), as well as the implicit ALF scheme (102) and the linear-implicit schemes BIM and BMM. The choice of the parameters is the same as in [15, Figure 6] with k3=0.4. In particular (x0,k1,k2,k3,q,T)=(1/16,1/16,1,0.4,3/4,1).Furthermore, we would also like to reveal the dependence of the order of the semidiscrete methods onq; that is, we want to verify our theoretical results and in particular the order shown in Theorem 2. We take the level of implicitness of the SD method (101) to be θ=1; that is, we consider the fully implicit scheme. We also discuss the fully explicit scheme, that is, when θ=0, but also an intermediate scheme θ=1/2, in Section 7.We want to estimate the endpointL2-norm ϵ=E|y(Δ)(T)-xT|2 of the difference between the numerical scheme evaluated at step size Δ and the exact solution of (2). For that purpose, we compute M batches of L simulation paths, where each batch is estimated by ϵ^j=1/L∑i=1L|yi,j(Δ)(T)-yi,j(ref)(T)|2 and the Monte Carlo estimator of the error is (111)ϵ^=1ML∑j=1M∑i=1Lyi,jΔT-yi,jrefT2, and requires M·L Monte Carlo sample paths. The reference solution is evaluated at step size 2-14 of the numerical scheme. For the SD case, we have shown in Theorems 1 and 2 that it strongly converges to the exact solution. We simulate 100·100=10000 paths, where the choice for L=100 is as in [28, page 118]. The choice of the number of trajectories M·L=104 is also considered in [26, Section 5] where a fundamental mean-square theorem is proved for SDEs with superlinear growing coefficients satisfying a one-side Lipschitz condition, but unfortunately it is not positivity preserving. Of course, the number of Monte Carlo paths has to be sufficiently large, so as not to significantly hinder the mean-square errors.We plot in alog2-log2 scale and error bars represent 98% confidence intervals. The results are shown in Table 1 and Figure 1. Table 1 does not present the computed Monte Carlo errors with 98% confidence, since they were at least 9 times smaller that the mean-square errors.Table 1
98% error and step size of fully implicit SD, HAL, ALF, BIM, and BMM scheme for (2) with (x0,k1,k2,k3,q,T)=(1/16,1/16,1,0.4,3/4,1).
StepΔ
SD(θ=1)
Rate
HAL
Rate
ALF
Rate
BIM
Rate
BMM
Rate
2
-
5
0.035161
—
0.03568
—
0.044343
—
0.033275
—
0.035819
−
2
-
7
0.035082
0.0016
0.035034
0.0131
0.044539
-
0.0032
0.033541
-
0.0057
0.035737
0.0017
2
-
9
0.034565
0.1071
0.035109
-
0.0015
0.020261
0.5682
0.032988
0.012
0.035192
0.0111
2
-
11
0.033205
0.029
0.033709
0.029
0.019532
0.0264
0.031724
0.0282
0.033731
0.0306
2
-
13
0.025078
0.2025
0.024978
0.2162
0.014600
0.2099
0.024998
0.1719
0.025292
0.2077Figure 1
Convergence of fully implicit SD, HAL, ALF, BIM, and BMM schemes applied to SDE (2) with parameters (x0,k1,k2,k3,q,T)=(1/16,1/16,1,0.4,3/4,1).In Table2 we present the computational times (we simulate with 3.06GHz Intel Pentium, 1.49GB of RAM in MATLAB R2014b Software. The random number generator is Mersenne Twister. The evaluated times do not include the random number generation time, since all the methods we compare involve the same amount of random numbers) of fully implicit SD, HAL, ALF, BIM, and BMM, for the same problem. Figure 2 shows the relation between the error and computer time consumption. As one can see from Table 2 the CPU times for ALF are at least 1000 times bigger than the other schemes; thus we choose in Figure 2 to restrict our attention to the rest of the methods.Table 2
Average computational time (in seconds) for a path, for different discretizations, for all considered positivity preserving methods for the mean-reverting CEV process (2) with q=3/4.
Step Δ
Time/path (in sec): Implicit SD
HAL
ALF
BIM
BMM
2
-
5
0.000013
0.0000164
0.0221883
0.0000174
0.0000196
2
-
7
0.0000422
0.0000558
0.0841705
0.0000584
0.0000657
2
-
9
0.0001586
0.0002137
0.2453943
0.0002207
0.0002482
2
-
11
0.0006243
0.0008437
0.9768619
0.0008703
0.0009795
2
-
13
0.0024975
0.0033977
3.9096332
0.0034785
0.0039143Figure 2
Strong convergence error of the mean-reverting CEV process (2) as a function of CPU time (in sec) using positivity preserving schemes SD, HAL, ALF, BIM, and BMM with (x0,k1,k2,k3,q,T)=(1/16,100,0.05,1/16,1,0.4,3/4,1).We show, in Table3, the L2-distance between our proposed method and the other methods for the numerical approximation of (2). We work as before and estimate the distance, (112)dG,H=1ML∑j=1M∑i=1Lyi,jΔ,GT-yi,jΔ,HT2, between method G and H, by considering sufficient small Δ, and in particular for Δ=10-2,10-3,10-4.Table 3
TheL2-distance between all the considered numerical schemes applied to SDE (2) with parameter set (x0,k1,k2,k3,q,T)=(1/16,1/16,1,0.4,3/4,1).
Step Δ
98
%-d (SD, HAL)
98
%-d (SD, ALF)
98
%-d (SD, BIM)
98
%-d (SD, BMM)
1
0
-
2
0.0005727
0.0716140
0.0038373
0.0005312
1
0
-
3
0.0001577
0.0286630
0.0013460
0.0001564
1
0
-
4
0.0000498
0.0283117
0.0004448
0.0000498Finally, we examine the behavior of all the methods for a value of the parameterq close to 1/2. The results are shown in Table 4.Table 4
98% error and step size of fully implicit SD, HAL, ALF, BIM, and BMM scheme for (2) with (x0,k1,k2,k3,q,T)=(1/16,1/16,1,0.4,0.55,1).
Step Δ
SD-error(θ=1)
Rate
HAL-error
Rate
BIM-error
Rate
BMM-error
Rate
2
-
10
0.059129
—
0.059666
—
0.055819
—
0.059661
—
2
-
11
0.058178
0.0234
0.058669
0.0243
0.055033
0.0205
0.058667
0.0242
2
-
12
0.054582
0.092
0.054989
0.0935
0.051796
0.0875
0.054987
0.0935
2
-
13
0.044065
0.3088
0.044288
0.3122
0.042444
0.2873
0.044287
0.3122The following points of discussion are worth mentioning:(i)
The performance of all methods, as shown in Table1 and Figure 1, implies, in terms of error estimates, that the implicit ALF scheme performs better, for values of discretization steps Δ≤2-9. Actually for these step sizes the ALF method starts to converge, and the same is true for the HAL and BIM methods. All the methods except ALF, that is, the semidiscrete SD and HAL, the BIM, and the BMM have a similar behavior for all values of Δ in the sense of error estimation as Figure 1 shows. The similarity of SD, HAL, BIM, BMM, and especially between SD and BMM is also indicated in Table 3, where we see how close they are with respect to the L2-norm. Nevertheless, Table 3 also shows that in order to get an accuracy to at least two decimal digits, which in practice may be adequate concerning that we want, for example, to evaluate an option and thus our results are in euros, we are free to use any of the above available methods. We may then choose the fastest one, as will be discussed later on.(ii)
The experimental strong order of convergence of implicit SD for problem (2) is 1/5 (at least 1/2(q-1/2)=1/8 as shown theoretically and presented in Table 1). We also see that all methods converge with similar orders and the theoretically rate 1 of the ALF method [4] does not hold for these values of Δ. Thus, again we see that the rate in practical situations does not necessarily matter, if one has to consider very small values of Δ to achieve it. Moreover, we present in Table 5 the performance of the explicit SD method and see that it is very close to the implicit, which is of course natural to happen.(iii)
Table4 concerns the case where the parameter q is 0.55. We do not present the ALF method since it required smaller values of Δ. All the methods again behave quite the same, with the BIM performing better with respect to error estimation.(iv)
In practice, the computer time consumed to provide a desired level of accuracy is of great importance. In particular, in financial applications, a scheme is considered better when, except of its accuracy, it is implemented faster. As mentioned before, the SD method as well as the HAL method performs well in that aspect, compared to the implicit ALF method, which requires the estimation of a root of a nonlinear equation in each step and is therefore time consuming. This is presented in Table2 and Figure 2 which illustrates the advantage of the semidiscrete method SD, performing slightly better than HAL and BMM, better than BIM, and of course a lot better compared with ALF (over 1000 times quicker to achieve an accuracy of almost two decimal digits). Moreover, the explicit SD performs slightly better in that aspect.(v)
A negative step of a numerical method appears when the computer-generated random variable exceeds a certain threshold, which tends to increase as the step sizeΔ decreases. Thus, the undesirable effect of negative values that are produced by some numerical schemes (such as the explicit Euler (EM) and standard Milstein (M)) tends to disappear, since, after a certain small step size, the threshold exceeds the maximum standard normal random number attainable by the computer system.Table 5
The performance of fully explicit SD scheme (101) applied to SDE (2) with parameter set (x0,k1,k2,k3,q,T)=(1/16,1/16,1,0.4,3/4,1).
Step Δ
98
% SD-error (θ=0)
Rate
2
-
7
0.0344244
—
2
-
9
0.0342415
0.0038
2
-
11
0.0331273
0.0239
2
-
13
0.0250195
0.2025


## 7. Approximation of Stochastic Model (1)

So far we have focused on the process(Vt), which is one part of system (1). Nevertheless, it can be treated independently, since the only way that it interacts with the process (St) is through the correlation ρ of the Wiener processes. First we apply Itô’s formula on ln⁡(St) to get (113)ln⁡St=ln⁡S0+∫0tμdu-12∫0tVu2pdu+∫0tVupdWu,t∈0,T.Then, we consider two different schemes for the integration of (113) (the reason for not considering other schemes such as the two-dimensional Milstein is that they generally are time consuming, since they involve additional random number generation for the approximation of double Wiener integrals). The first is the EM scheme which reads (114)ln⁡Stn+1EM=ln⁡Stn+μΔ-12Vtn2pΔ+VtnpΔWn,which has strong convergence order 1/2 and is easy to implement. The second scheme, which is based on an interpolation of the drift term and an interpolation of the diffusion term, considering decorrelation of the diffusion term, including a higher order Milstein term [15, Section 4.2], is denoted IJK and is given by [15, (137)](115)ln⁡Stn+1IJK=ln⁡Stn+μΔ-14Vtn2p+Vtn+12pΔ+ρVtnpΔW~n+12Vtnp+Vtn+1pΔWn-ρΔW~n+12ρpk3Vtnq+p-1ΔW~n2-Δ.We therefore consider the EM scheme (114) combined with SD (101) and the IJK scheme (115) combined with SD (101) and compare with the case where the stochastic variance (p=1/2) is integrated with BMM scheme (110), for three different correlation parameters, ρ=0,ρ=-0.4, and ρ=-0.8 with S0=100,μ=0.05, as in [15, Section 5]. We present in Tables 6, 7, and 8 the errors, in the sense of distance (112), for all the above considered ways of numerical integration of process (St), for different step sizes, as well as the average computational time (in seconds) consumed for each discretization. We also give an illustrative representation just for one case (ρ=-0.8) in Figure 3.Table 6
98% error, step size, and average computational time of numerical integration of process(St) using log-Euler or IJK method with SD or BMM scheme for (1) with (x0,S0,μ,k1,k2,k3,q,T)=(1/16,100,0.05,1/16,1,0.4,3/4,1) and correlation ρ=0.
Step Δ
EM and SD-error(θ=0.5)
IJK and SD-error(θ=0.5)
EM and BMM-error(Θ=0.5)
IJK and BMM-error(Θ=0.5)
2
-
5
26.901
(
0.0000261
)
26.901
(
0.0000159
)
26.891
(
0.00002
)
26.890
(
0.0000294
)
2
-
7
27.288
(
0.0000919
)
27.288
(
0.0000492
)
27.277
(
0.0000676
)
27.277
(
0.0001043
)
2
-
9
27.298
(
0.0003595
)
27.297
(
0.0001843
)
27.289
(
0.0002610
)
27.288
(
0.0004081
)
2
-
11
25.057
(
0.0014255
)
25.058
(
0.0007309
)
25.051
(
0.0010309
)
25.051
(
0.0016191
)
2
-
13
19.441
(
0.0057322
)
19.441
(
0.0028928
)
19.442
(
0.0041177
)
19.442
(
0.0064721
)Table 7
98% error, step size, and average computational time of numerical integration of process(St) using log-Euler or IJK method with SD or BMM scheme for (1) with (x0,S0,μ,k1,k2,k3,q,T)=(1/16,100,0.05,1/16,1,0.4,3/4,1), correlation ρ=-0.4.
Step Δ
EM and SD-error(θ=0.5)
IJK and SD-error(θ=0.5)
EM and BMM-error(Θ=0.5)
IJK and BMM-error(Θ=0.5)
2
-
5
26.382
(
0.0000266
)
26.331
(
0.0000161
)
26.372
(
0.0000202
)
26.324
(
0.00003
)
2
-
7
26.448
(
0.0000951
)
26.396
(
0.000005
)
26.439
(
0.0000691
)
26.389
(
0.0001081
)
2
-
9
25.951
(
0.0003631
)
25.909
(
0.000184
)
25.944
(
0.0002606
)
25.904
(
0.0004131
)
2
-
11
24.540
(
0.0014506
)
24.494
(
0.0007355
)
24.531
(
0.0010378
)
24.486
(
0.0016495
)
2
-
13
18.738
(
0.0060748
)
18.749
(
0.0030185
)
18.735
(
0.0042868
)
18.747
(
0.0068395
)Table 8
98% error, step size, and average computational time of numerical integration of process(St) using log-Euler or IJK method with SD or BMM scheme for (1) with (x0,S0,μ,k1,k2,k3,q,T)=(1/16,100,0.05,1/16,1,0.4,3/4,1), correlation ρ=-0.8.
Step Δ
EM and SD-error(θ=0.5)
IJK and SD-error(θ=0.5)
EM and BMM-error(Θ=0.5)
IJK and BMM-error(Θ=0.5)
2
-
5
25.552
(
0.0000263
)
25.455
(
0.0000159
)
25.541
(
0.0000199
)
25.449
(
0.0000296
)
2
-
7
25.670
(
0.0000932
)
25.569
(
0.0000494
)
25.659
(
0.0000678
)
25.564
(
0.0001059
)
2
-
9
25.217
(
0.0003622
)
25.137
(
0.0001835
)
25.208
(
0.0002595
)
25.132
(
0.0004111
)
2
-
11
23.743
(
0.0014407
)
23.711
(
0.0007306
)
23.734
(
0.0010307
)
23.707
(
0.0016376
)
2
-
13
18.082
(
0.005871
)
18.316
(
0.0029312
)
18.078
(
0.0041637
)
18.312
(
0.0066239
)Figure 3
Strong convergence error of the financial underlying process(St), as a function of CPU time (in sec) using log-Euler or IJK method with SD or BMM scheme for (1) with (x0,S0,μ,k1,k2,k3,q,T)=(1/16,100,0.05,1/16,1,0.4,3/4,1), correlation ρ=-0.8.Tables6, 7, and 8 indicate that in all cases the favorable choice is to integrate (St) using IJK method combined with the SD scheme for (Vt) in model (1). The IJK-SD approximation of system (1) seems to be the better one, with respect to CPU time, for every correlation coefficient considered.


## 8. Conclusion

In this paper, we exploit further the semidiscrete method (SD), which originally appeared in [8], to numerically approximate stochastic processes that appear in financial mathematics and are meant to be nonnegative. In [23] we examined the Heston 3/2-model, that is, a mean-reverting process with superlinear diffusion, described by a SDE of the form (2) with q=3/2. Now, we deal with SDEs with sublinear diffusion coefficients of the type (xt)q with 1/2<q<1. These kinds of SDEs, called mean-reverting CEV processes, appear in stochastic models, where they represent the instantaneous volatility-variance of an underlying financially observable variable. We prove theoretically the strong convergence of our proposed SD scheme, revealing the order of convergence. The resulting polynomial rate is shown in Theorem 1. We present a comparative study between various positivity preserving schemes and the SD method seems to be the best with respect to CPU time consumption. The advantage of the SD method here is that although implicit, it has an explicit formula and thus requires fewer arithmetic operations and consequently less computational time. Moreover, our method can cover cases where (2) has time varying coefficients, that is, k1(t),k2(t),k3(t).We also treat the two-dimensional stochastic volatility model (1). In order to do that, we actually integrate the process ln(St) which satisfies a SDE of the form (113) and in the end transform back for (St). We only consider two different schemes for the integration of ln(St), namely, the Euler Maruyama (EM) scheme, which is easy to implement and the IJK scheme [15, (137)] which is shown to be the most efficient method, robust and simple as EM [15]. We do not apply other two-dimensional schemes, such as, for example, the Milstein scheme, since they are in general time consuming, as they involve approximations of double Wiener integrals which require additional random number generation. We therefore combine the EM scheme with SD ((114) and (101)) and the IJK scheme with SD ((115) and (101)) and compare with the case where the stochastic variance (p=1/2) is integrated with BMM scheme (110), for three different correlation parameters, ρ=0,ρ=-0.4, and ρ=-0.8 with S0=100,μ=0.05, as in [15, Section 5]. The combination IJK with SD seems to be the most favorable with respect to CPU time, for all the cases.



---

*Source: 513137-2015-11-23.xml*