Skip to content
BY 4.0 license Open Access Published by De Gruyter February 15, 2022

A note on the asymptotic stability of the semi-discrete method for stochastic differential equations

Nikolaos Halidias ORCID logo and Ioannis S. Stamatiou ORCID logo

Abstract

We study the asymptotic stability of the semi-discrete (SD) numerical method for the approximation of stochastic differential equations. Recently, we examined the order of 2 -convergence of the truncated SD method and showed that it can be arbitrarily close to 1 2 ; see [I. S. Stamatiou and N. Halidias, Convergence rates of the semi-discrete method for stochastic differential equations, Theory Stoch. Process. 24 2019, 2, 89–100]. We show that the truncated SD method is able to preserve the asymptotic stability of the underlying SDE. Motivated by a numerical example, we also propose a different SD scheme, using the Lamperti transformation to the original SDE. Numerical simulations support our theoretical findings.

1 Introduction

We study the following class of scalar stochastic differential equations (SDE):

(1.1) d x t = a ( t , x t ) d t + b ( t , x t ) d W t , t [ 0 , T ] ,

where a , b : [ 0 , T ] × are measurable functions such that (1.1) has a unique solution and x 0 is independent of all { W t } t 0 . We assume that the coefficients a ( t , x ) , b ( t , x ) in (1.1) depend explicitly on t. SDEs of the type (1.1) rarely have explicit solutions, and therefore the need for numerical approximations for simulations of the solution process x t is apparent. In the case of non-linear drift and diffusion coefficients, classical methods may fail to strongly approximate (in the mean-square sense) the solution of (1.1); cf. [8], where the Euler method may explode in finite time.

The semi-discrete (SD) method, originally proposed in [1], is a strongly converging numerical method with the qualitative property of domain preservation; if for instance the solution process x t is nonnegative, then the approximation process y t is also nonnegative (see [6, 2, 4, 3, 5, 12, 13]). The 2 -convergence of the truncated SD method, see [14], was recently shown to be arbitrarily close to 1 2 .

The main idea behind the semi-discrete method is freezing on each subinterval appropriate parts of the drift and diffusion coefficients of the solution at the boundaries of the subinterval, ending up with an explicitly solvable SDE.

Our main goal is to further examine qualitative properties of the SD method relevant to the stability of the method and answer questions of the following type: Is the SD method able to preserve the asymptotic stability of the underlying SDE?

The answer of the question above is in the positive, and is given in our main result, Theorem 3.3. In Section 2, we give all necessary information about the truncated version of the semi-discrete method; the way of construction of the numerical scheme and some useful properties. Section 3 contains the main result with the proof. Section 4 provides a numerical example. Motivated by the SDE appearing in the example, we also propose a different SD scheme, using the Lamperti transformation to the original SDE. Numerical simulations support our theoretical findings. Finally, Section 5 contains concluding remarks.

2 Setting and assumptions

We will use some standard notation and setting that has already appeared in previous works of the SD method in a condensed way; we refer the interested reader mainly to [14]. Let T > 0 and ( Ω , , { t } 0 t T , ) be a complete probability space. Let also W t , ω : [ 0 , T ] × Ω be a one-dimensional Wiener process adapted to the filtration { t } 0 t T . We rewrite SDE (1.1) in its integral form

(2.1) x t = x 0 + 0 t a ( s , x s ) 𝑑 s + 0 t b ( s , x s ) 𝑑 W s , t [ 0 , T ] ,

and assume that it admits a unique strong solution.

We recall the SD scheme. We use the auxiliary real functions f ( s , r , x , y ) and g ( s , r , x , y ) with s , r [ 0 , T ] and x , y with the property f ( s , s , x , x ) = a ( s , x ) and g ( s , s , x , x ) = b ( s , x ) . Suppose the equidistant partition 0 = t 0 < t 1 < < t N = T with step-size Δ = T N and assume that for every n N - 1 the SDE

(2.2) y t = y t n + t n t f ( t n , s , y t n , y s ) 𝑑 s + t n t g ( t n , s , y t n , y s ) 𝑑 W s , t ( t n , t n + 1 ] ,

with y 0 = x 0 a.s., has a unique strong solution. The first and third variable in f , g denote the discretized part of the original SDE. Now, we proceed to the construction of the truncated SD method. We define the truncated functions f Δ and g Δ by

f Δ ( s , r , x , y ) := f ( s , r , ( | x | μ - 1 ( h ( Δ ) ) ) x | x | , y ) and g Δ ( s , r , x , y ) := g ( s , r , ( | x | μ - 1 ( h ( Δ ) ) ) x | x | , y )

for x , y , where we set x / | x | = 0 when x = 0 , the strictly increasing function μ : + + is such that

(2.3) sup | x | u ( | f ( s , r , x , y ) | | g ( s , r , x , y ) | ) μ ( u ) ( 1 + | y | ) , u 1 ,

for every s , r T , and the strictly decreasing function h : ( 0 , 1 ] [ μ ( 1 ) , ) satisfies

(2.4) lim Δ 0 h ( Δ ) = and Δ 1 / 6 h ( Δ ) h ^    for every  Δ ( 0 , 1 ] ,

for a constant h ^ 1 μ ( 1 ) . Let

(2.5) y t Δ = y t n Δ + t n t f Δ ( t n , s , y t n Δ , y s Δ ) 𝑑 s + t n t g Δ ( t n , s , y t n Δ , y s Δ ) 𝑑 W s , t ( t n , t n + 1 ] ,

with y 0 = x 0 a.s., admit a unique strong solution for every n N - 1 .

Assumption 2.1.

Let f Δ ( s , r , x , y ) , g Δ ( s , r , x , y ) satisfy the following condition:

| ϕ Δ ( s 1 , r 1 , x 1 , y 1 ) - ϕ Δ ( s 2 , r 2 , x 2 , y 2 ) | h ( Δ ) ( | s 1 - s 2 | + | r 1 - r 2 | + | x 1 - x 2 | + | y 1 - y 2 | )

for all 0 < Δ 1 and x 1 , x 2 , y 1 , y 2 , where h ( Δ ) is as in (2.4) and ϕ stands for f and g, respectively.

In the statement of the main result in [14], which we rewrite below, we use the compact form of the approximation process:

(2.6) y t Δ = y 0 + 0 t f Δ ( s ^ , s , y s ^ Δ , y s Δ ) 𝑑 s + 0 t g Δ ( s ^ , s , y s ^ Δ , y s Δ ) 𝑑 W s ,

where s ^ = t n when s [ t n , t n + 1 ) .

Theorem 2.2 (Order of strong convergence).

Let the coefficients of (2.1) satisfy the Khasminskii-type condition and suppose Assumption 2.1 holds and (2.5) has a unique strong solution for every n N - 1 , where x 0 L p ( Ω , R ) for some p 14 + 2 γ . Let ϵ ( 0 , 1 3 ) and define for γ > 0 ,

μ ( u ) = C ¯ u 1 + γ , u 0 , h ( Δ ) = C ¯ + ln Δ - ϵ , Δ ( 0 , 1 ] ,

where Δ 1 and h ^ are such that (2.4) holds. Then the semi-discrete numerical scheme (2.6) converges to the true solution of (2.1) in the L 2 -sense with order arbitrarily close to 1 2 , that is,

𝔼 sup 0 t T | y t Δ - x t | 2 C Δ 1 - ϵ .

3 Asymptotic stability

Now we are ready to study the ability of the truncated SD method to preserve the asymptotic stability of (2.1). For that reason, we also assume that a ( 0 , 0 ) = 0 and b ( 0 , 0 ) = 0 . Moreover, to guarantee the asymptotic stability of (2.1), we use an assumption similar to [7, Assumption 5.1].

Assumption 3.1.

We assume the existence of a continuous non-decreasing function κ : + + with κ ( 0 ) = 0 and κ ( u ) > 0 for all u > 0 such that

(3.1) 2 x a ( s , x ) + | b ( x ) | 2 - κ ( | x | )

for all x and s [ 0 , ) .

Now, we state a result without proof concerning the asymptotic stability of (2.1); see also [7, Theorem 5.2], where autonomous coefficients are assumed.

Theorem 3.2 (Asymptotic stability of underlying process).

Let Assumption 3.1 hold. Then the solution process of SDE (2.1) is asymptotically stable, that is,

lim t x t = 0 a.s.

for any x 0 R .

Recall equation (2.5), which defines the truncated SD numerical scheme. We rewrite our proposed scheme, that is, the solution of (2.5) at the discrete points 0 , t 1 , , t n + 1 , in the following way:

(3.2) y n + 1 Δ = Φ Δ ( y n Δ , t n , Δ , Δ W n ) ,

where Δ W n are the Wiener increments, Δ = t n + 1 - t n is the step-size and y n stands for y t n . We assume the following decomposition of Φ Δ for the above representation (3.2):

(3.3) Φ Δ ( y n Δ , t n , Δ , Δ W n ) 2 = ( y n Δ ) 2 + Φ 1 Δ ( y n Δ , t n , Δ ) + Φ 2 Δ ( y n Δ , t n , Δ , Δ W n ) ,

where

𝔼 ( Φ 2 Δ ( y n Δ , t n , Δ , Δ W n ) | t n ) = 0 .

The following theorem shows that the truncated SD method is able to preserve the asymptotic stability property of the underlying SDE.

Theorem 3.3 (Asymptotic numeric stability).

Let the auxiliary function Φ 1 Δ from (3.3) satisfy

(3.4) Φ 1 Δ ( y n Δ , t n , Δ ) - κ 1 ( | ( | y n Δ | μ - 1 ( h ( Δ ) ) ) y n Δ | y n Δ | | )

for any 0 < Δ Δ * , where κ 1 has the same properties as κ in (3.1) with κ 1 κ . Let also Assumption 3.1 hold. Then the solution of the truncated SD method (3.2) is numerically asymptotically stable, that is,

(3.5) lim n y n Δ = 0 a.s.

for all x 0 R and 0 < Δ Δ * .

Proof of Theorem 3.3.

Let us first fix a Δ ( 0 , Δ * ] . Set

π Δ ( x ) := ( | x | μ - 1 ( h ( Δ ) ) ) x | x | .

Then, combining (3.2), (3.3) and (3.4), we get

( y n + 1 Δ ) 2 ( y n Δ ) 2 - κ 1 ( | π Δ ( y n Δ ) | ) + Φ 2 Δ ( y n Δ , t n , Δ , Δ W n ) ( x 0 ) 2 - j = 0 n κ 1 ( | π Δ ( y j Δ ) | ) + M n ,

where

M n := j = 0 n Φ 2 Δ ( y j Δ , t j , Δ , Δ W j ) .

Recall that for the conditional expectation holds

𝔼 ( Φ 2 Δ ( y n Δ , t n , Δ , Δ W n ) | t n ) = 0

to find that M n , n = 0 , 1 , , is a martingale. Application of the nonnegative semi-martingale convergence theorem, cf. [10, Theorem 7, p. 139], implies

j = 0 κ 1 ( | π Δ ( y j Δ ) | ) < a.s. ,

which in turn yields

lim j κ 1 ( | π Δ ( y j Δ ) | ) = 0 a.s.

By the property of the function κ 1 , we get that

lim j ( | y j Δ | μ - 1 ( h ( Δ ) ) ) y j Δ | y j Δ | = 0 a.s.

Assertion (3.5) follows. ∎

Remark 3.4.

We would like to mention here that representation (3.2), which in turn defines decomposition (3.3) and consequently the function Φ 1 , depends heavily on the discretization (2.5), which enjoys the freedom of choice. Depending on the SDE at hand, we may discretize the coefficients in (2.5) in an additive or multiplicative way; cf. [1, 2, 3, 4, 5, 6, 12, 14]. We also present an example in the following section.

Remark 3.5.

One of the main advantages of the SD method is the domain preservation. Therefore, even if we may require no weaker conditions than the existing conditions to guarantee the stability of the SD method, we prefer to use it in cases where the solution process of the original SDE stays at a specific domain a.s.

4 Example

We will use the numerical example of [7, Example 5.4], that is, we consider an autonomous SDE of the form (2.1) with a ( x ) = - 10 x 3 and b ( x ) = x 2 , with initial condition x 0 , that is,

(4.1) x t = x 0 - 10 0 t x s 3 𝑑 s + 0 t x s 2 𝑑 W s , t 0 .

Using standard arguments, one may show that the solution process of SDE (4.1) is positive; see Section B. Assumption 3.1 holds with κ ( u ) = 19 u 4 . Therefore, by Theorem 3.2, SDE (4.1) is almost surely asymptotically stable. The classical Euler–Maruyama method is not able to reproduce this asymptotic stability; see [7, Appendix]. In the following, we show that the truncated SD method can reproduce this asymptotic stability. Since in the construction of the semi-discrete method the way of discretizing is not unique (but rather indicated by the equation itself), we will try two versions of the SD method by freezing different parts of the diffusion coefficient. We first choose the auxiliary functions f , g 1 and g 2 in the following way:

f ( s , r , x , y ) = - 10 x 2 y , g 1 ( s , r , x , y ) = x 2 , g 2 ( s , r , x , y ) = x y .

Thus, (2.2) becomes

(4.2) y t = y t n - 10 y t n 2 t n t y s 𝑑 s + y t n 2 t n t 𝑑 W s , t ( t n , t n + 1 ] ,

and

(4.3) y ^ t = y ^ t n - 10 y ^ t n 2 t n t y ^ s 𝑑 s + y ^ t n t n t y ^ s 𝑑 W s , t ( t n , t n + 1 ] ,

respectively, with y 0 = y ^ 0 = x 0 a.s. SDEs (4.2) and (4.3) are linear equations ((4.2) is linear in the narrow sense and is known as Langevin equation) with variable coefficients which admit a unique strong solution; cf. [9, Chapter 4.4] and Section A. In particular,

(4.4) y n + 1 = e - 10 y n 2 Δ ( y n + y n 2 t n t n + 1 e 10 y n 2 ( s - t n ) 𝑑 W s ) , n ,

and

(4.5) y ^ n + 1 = y ^ n exp { - 21 2 y ^ n 2 Δ + y ^ n Δ W n } , n .

Note that (2.3) holds with μ ( u ) = 10 | u | 2 since

sup | x | u ( | - 10 x 2 y | | x | | x 2 | ) 10 | u | 2 ( 1 + | y | ) , u 1 .

Therefore, in the notation of Theorem 2.2, γ = 1 and C ¯ = 10 . Finally,

h ( Δ ) = 10 + ln Δ - ϵ 1

for any Δ ( 0 , 1 ] . Clearly, h ( 1 ) μ ( 1 ) and

Δ 1 / 6 h ( Δ ) 10 Δ 1 / 6 + Δ 1 / 3 ln Δ - ϵ 1 11

for any Δ ( 0 , 1 ] and 0 < ϵ 1 1 3 . Therefore, we take h ^ = 11 . The truncated versions of the semi-discrete method (TSD) read

(4.6) y n + 1 Δ = e - 10 π Δ 2 ( y n Δ ) Δ ( y n Δ + π Δ 2 ( y n Δ ) t n t n + 1 e 10 π Δ 2 ( y n Δ ) ( s - t n ) 𝑑 W s )

and

(4.7) y ^ n + 1 Δ = y ^ n Δ exp { - 21 2 π Δ 2 ( y ^ n Δ ) Δ + π Δ ( y ^ n Δ ) Δ W n }

for n , where

π Δ ( x ) = ( | x | h ( Δ ) 10 ) x | x | ,

and therefore

π Δ 2 ( x ) = | x | 2 h ( Δ ) 10 .

4.1 Asymptotic stability of truncated semi-discrete method

Now, we compute ( y n + 1 Δ ) 2 taking the square of (4.6) and making some rearrangements to show that it admits representation (3.3):

( y n + 1 Δ ) 2 = ( y n Δ ) 2 e - 20 π Δ 2 ( y n Δ ) Δ + 2 y n Δ e - 10 π Δ 2 ( y n Δ ) t n + 1 π Δ 2 ( y n Δ ) t n t n + 1 e 10 π Δ 2 ( y n Δ ) s 𝑑 W s
+ π Δ 4 ( y n Δ ) e - 20 π Δ 2 ( y n Δ ) t n + 1 ( t n t n + 1 e 10 π Δ 2 ( y n Δ ) s 𝑑 W s ) 2
= ( y n Δ ) 2 - ( 1 - e - 20 π Δ 2 ( y n Δ ) t n + 1 ) ( y n Δ ) 2 + 1 20 π Δ 2 ( y n Δ ) e - 20 π Δ 2 ( y n Δ ) t n + 1 ( e 20 π Δ 2 ( y n Δ ) t n + 1 - e 20 π Δ 2 ( y n Δ ) t n )
+ 2 y n Δ e - 10 π Δ 2 ( y n Δ ) t n + 1 π Δ 2 ( y n Δ ) t n t n + 1 e 10 π Δ 2 ( y n Δ ) s 𝑑 W s + π Δ 4 ( y n Δ ) e - 20 π Δ 2 ( y n Δ ) t n + 1 ( t n t n + 1 e 10 π Δ 2 ( y n Δ ) s 𝑑 W s ) 2
- 1 20 π Δ 2 ( y n Δ ) e - 20 π Δ 2 ( y n Δ ) t n + 1 ( e 20 π Δ 2 ( y n Δ ) t n + 1 - e 20 π Δ 2 ( y n Δ ) t n ) .

Set

I := I ( y n Δ , t n , Δ , Δ W n ) = t n t n + 1 e 10 π Δ 2 ( y n Δ ) s 𝑑 W s

and set

Φ 2 Δ ( y n Δ , t n , Δ , Δ W n ) := 2 y n Δ e - 10 π Δ 2 ( y n Δ ) t n + 1 π Δ 2 ( y n Δ ) I + π Δ 4 ( y n Δ ) e - 20 π Δ 2 ( y n Δ ) t n + 1 I 2
- 1 20 π Δ 2 ( y n Δ ) e - 20 π Δ 2 ( y n Δ ) t n + 1 ( e 20 π Δ 2 ( y n Δ ) t n + 1 - e 20 π Δ 2 ( y n Δ ) t n )

to see that

𝔼 ( Φ 2 Δ ( y n Δ , t n , Δ , Δ W n ) | t n ) = 0 .

Moreover,

Φ 1 Δ ( y n Δ , t n , Δ ) := - ( 1 - e - 20 π Δ 2 ( y n Δ ) Δ ) ( y n Δ ) 2 + 1 20 π Δ 2 ( y n Δ ) ( 1 - e - 20 π Δ 2 ( y n Δ ) Δ )
- 19 20 ( 1 - e - 20 π Δ 2 ( y n Δ ) Δ ) π Δ 2 ( y n Δ ) ,

implying that we may choose κ 1 in the following way:

κ 1 ( u ) := 19 20 ( 1 - e - 20 u 2 Δ ) u 2 ,

so that condition (3.4) holds and therefore Theorem 3.3 applies. Note that κ 1 ( 0 ) = 0 and κ 1 ( u ) > 0 for any Δ > 0 . We conclude that the truncated SD scheme (4.6) preserves the asymptotic stability perfectly in the sense that lim n y n Δ = 0 a.s. for any 0 < Δ 1 .

4.2 Asymptotic stability of exponential truncated semi-discrete method

We examine ( y ^ n + 1 Δ ) 2 . We take the square of (4.7) and get that

( y ^ n + 1 Δ ) 2 = ( y ^ n Δ ) 2 e - 21 π Δ 2 ( y ^ n Δ ) Δ + 2 π Δ ( y ^ n Δ ) Δ W n
= ( y ^ n Δ ) 2 - ( 1 - e - 19 π Δ 2 ( y ^ n Δ ) Δ ) ( y ^ n Δ ) 2 + ( y ^ n Δ ) 2 e - 19 π Δ 2 ( y ^ n Δ ) Δ ( 1 - e - 2 π Δ 2 ( y ^ n Δ ) Δ + 2 π Δ ( y ^ n Δ ) Δ W n ) .

Set the last term of the above equality to Φ 2 Δ , that is,

Φ 2 Δ ( y ^ n Δ , t n , Δ , Δ W n ) := ( y ^ n Δ ) 2 e - 19 π Δ 2 ( y ^ n Δ ) Δ ( 1 - e - 2 π Δ 2 ( y ^ n Δ ) Δ + 2 π Δ ( y ^ n Δ ) Δ W n )

to see that

𝔼 ( Φ 2 Δ ( y n Δ , t n , Δ , Δ W n ) | t n ) = 0 ,

since

n := e - 2 π Δ 2 ( y ^ n Δ ) Δ + 2 π Δ ( y ^ n Δ ) Δ W n

is an exponential martingale. Moreover,

Φ 1 Δ ( y ^ n Δ , t n , Δ ) := - ( 1 - e - 19 π Δ 2 ( y ^ n Δ ) Δ ) ( y ^ n Δ ) 2 - ( 1 - e - 19 π Δ 2 ( y ^ n Δ ) Δ ) π Δ 2 ( y ^ n Δ ) ,

implying that we may choose κ 1 in the following way:

κ 1 ( u ) := ( 1 - e - 19 u 2 Δ ) u 2 ,

so that once more condition (3.4) holds and consequently Theorem 3.3 applies. We conclude that the truncated exponential SD scheme (4.7) preserves the asymptotic stability perfectly in the sense that lim n y ^ n Δ = 0 a.s. for any 0 < Δ 1 .

4.3 Semi-discrete method and Lamperti transformation

Instead of approximating directly (4.1), we first study a transformation of it, which produces a new SDE with constant diffusion coefficient; in other words, we use the Lamperti transformation of (4.1). In particular, consider z = - 1 x . Itô’s formula implies the following dynamics for ( z t ) (see Section C):

(4.8) z t = z 0 + 11 0 t ( z s ) - 1 𝑑 s + 0 t 𝑑 W s , t 0 .

Let t ( t n , t n + 1 ] and

(4.9) y ~ t = Δ W n + y ~ t n + 11 t n t ( y ~ s ) - 1 𝑑 s ,

with y ~ 0 = z 0 . Equation (4.9) is a Bernoulli-type equation with solution satisfying

(4.10) ( y ~ t ) 2 = ( Δ W n + y ~ t n ) 2 + 22 ( t - t n ) .

Recall that when x 0 > 0 , for the solution process holds x t > 0 a.s., which implies z t < 0 a.s., which in turn suggests that we take the negative root of (4.10) as the solution. Therefore, we propose the following semi-discrete method for the approximation of (4.8):

(4.11) y ~ t n + 1 = - ( Δ W n + y ~ t n ) 2 + 22 Δ ,

which suggests ( z ~ n ) n for the approximation of (4.1):

(4.12) z ~ t n + 1 = 1 ( Δ W n + y ~ t n ) 2 + 22 Δ .

We also examine another modification of the semi-discrete method, ( y ˘ t ) , where in each subinterval ( t n , t n + 1 ] we do not need to solve a new differential equation, but only an algebraic equation. For t ( t n , t n + 1 ] consider

(4.13) y ˘ t = W t - W t n + y ˘ t n + 11 ( y ˘ t ) - 1 Δ ,

with y ˘ 0 = z 0 . The solution of (4.13) satisfies

(4.14) ( y ˘ t ) 2 - ( W t - W t n + y ˘ t n ) y ˘ t - 11 Δ = 0 .

We propose the following version of the semi-discrete method for the approximation of (4.8):

(4.15) y ˘ t n + 1 = Δ W n + y ˘ t n - ( Δ W n + y ˘ t n ) 2 + 44 Δ 2 ,

which suggests ( z ~ n 2 ) n for the approximation of (4.1):

(4.16) z ~ t n + 1 2 = 2 - Δ W n - y ˘ t n + ( Δ W n + y ˘ t n ) 2 + 44 Δ .

4.4 Simulation paths

We present simulations for the numerical approximation of (4.1) with x 0 = 10 and compare with the truncated Euler–Maruyama method (TEM), which reads

(4.17) y n + 1 T E M = y n - 10 ( | y n | μ ¯ - 1 ( h ¯ ( Δ ) ) y n | y n | ) 3 Δ + ( | y n | μ ¯ - 1 ( h ¯ ( Δ ) ) ) 2 Δ W n

for n , where h ¯ ( Δ ) = Δ - 1 / 4 and μ ¯ ( u ) = 10 u 3 . According to the results in [7], it is shown that method (4.17) is asymptotically stable for any Δ 0.095 . Therefore, for such small step sizes we compare all methods presented here and for bigger Δ only the SD schemes (4.4), (4.5), (4.6) and (4.7). We also present the Lamperti semi-discrete schemes (LSD) (4.12) and (4.16). Moreover, the TEM method does not preserve positivity. Figures 1, 2 and 3 show sample simulations paths of TSD and TEM, respectively, for various stepsizes. Note that the truncated TSD, exponential truncated expTSD and the Lamperti semi-discrete schemes LSD1 (4.12) and LSD2 (4.16) work for all Δ < 1 .

Figure 1

Trajectories of (4.4)–(4.7), (4.12), (4.16) and (4.17) for the approximation of (4.1) with Δ < 0.095 . In (b) we find a zoom of Figure (a).

(a)

(a)

(b)

(b)

Figure 2

Trajectories of (4.4)–(4.7) and (4.12), (4.16) for the approximation of (4.1) with Δ = 0.25 . In (b) we find a zoom of Figure (a).

(a)

(a)

(b)

(b)

In the numerical simulation of the stochastic integral of the (truncated) TSD methods (4.4) and (4.6), we used the approximation

t n t n + 1 e 10 π Δ 2 ( y n Δ ) s 𝑑 W s e 10 π Δ 2 ( y n Δ ) t n Δ W n ,

that is, we calculated the integrand in the lower limit of integration. The above inequality is of the order Δ r , with 0 < r < 1 2 ; see Section D.

Figure 3

Trajectories of (4.4)–(4.7) and (4.12), (4.16) for the approximation of (4.1) with Δ = 0.5 . In (b) we find a zoom of Figure (a).

(a)

(a)

(b)

(b)

We also present in Figure 4 the difference between the Lamperti schemes and the truncated Euler–Maruyama scheme, y TEM - y LSD , for small enough Δ such that TEM works.

Figure 4

Difference of (4.17) with (4.12) and (4.16) for the approximation of (4.1) with various step sizes.

(a)

(a)

(b)

(b)

5 Conclusion and future work

In this paper, we study the asymptotic stability of the semi-discrete (SD) numerical method for the approximation of stochastic differential equations. Recently, we examined the order of 2 -convergence of the truncated SD method and showed that it can be arbitrarily close to 1 2 ; see [14]. We show that the truncated SD method is able to preserve the asymptotic stability of the underlying SDE. Motivated by a numerical example, we also propose a different SD scheme, where we actually approximate first the Lamperti transformation of the original SDE. This scheme preserves positivity (in this case) of the solution, has similar asymptotic properties as the other versions of the SD method and seems promising, since there is no need for an exponential calculation. We will study this numerical method and its properties in a forthcoming paper.

A Solution of linear SDEs in the narrow sense

Consider the following linear in the narrow sense SDE:

x t = x t 0 + t 0 t a x s 𝑑 s + t 0 t b 𝑑 W s

for t t 0 , where a , b are constants. Applying the Itô formula to the transformation U ( t , x ) = e - a ( t - t 0 ) x , we obtain

d U ( t , x ) = ( d e - a ( t - t 0 ) d t x t + a x t e - a ( t - t 0 ) ) d t + b e - a ( t - t 0 ) d W t
= b e - a ( t - t 0 ) d W t ,

or

e - a ( t - t 0 ) x t = x t 0 + b t 0 t e - a ( s - t 0 ) 𝑑 W s ,
x t = e a ( t - t 0 ) x t 0 + b e a ( t - t 0 ) t 0 t e - a ( s - t 0 ) 𝑑 W s .

B Positivity of (4.1)

In order to prove that x t > 0 a.s., we first show moment bounds of SDE (4.1).

Lemma B.1 (Uniform moment bounds for ( x t ) ).

The solution process ( x t ) of SDE (4.1) satisfies

𝔼 ( sup 0 t T ( x t ) p ) < A

for some A > 0 and any integer p with 2 p 19 2 .

Proof of Lemma B.1.

For all | x | > R with R > 1 , we have that

J ( x ) := x a ( x ) + ( p - 1 ) b 2 ( x ) / 2 1 + x 2
= x ( - 10 x 3 ) + ( p - 1 ) x 4 / 2 1 + x 2
= - 21 + p 2 x 4 1 + x 2 0 ,

where the last inequality is valid for all p such that p 21 . Thus J ( x ) is bounded for all x , since this clearly holds for all | x | R . Application of [11, Theorem 2.4.1] implies

𝔼 ( x t ) p 2 ( p - 2 ) / 2 ( 1 + ( x 0 ) p )

for any 2 p 21 , since x 0 . Using Itô’s formula on ( x t ) p , with p 19 2 (in order to use Doob’s martingale inequality later), we have that

( x t ) p = ( x 0 ) p + 0 t ( p ( x s ) p - 1 ( - 10 x s 3 ) + p ( p - 1 ) 2 ( x s ) p - 2 x s 4 ) 𝑑 s + 0 t p ( x s ) p - 1 x s 2 𝑑 W s
( x 0 ) p + p 0 t ( - 10 + p - 1 2 ) ( x s ) p + 2 𝑑 s + M t
( x 0 ) p + M t

for any even p with 2 p 21 , or p = 21 , where M t = 0 t p ( x s ) p + 1 𝑑 W s . Taking the supremum and then expectations in the above inequality, we get

𝔼 ( sup 0 t T ( x t ) p ) 𝔼 ( x 0 ) p + 𝔼 sup 0 t T M t
( x 0 ) p + 𝔼 sup 0 t T M t 2
( x 0 ) p + 4 𝔼 M T 2 ,

where in the last step we have used Doob’s martingale inequality to the diffusion term M t . ∎

Lemma B.2 (Positivity of ( x t ) ).

The solution process ( x t ) of SDE (4.1) is positive in the sense that x t > 0 a.s.

Proof of Lemma B.2.

Set the stopping time θ R = inf { t [ 0 , T ] : x t - 1 > R } for some R > 0 , with the convention that inf = . Application of Itô’s formula on ( x t θ R ) - 2 implies

( x t θ R ) - 2 = ( x 0 ) - 2 + 0 t θ R ( - 2 ) ( ( x s ) - 3 ( - 10 ) ( x s ) 3 + 3 ( x s ) - 4 ( x s ) 4 ) 𝑑 s + 0 t θ R ( - 2 ) ( x s ) - 3 x s 2 𝑑 W s
( x 0 ) - 2 + 0 t θ R 23 𝑑 s + 0 t ( - 2 ) x s - 1 𝑑 W s
( x 0 ) - 2 + 23 T + M t ,

where

M t := 0 t ( - 2 ) x s - 1 𝕀 ( 0 , t θ R ) ( s ) 𝑑 W s .

Taking expectations in the above inequality and using the fact that 𝔼 M t = 0 , we get that

𝔼 ( x t θ R ) - 2 ( x 0 ) - 2 + 23 T < C ,

with C independent of R. Therefore,

𝔼 ( 1 x t θ R 2 ) = R 2 ( θ R t ) + 𝔼 ( 1 x t 2 𝕀 ( t < θ R ) ) < C ,

implying that

( x t 0 ) = ( R = 1 { x t < 1 R } ) = lim R ( { x t < 1 R } ) lim R ( θ R t ) = 0 .

We conclude that x t > 0 a.s. ∎

C Lamperti transformation of (4.1)

Applying Itô’s formula to the transformation z ( x ) = - 1 x , we obtain

d z t = ( ( x t ) - 2 ( - 10 ) ( x t ) 3 + 1 2 ( - 2 ) ( x t ) - 3 ( x t ) 4 ) d t + ( x t ) - 2 ( x t ) 2 d W t
= - 11 x t + d W t
= 11 ( z t ) - 1 d t + d W t

or, for t t 0 ,

z t = z t 0 + 11 t 0 t ( z s ) - 1 + t 0 t 𝑑 W s
= z t 0 + 11 t 0 t ( z s ) - 1 + W t - W t 0 .

D Stochastic integral approximation

We want to estimate the stochastic integral appearing in the proposed truncated semi-discrete method (4.6) for the approximation of SDE (4.1). In a similar way, we calculate the integral appearing in the exponential truncated semi-discrete scheme (4.4).

In the numerical simulations, we used the following relation:

t n t n + 1 e 10 π Δ 2 ( y n Δ ) s 𝑑 W s e 10 π Δ 2 ( y n Δ ) t n Δ W n .

We show the following estimation:

(D.1) ( | t n t n + 1 e 10 π Δ 2 ( y n Δ ) s d W s - e 10 π Δ 2 ( y n Δ ) t n Δ W n | Δ r ) 2 e 20 π Δ 2 ( y n Δ ) t n + 1 Δ 1 - 2 r ,

suggesting that the probability of the absolute difference of these two random variables being of order Δ r , with 0 < r < 1 2 , approaches unity as Δ goes to zero. First, we write the difference of the two local martingales as

t n t n + 1 e 10 π Δ 2 ( y n Δ ) s 𝑑 W s - e 10 π Δ 2 ( y n Δ ) t n Δ W n = t n t n + 1 ( e 10 π Δ 2 ( y n Δ ) s - e 10 π Δ 2 ( y n Δ ) t n ) 𝑑 W s

and then use the martingale inequality to get for any ϵ > 0 that

( | t n t n + 1 ( e 10 π Δ 2 ( y n Δ ) s - e 10 π Δ 2 ( y n Δ ) t n ) d W s | ϵ ) ϵ - 2 𝔼 ( ( t n t n + 1 ( e 10 π Δ 2 ( y n Δ ) s - e 10 π Δ 2 ( y n Δ ) t n ) d W s ) 2 | t n )
ϵ - 2 t n t n + 1 𝔼 ( ( e 10 π Δ 2 ( y n Δ ) s - e 10 π Δ 2 ( y n Δ ) t n ) 2 | t n ) d s
ϵ - 2 [ 1 20 π Δ 2 ( y n Δ ) ( e 20 π Δ 2 ( y n Δ ) t n + 1 - e 20 π Δ 2 ( y n Δ ) t n ) + e 20 π Δ 2 ( y n Δ ) t n Δ ]
- 2 ϵ - 2 e 10 π Δ 2 ( y n Δ ) t n 1 10 π Δ 2 ( y n Δ ) ( e 10 π Δ 2 ( y n Δ ) t n + 1 - e 10 π Δ 2 ( y n Δ ) t n )
ϵ - 2 e 20 π Δ 2 ( y n Δ ) t n + 1 ( 1 20 π Δ 2 ( y n Δ ) ( 1 - e - 20 π Δ 2 ( y n Δ ) Δ ) + e - 20 π Δ 2 ( y n Δ ) Δ Δ )
2 ϵ - 2 e 20 π Δ 2 ( y n Δ ) t n + 1 Δ ,

where in the last step we used the inequality 1 - e - x x for any x > 0 . We apply the above inequality for ϵ = Δ r , with 0 < r < 1 2 to get (D.1).

References

[1] N. Halidias, Semi-discrete approximations for stochastic differential equations and applications, Int. J. Comput. Math. 89 (2012), no. 6, 780–794. 10.1080/00207160.2012.658380Search in Google Scholar

[2] N. Halidias, A novel approach to construct numerical methods for stochastic differential equations, Numer. Algorithms 66 (2014), no. 1, 79–87. 10.1007/s11075-013-9724-9Search in Google Scholar

[3] N. Halidias, Constructing positivity preserving numerical schemes for the two-factor CIR model, Monte Carlo Methods Appl. 21 (2015), no. 4, 313–323. 10.1515/mcma-2015-0109Search in Google Scholar

[4] N. Halidias, Construction of positivity preserving numerical schemes for some multidimensional stochastic differential equations, Discrete Contin. Dyn. Syst. Ser. B 20 (2015), no. 1, 153–160. 10.3934/dcdsb.2015.20.153Search in Google Scholar

[5] N. Halidias and I. S. Stamatiou, Approximating explicitly the mean-reverting CEV process, J. Probab. Stat. 2015 (2015), Article ID 513137, 20 pages. 10.1155/2015/513137Search in Google Scholar

[6] N. Halidias and I. S. Stamatiou, On the numerical solution of some non-linear stochastic differential equations using the semi-discrete method, Comput. Methods Appl. Math. 16 (2016), no. 1, 105–132. 10.1515/cmam-2015-0028Search in Google Scholar

[7] L. Hu, X. Li and X. Mao, Convergence rate and stability of the truncated Euler–Maruyama method for stochastic differential equations, J. Comput. Appl. Math. 337 (2018), 274–289. 10.1016/j.cam.2018.01.017Search in Google Scholar

[8] M. Hutzenthaler, A. Jentzen and P. E. Kloeden, Strong and weak divergence in finite time of Euler’s method for stochastic differential equations with non-globally Lipschitz continuous coefficients, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 467 (2011), no. 2130, 1563–1576. 10.1098/rspa.2010.0348Search in Google Scholar

[9] P. E. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations, Appl. Math. (New York) 23, Springer, Berlin, 1992. 10.1007/978-3-662-12616-5Search in Google Scholar

[10] R. S. Liptser and A. N. Shiryayev, Theory of Martingales, Math. Appl. 49, Springer, Cham, 1989. 10.1007/978-94-009-2438-3Search in Google Scholar

[11] X. Mao, Stochastic Differential Equations and Applications, 2nd ed., Horwood, Chichester, 2008. 10.1533/9780857099402Search in Google Scholar

[12] I. S. Stamatiou, A boundary preserving numerical scheme for the Wright–Fisher model, J. Comput. Appl. Math. 328 (2018), 132–150. 10.1016/j.cam.2017.07.011Search in Google Scholar

[13] I. S. Stamatiou, An explicit positivity preserving numerical scheme for CIR/CEV type delay models with jump, J. Comput. Appl. Math. 360 (2019), 78–98. 10.1016/j.cam.2019.04.005Search in Google Scholar

[14] I. S. Stamatiou and N. Halidias, Convergence rates of the semi-discrete method for stochastic differential equations, Theory Stoch. Process. 24 (2019), no. 2, 89–100. Search in Google Scholar

Received: 2021-08-01
Revised: 2022-01-23
Accepted: 2022-01-24
Published Online: 2022-02-15
Published in Print: 2022-03-01

© 2022 Walter de Gruyter GmbH, Berlin/Boston

This work is licensed under the Creative Commons Attribution 4.0 International License.