Stationary distribution of the stochastic theta method for nonlinear stochastic differential equations
Yanan Jiang Wei Liu Lihui Weng
math.NA math.PR









F
G
P


theoremTheorem[section]
lemma[theorem]Lemma
coro[theorem]Corollary
defn[theorem]Definition
assp[theorem]Assumption
cond[theorem]Condition
expl[theorem]Example
prop[theorem]Proposition
rmk[theorem]Remark
conj[theorem]Conjecture










α
ε ζ η øθ
ϑ κ̨ łλ μ ν
ξ  ρ̊ σ
ϕ φ   ω
√ ı ∀ ȷ∅

β δ υ ≡
∨ ∧


Δ ØΘ ŁΛ
Ξ Σ Ω
∂ ∇ ∃ ×
⋁ ⋀

⊘ ⊕ ⊗ ⊖
∘ ⊙ \ ∞
⋂ ⋃ ± ­∓
⟨ ⟩

Proof. 

trace
diag
  for 
=0.3in=1











Stationary distribution of the stochastic theta method for nonlinear stochastic differential equations
    Yanan Jiang, Wei Liu, Lihui Weng[Corresponding author, Email: 1151022187@qq.com]


Department of Mathematics, 
 Shanghai Normal University, Shanghai, China

    
================================================================================================================================================================




The existence and uniqueness of the stationary distribution of the numerical solution generated by the stochastic theta method is studied. When the parameter θ takes different values, the requirements on the drift and diffusion coefficients are different. The convergence of the numerical stationary distribution to the true counterpart is investigated. Several numerical experiments are presented to demonstrate the theoretical results.

Key words: stochastic theta method, nonlinear stochastic differential equations,

 numerical stationary distribution.





§ INTRODUCTION
 

The classical method to find the stationary distributions of some stochastic differential equations (SDEs) is to solve the corresponding Kolmogorov-Fokker-Planck equations. However, it is not trivial to find the solution to those partial differential equations when some nonlinearity appears in the drift or the diffusion coefficient of the SDEs. In this paper, the alternative path that the stationary distributions generated by some numerical methods for SDEs are used as the approximates to those of the underlying equations is investigated.
In the series papers <cit.>, the authors studied the approximates to stationary distributions of SDEs and SDEs with Markovian switching by using the Euler-Maruyama method. In <cit.>, the approximations of invariant measures of SDEs with different sorts of Markovian switchings were investigated using the Euler-Maruyama method. Both of the drift and diffusion coefficients of the SDEs in those papers above need to satisfy the global Lipschitz condition. As indicated in <cit.>, the classical Euler-Maruyama fails to convergence when either the drift or the diffusion coefficient grows super-linearly. To tackle this drawback, the backward Euler-Maruyama method was employed in <cit.> for those SDEs with the super-linear drift coefficient. Higher order methods were also discussed for sampling the invariant measures <cit.>.
In this paper, we study the numerical stationary distributions of the stochastic theta (ST) method and discuss the effect of the choice of the theta on the conditions of the coefficients. Different types of asymptotic properties of the stochastic theta method for SDEs have been widely investigated. In <cit.>, the stability of the ST method was studied in both the mean-square and almost sure senses. The stability of the ST method with nonrandom variable step sizes for bilinear, nonautonomous, homogenous test equations was investigated in <cit.>. The ST method was applied to a test system with stabilising and destabilising stochastic perturbations and almost sure asymptotic stability was analysed in <cit.>. The abilities to preserve the almost sure and the mean square exponential stabilities were discussed for different choices of the theta in <cit.> and <cit.>, respectively. The asymptotic boundedness of the ST method was studied in <cit.>. The results presented in this paper could be regarded as a complement to the existing study of the asymptotic behaviours of the ST method.
This paper is constructed in the following way. The necessary mathematical preliminaries are presented in Section <ref>. Section <ref> contains the main results. Some numerical examples are used to illustrate the theoretical results in Section <ref>. We summarize the paper by Section <ref>.



§ MATHEMATICAL PRELIMINARIES
 
In this paper, let (Ω , , ) be a complete probability space with a filtration {_t}_t ≥ 0 satisfying the usual conditions that it is right continuous and increasing while _0 contains all -null sets. Let |·| denote the Euclidean norm in ^d. The transpose of a vector or matrix, M, is denoted by M^T and the trace norm of a matrix, M, is denoted by |M| = √((M^T M)).
Let f,g : ^d →^d. To keep symbols simple, let B(t) be a scalar Brownian motion. The results in this paper can be extended to the case of multi-dimensional Brownian motions.
We consider the d-dimensional stochastic differential equation of the Itô type

    dx(t) = f(x(t))dt + g(x(t))dB(t)

with initial value x(0) = x_0.
Now, we present the conditions on the coefficients.



 Assume there exists a constant K_1>0 such that for any x,y ∈^d

    |f(x)-f(y)|^2∨|g(x)-g(y)|^2≤ K_1 |x-y|^2.




 Assume there exists a constant K_2<0 such that for any x,y ∈^d

    ⟨ x-y, f(x) - f(y)  ⟩≤ K_2 |x-y|^2.


In addition, we require that

    2K_2+K_1<0.

The next two conditions can be derived from Conditions <ref> and <ref> but with a little bit complicated coefficients. For the simplicity, we give two new conditions as follows.

 
There exist constants μ< 0 and a>0 such that for any x ∈^d

    ⟨ x,f(x)⟩≤μ|x|^2+a.





There exist positive constants σ, κ, b and c such that
for any x∈ℝ^d

    |g(x)|^2 ≤σ|x|^2+b,

and

    |f(x)|^2 ≤κ|x|^2+c,


In addition, we require that

    2μ+σ<0.

The existence and uniqueness of the underlying SDE (<ref>) has been broadly studied. We refer the readers to Theorem 3.1 in <cit.> for a quite general theory. There are other more general theories, the reason we refer the readers to this one is that the structure of it is similar to the following theory, Theorem <ref>.
The stochastic theta method to SDE (<ref>) is defined by

    X_k+1 = X_k + θ f(X_k+1)h +(1-θ) f(X_k)h+ g(X_k)  B_k,   X_0 = x(0)= x_0,

where B_k = B(t_k+1) - B(t_k) is the Brownian motion increment and t_k = kh, for k = 1,2,...
The proof of the next lemma is similar to those in <cit.>.


Let Condition <ref> <ref> hold and θ hK_2<1, the ST method (<ref>) is well defined.


It is useful to write (<ref>) as

    X_k+1-θ f(X_k+1)h  = X_k +(1-θ) f(X_k)h+ g(X_k)  B_k.

Define a function G:ℝ^d →ℝ^d by G(x)=x-f(x)θ h. Since

    ⟨ x-y,G(x)-G(y)⟩   ≥⟨ x-y,x-y-θ h(f(x)-f(y))⟩
       ≥|x-y|^2-θ hK_2|x-y|^2
       =(1-θ hK_2)|x-y|^2>0,

for  θ hK_2<1, we know that G has the
 inverse function G^-1:ℝ^d →ℝ^d. And G(x) is monotone. The ST method (<ref>) can be written as

    X_k+1 = G^-1(X_k +(1-θ) f(X_k)h+ g(X_k)  B_k).

Thus, the ST method (<ref>) is well defined. 


Let Conditions <ref> to <ref> hold, then

    ℙ(X_k+1∈ B|X_k=x)=ℙ(X_1∈ B|X_0=x),

for any Borel set B⊂ℝ^d.


 If X_k = x and X_0 = x, by (<ref>) we see

    X_k+1 - θ f(X_k+1)h = x + (1-θ)f(x)h+ g(x)  B_k,

and

    X_1 - θ f(X_1)h = x + (1-θ)f(x)h+ g(x)  B_0.

Because B_k and B_0 are identical in probability law, comparing the two equations above, we know that X_k+1 - θ f(X_k+1)h and X_1 - θ f(X_1)h
have the identical probability law. Then, due to Lemma <ref>, we have that X_k+1 and X_1 are identical in probability law under X_k = x and X_0 = x.
Therefore, the assertion holds.   


To prove Theorem <ref>, we cite the following classical result (see, for example, Lemma 9.2 on page 87 of <cit.>).


 Let h(x,ω) be a scalar bounded measurable random function of x, independent of _s. Let ζ be an _s-measurable random variable. Then

    (h(ζ,ω) |_s) = H(ζ),

where H(x) =  h(x,ω).

For any x ∈^d and any Borel set B ⊂^d, define

    (x,B) :=  (X_1 ∈ B | X_0 = x) and _k(x,B) :=  (X_k ∈ B | X_0 = x).



 
The solution generated by the ST method (<ref>) is a homogeneous Markov process with transition probability kernel (x,B).


 The homogeneous property follows Lemma <ref>, so we only need to show the Markov property. Define

    Y_k+1^x = G^-1 (x + (1 - θ)f(x)h + g(x)  B_k),

for x ∈^d and k ≥ 0. By (<ref>) we know that X_k+1 = Y_k+1^X_k. Let _t_k+1 = σ{ B(t_k+1) - B(t_k) }.
Clearly, _t_k+1 is independent of _t_k. Moreover, Y_k+1^x depends completely on the increment B(t_k+1) - B(t_k), so is _t_k+1-measurable.
Hence, Y_k+1^x is independent of _t_k. Applying Lemma <ref> with h(x,ω) = I_B(Y_k+1^x), we compute that

    (X_k+1∈ B |_t_k) = (I_B(X_k+1) |_t_k) = ( I_B(Y_k+1^X_k) |_t_k)
    = ( I_B(Y_k+1^x)) |_x=X_k
        = (x,B) |_x=X_k = (X_k,B) =  (X_k+1∈ B | X_k).

The proof is complete.   
Therefore, we see that (·,·) is the one-step transition probability and _k(·,·) is the k-step transition probability, both of which are induced by the BEM solution.
We state a simple version of the discrete-type Gronwall inequality in the next Lemma (see, for example, <cit.>).

Let {u_n} and {w_n} be nonnegative sequences, and α be a nonnegative constant. If

    u_n ≤α + ∑_k=0^n-1 u_k w_k   for n≥ 0,

then

    u_n ≤αexp(∑_k=0^n-1w_k).


Denote the family of all probability measures on ^d by (^d). Define by  the family of mappings F: ^d → satisfying

    |F(x) - F(y)| ≤ |x-y|    and    |F(x)| ≤ 1,

for any x,y ∈^d. For _1,_2 ∈(^d), define metric d_ by

    d_(_1,_2) = sup_F ∈| ∫_^d F(x)_1(dx) -  ∫_^d F(x)_2(dx) |.

The weak convergence of probability measures can be illustrated in terms of metric d_ <cit.>. That is,
a sequence of probability measures {_k}_k ≥ 1 in (^d) converge weakly to a probability measure ∈(^d) if and only if

    lim_k →∞ d_(_k,) = 0.

Then we define the stationary distribution for {X_k}_k ≥ 0 by using the concept of weak convergence.

 For any initial value x ∈^d and a given step size t > 0, {X_k}_k ≥ 0 is said to have a stationary distribution Π_ t∈(^d) if
the k-step transition probability measure _k(x,·) converges weakly to Π_ t(·) as k →∞ for every x ∈^d, that is

    lim_k →∞( sup_F ∈|  (F(X_k)) - E_Π_ t (F) | ) = 0,

where

    E_Π_ t (F) = ∫_^d F(y) Π_ t(dy).


In <cit.>, the authors presented a very general theory, Theorem 3.1, on the existence and uniqueness of the stationary distribution for any one step numerical methods. We adapt it here and state the theory for the stochastic theta method as follows.


Assume that the following three requirements are fulfilled.


  * For any >0 and x_0 ∈^d, there exists a constant R = R(,x_0) > 0 such that

    (|X_k^x_0| ≥ R) < ,   for any k ≥ 0.



  * For any >0 and any compact subset K of ^d, there exists a positive integer k^* = k^*(,K) such that

    (|X_k^x_0 - X_k^y_0|<) ≥ 1 - ,   for any  k ≥ k^*  and any (x_0,y_0) ∈ K  K.



  * For any >0, n ≥ 1 and any compact subset K of ^d, there exists a R = R(,n,K) > 0 such that

    (sup_0 ≤ k ≤ n |X_k^x_0| ≤ R) > 1 - ,    for any x_0 ∈ K.


Then the numerical solution generated by the stochastic theta method {X_k}_k ≥ 0 has a unique stationary distribution Π_ t.


Although the theory is very general, the conditions in it are in the sense of probability which are not easy to check. In this paper, we give some coefficients related conditions, i.e Conditions <ref> to <ref>, and prove the existence and uniqueness of the stationary distribution of the solution generated by the ST method under those conditions.






§ MAIN RESULTS


In this section, we present the main results of this paper. Since different choices of the parameter θ in (<ref>) require different requirements on the coefficients, f and g, we divide this section into three parts. We discuss the case when θ∈ [0,1/2) in Section <ref> and the situation when θ∈ [1/2,1] is presented in Section <ref>. The convergence of the numerical stationary distribution to the underlying counterpart is discussed in Section <ref>.





 §.§ θ∈ [0,1/2)
 






 Assume Conditions <ref> and <ref> hold, then for
 h<-(2μ+σ)/(1-θ)^2κ, the solution generated by the ST method (<ref>) obeys

 
    |X_k|^2≤ C_1,

 where C_1 is a constant that does not rely on k.


Applying Conditions <ref> and <ref>, we have

    |X_k+1|^2=   ⟨ X_k+1,X_k+(1-θ)f(X_k)h+g(X_k)Δ B_k⟩+θ h⟨ X_k+1,F(X_k+1)⟩
    ≤   1/2|X_k+1|^2+1/2[X_k+(1-θ)hf(X_k)+g(X_k)Δ B_K]^2+θ h μ|X_k+1|^2+aθ h
    
    =   1/2|X_k+1|^2+1/2[|X_k|^2+(1-θ)^2 f(X_k)^2 h^2+g(X_K)^2Δ B_K^2+2(1-θ)h X_kf(X_k)]
       +θ hμ|X_k+1|^2+aθ h+A_1,

where A_1=⟨ X_k,g(X_k)⟩Δ B_k+2(1-θ)h⟨ f(X_k),g(X_k)⟩Δ B_k.
Since 𝔼Δ B_k=0, we have 𝔼A_1=0. By iteration, we have

    𝔼|X_k+1|^2   ≤ A_2𝔼|X_k|^2+A_3
       ≤ A_2(A_2𝔼|X_k-1|^2+A_3)+A_3
       ≤ A_2^k+1𝔼|X_0|^2+A_3+A_2A_3+A_2^2A_3+⋯+A_2^k A_3,

where

    A_2=1+(1-θ)^2 h^2 κ+Δ B_k^2 σ+2(1-θ)hμ/1-2μθ h,

and

    A_3=(1-θ)^2 h^2c+bh+2(1-θ)ha+2aθ h/1-2μθ h.

Due to the facts h<-(2μ+σ)/(1-θ)^2κ, 1-2μθ h>0 and κ>0 and (<ref>),we have 0<A_2<1. This complete the proof.     



 Let Condition <ref> and <ref> hold. Then, for h<-(2K_2+K_1)/(1-θ)^2K_1 and any two initial values x,y∈ℝ^d with x ≠ y the solutions generated by the ST method (<ref>) satisfy
 
    |X_k^x-X_k^y|^2≤ C_3𝔼|x-y|^2,

 where

    C_3=[1+(1-θ)^2h^2K_1+hK_1+2K_2(1-θ)h/1-2K_2θ h]^k+1,

and lim_k→+∞C_3=0.


From (<ref>), we have

    |X_k+1^x-X_k+1^y|=   X_k^x-X_k^y+θ[f(X_k+1^x)-f(X_k+1^y)]h+(1-θ)[f(X_k^x)-f(X_k^y)]h
       +[g(X_k^x)-g(X_k^y)]Δ B_k.

Applying Conditions <ref> and <ref>, we have

    |X_k+1^x-X_k+1^y|^2=   ⟨ X_k+1^x-X_k+1^y,X_k^x-X_k^y+(1-θ)[f(X_k^x)-f(X_k^y)]h
       +[g(X_k^x)-g(X_k^y)]Δ B_k⟩+⟨ X_k+1^x-X_k+1^y,f(X_k+1^x)-f(X_k+1^y)⟩θ h
    ≤   1/2|X_k+1^x-X_k+1^y|^2+1/2[X_k^x-X_k^y+(1-θ)[f(X_k^x)-f(X_k^y)]h
       +[g(X_k^x)-g(X_k^y)]Δ B_k]^2+θ hK_2|X_k+1^x-X_k+1^y|^2.

Then, we have

    (1/2-K_2θ h )|X_k+1^x-X_k+1^y|^2≤   1/2|X_k^x-X_k^y|^2+1/2(1-θ)^2[f(X_k^x)-f(X_k^y)]^2h^2
       +1/2[g(X_k^x)-g(X_k^y)]^2Δ B_k^2
       +⟨ X_k^x-X_k^y,f(X_k^x)-f(X_k^y)⟩ (1-θ)h+Q_1,

where

    Q_1=⟨ X_k^x-X_k^y,g(X_k^x)-g(X_k^y)⟩Δ B_k+(1-θ)hΔ B_k⟨ f( X_k^x)-f(X_k^y),g(X_k^x)-g(X_k^y)⟩.

It is not difficult to show that

    |X_k+1^x-X_k+1^y|^2≤   1+(1-θ)^2h^2K_1+Δ B_k^2K_1+2K_2(1-θ)h/1-2K_2θ h|X_k^x-X_k^y|^2+Q_1/1-2K_2θ h.

Since 𝔼Δ B_k=0, we have 𝔼Q_1=0. Then, we obtain

    𝔼|X_k+1^x-X_k+1^y|^2≤C̅_3𝔼|X_k^x-X_k^y|^2.

By iteration, we have

    𝔼|X_k+1^x-X_k+1^y|^2≤C̅_3^k+1𝔼|x-y|^2,

where

    C̅_3=1+(1-θ)^2h^2K_1+hK_1+2K_2(1-θ)h/1-2K_2θ h.

Since h<-(2K_2+K_1)/(1-θ)^2K_1, 1-2K_2θ h>0 and (<ref>), we have
0<C̅_3<1.
This complete the proof. 




Given Conditions <ref> and <ref>, the solution generated by the ST method (<ref>) obeys

    ( sup_0 ≤ k ≤ n |X_k|^2 )≤ C_2,

where C_2 is a constant that can rely on k.


From (<ref>), we have

    |X_k+1|^2=   ⟨ X_k+(1-θ) f(X_k)h+g(X_k)Δ B_k,X_k+1⟩+⟨θ f(X_k)h, X_k+1⟩
    ≤   1/2|X_k+(1-θ) f(X_k)h +g(X_k)Δ B_k|^2 + 1/2|X_k+1|^2 + θ h(μ|X_k+1|^2+a).

Applying Conditions <ref> and <ref>, we get

    (1-1/2-θ h μ)|X_k+1|^2≤   1/2[|X_k|^2+(1-θ)^2 h^2 |f(X_k)|^2+|g(X_k)|^2|Δ B_k|^2
       + 2 h (1-θ)⟨ X_k,f(X_k)⟩+2 ⟨ X_k,g(X_k)⟩Δ B_k
       +2 (1-θ) h ⟨ f(X_k),g(X_k)⟩Δ B_k]+θ h a
    ≤   1/2[|X_k|^2+(1-θ)^2 h^2 (κ |X_k|^2+c)+|g(X_k)|^2|Δ B_k|^2
       +2 h(1-θ)(μ|X_k|^2 +a)+|X_k|^2+|g(X_k)|^2|Δ B_k|^2
       +Δ(1-θ)(|f(X_k)|^2+|g(X_k)|^2|Δ B_k|^2)]+θ h a
    ≤   1/2[(2+(1-θ)^2 h^2κ+2h(1-θ)μ+h(1-θ)κ)|X_k|^2
       +(2+h(1-θ))|g(X_k)|^2|Δ B_k|^2]+1/2h(1-θ)c
       +1/2ch^2(1-θ)^2+ah(1-θ)+θ h a.

That is

    |X_k+1|^2≤ D_1|X_k|^2+D_2|g(X_k)|^2|Δ B_k|^2+D_3,

where

    D_1=   [1+1/2κ h^2(1-θ)^2+hμ(1-θ)+1/2κ h(1-θ)]/(1/2-hμθ),
    
    D_2=   [1+1/2 h (1-θ)]/(1/2-h μθ),
    
    D_3=   [1/2h(1-θ)c+1/2ch^2(1-θ)^2+ah(1-θ)+ahθ ]/(1/2-hμθ).

Summarizing both sides yields

    ∑_i=1^k+1|X_i|^2=D_1∑_i=0^k|X_i|^2+D_2∑_i=0^k|g(X_i)|^2|Δ B_k|^2+(k+1)D_3.

Now we have

    |X_k+1|^2=(D_1-1)∑_i=0^k|X_i|^2+|X_0|^2+D_2∑_i=0^k|g(X_i)|^2|Δ B_k|^2+(k+1)D_3.

Taking the supreme and expectation on both sides gives

    ( sup_0 ≤ k ≤ n |X_k|^2 )
    ≤   (D_1-1)∑_i=0^k( sup_0 ≤ k ≤ n |X_i|^2 )
       +D_2( sup_0 ≤ k ≤ n(∑_i=0^k|g(X_i)|^2|Δ B_k|^2))+(k+1)D_3+|X_0|^2.

Then we have

    ( sup_0 ≤ k ≤ n |X_k|^2 )
    ≤   (D_1-1+D_2h^2σ)∑_i=0^k( sup_0 ≤ k ≤ n |X_k|^2 )
       +D_2h^2kb+(k+1)D_3+|X_0|^2,

where | B_k|^2 = h is used. Using the discrete version of the Gronwall inequality, Lemma <ref>, we have

    ( sup_0 ≤ k ≤ n |X_k|^2 ) ≤ (D_2h^2kb+(k+1)D_3+|X_0|^2)exp ((k+1)(D_1-1+D_2h^2σ)).

The proof is complete.    
Combining Lemmas <ref>, <ref> and <ref> and using Chebyshev's inequality, we derive the existence and uniqueness of the stationary distribution of the ST method with θ∈ [0, 1/2) from Theorem <ref>.



 §.§ θ∈ [1/2,1]
 

When θ∈ [1/2,1], we do not need the part for f(x) in Condition <ref> but only need that

    |g(x)-g(y)|^2≤ K_1 |x-y|^2,

for any x,y ∈^d.

To prove Lemmas <ref> and <ref>, let us present the following two lemmas and we refer the readers to <cit.> for the proof.


Let Condition <ref> hold, then for any β_1, β_2∈ℝ with β_2≥β_1≥0, the inequality

    |x-β_1f(x)|^2+2β_1a≤1-μβ_1/1-μβ_2(|x-β_2f(x)|^2+2β_2a)

holds.




Let Condition <ref> hold, then for any λ_1, λ_2∈ℝ with λ_2≥λ_1≥0, the inequality

    |x-y-λ_1[f(x)-f(y)]|≤1-K_2λ_1/1-K_2λ_2|x-y-λ_2[f(x)-f(y)]|

holds.

Now we are ready to present the three main lemmas in this subsection.


Given Condition <ref>, (<ref>) and (<ref>) hold, the solution generated by ST method (<ref>) obeys
 
    |X_k|^2≤ c_1,

 where c_1 is a constant that does not rely on k.


Denote

    θ^*=1+σ/4μ, λ=2μ+σ/2μ∧(2θ-1).

If θ∈[1/2,θ^*], by the definition of λ and Lemma <ref>, we have

    |X_k+1-(1-θ+λ)hf(X_k+1)|^2+2a(1-θ+λ)h
    ≤   |X_k+1-θ f(X_k+1)h+(2θ-1-λ)f(X_k+1)h|^2+2a(1-θ+λ)h
    ≤   |X_k+1-θ f(X_k+1)h|^2+2(2θ-1-λ)⟨ X_k+1,f(X_k+1)⟩ h
       +[(2θ-1-λ)^2-θ^2]|f(X_k+1)|^2h^2+2a(1-θ+λ)h
    ≤   |X_k-(1-θ)f(X_k)h|^2+4(1-θ)⟨ X_k,f(X_k)⟩ h+|g(X_k)|^2h
       +2(2θ-1-λ)⟨ X_k+1,f(X_k+1)⟩ h+2a(1-θ+λ)h+G_k
    ≤   |X_k-(1-θ)f(X_k)h|^2+2(1-θ)ah+4(1-θ)⟨ X_k,f(X_k)⟩ h+|g(X_k)|^2h
       +2(2θ-1-λ)⟨ X_k+1,f(X_k+1)⟩ h+2a(1-θ+λ)h-2(1-θ)ah+G_k
    ≤   1-μ(1-θ)h/1-μ(1-θ+λ)h(|X_k-(1-θ+λ)f(X_k)h|^2+2(1-θ+λ)ah)
       +4(1-θ)⟨ X_k,f(X_k)⟩ h+|g(X_k)|^2h+2(2θ-1-λ)⟨ X_k+1,f(X_k+1)⟩ h
       +2a(1-θ+λ)h-2(1-θ)ah+G_k,

where

    G_k=2⟨ X_k+(1-θ)hf(X_k),g(X_k)⟩Δ B_k+g(X_k)^2(Δ B_k^2-h).

Denote

    F_k=𝔼|X_k-(1-θ+λ)f(X_k)h|^2+2(1-θ+λ)ah,

and

    N_h(λ)=[1-μ(1-θ)h]/[1-μ(1-θ+λ)h].

It is clear that 𝔼G_k=0. By Condition <ref> and (<ref>), we have

    F_k+1≤    N_hF_k+[4(1-θ)μ h+σ h]𝔼|X_k|^2+2hμ(2θ-1-λ)𝔼|X_k+1|^2
       +(2a+b)h
    ≤   N_h^k+1[F_0-2(2θ-1-λ)μ h|X_0|^2]+ψ_λ(h)h∑_i=0^kN_h^k-i𝔼|X_i|^2
       +2(2θ-1-λ)μ h𝔼|X_k+1|^2+1/1-N_h(2a+b)h,

where

    ψ_λ(h)=4(1-θ)+σ+2N_hμ(2θ-1-λ).

For θ∈(1/2,θ^*], we have ψ_λ(h)<0, 2(2θ-1-λ)μ h≤0.

From Lemma <ref>, we have 𝔼|X_k|^2≤ F_k and 0<N_h(λ)<1. Therefore, we get

    𝔼|X_k|^2≤    N_h^k[F_0-2(2θ-1-λ)μ h|X_0|^2]+1/1-N_h(2a+b)h
    ≤   [F_0-2(2θ-1-λ)μ h|X_0|^2]+1/1-N_h(2a+b)h.

Denote

    c_1=[F_0-2(2θ-1-λ)μ h|X_0|^2]+1/1-N_h(2a+b)h,

then, we have

    |X_k|^2≤ c_1.

For θ∈(1+σ/(4μ),1), choosing λ'<λ sufficiently small such that ψ_λ'(h)<0 for any h>0, then using the same arguments above, we complete the proof.   



Given Condition <ref> and (<ref>) , for any two initial values x,y∈ℝ_d with x ≠ y the solutions generated by the ST method (<ref>) obey
 
    |X_k^x-X_k^y|^2≤ c_3

with lim_i→+∞c_3=0.


Denote

    θ^*=1+K_1/4K_2,  λ=2K_2+K_1/2K_2∧(2θ-1).

If θ∈[1/2,θ^*], by the definition of λ and Lemma <ref>, we have

    ||X_k+1^x-X_k+1^y|-(1-θ+λ)h|f(X_k+1^x)-f(X_k+1^y)||^2 
    ≤   ||X_k+1^x-X_k+1^y|-θ h|f(X_k+1^x)-f(X_k+1^y)||^2
       +2(2θ-1-λ)h⟨ X_k+1^x-X_k+1^y,f(X_k+1^x)-f(X_k+1^y)⟩
    ≤   ||X_k^x-X_k^y|-(1-θ)h|f(X_k^x)-f(X_k^y)||^2
       +4(1-θ)h⟨ X_k^x-X_k^y,f(X_k^x)-f(X_k^y)⟩
       +|g(X_k^x)-g(X_k^y)|^2h+2(2θ-1-λ)h⟨ X_k+1^x-X_k+1^y,f(X_k+1^x)-f(X_k+1^y)⟩+M_k
    ≤   |1-K_2(1-θ)h/1-K_2(1-θ+λ)h|^2||X_k^x-X_k^y|-(1-θ+λ)h|f(X_k^x)-f(X_k^y)||^2
       +4(1-θ)h⟨ X_k^x-X_k^y,f(X_k^x)-f(X_k^y)⟩+|g(X_k^x)-g(X_k^y)|^2h
       +2(2θ-1-λ)h⟨ X_k+1^x-X_k+1^y,f(X_k+1^x)-f(X_k+1^y)⟩+M_k,

where

    M_k=   2⟨|X_k^x-X_k^y|+(1-θ)h|f(X_k^x)-f(X_k^y)|,|g(X_k^x)-g(X_k^y)|⟩Δ B_k
       +|g(X_k^x)-g(X_k^y)|^2(|Δ B_k^2-h).

Denote

    W_k=𝔼||X_k^x-X_k^y|-(1-θ)h|f(X_k^x)-f(X_k^y)||^2,

and

    L_h(λ)=|(1-K_2(1-θ)h)/(1-K_2(1-θ+λ)h)|^2.

It is not hard to see that 𝔼M_k=0. By Condition <ref> and (<ref>), we have

    W_k+1≤    L_h W_k+[4(1-θ)K_2h+K_1 h]𝔼|X_k^x-X_k^y|^2+2K_2 h(2θ-1-λ)𝔼|X_k+1^x-X_k+1^y|^2
    ≤   L_h^k+1[A_0-2(2θ-1-λ)K_2 h𝔼|x-y|^2]+φ_λ(h)h∑_i=0^kL_h^k-i𝔼|X_i^x-X_i^y|^2,

where φ_λ^h=4(1-θ)K_2+K_1+2(2θ-1-λ)K_2L_h.

For θ∈[1/2,θ^*], we have φ_λ^h<0, 𝔼|X_k|^2≤ W_k, and |L_h|<1. Then, we get

    𝔼|X_k^x-X_k^y|^2≤ L_h^k[W_0-2(2θ-1-λ)K_2h|x-y|^2].

For θ∈(1+K_1/(4K_2),1), choosing λ'<λ sufficiently small such that ψ_λ'(h)<0 for any h>0 and using the same arguments above, we complete the proof. 

 
Assume that Conditions <ref>, (<ref>) and (<ref>) hold, then

    ( sup_0 ≤ k ≤ n |X_k|^2 )≤ c_2

where c_2 is a constant that can rely on k.


From (<ref>), we have

    |X_k+1 - θ f(X_k+1) h |^2 = |X_k - θ f(X_k)h + f(X_k)h + g(X_k) Δ B_k|^2.

Rewriting the right hand side, we have

    |X_k+1 - θ f(X_k+1) h |^2
       = |X_k - θ f(X_k) h |^2 + 2 ⟨ X_k, f(X_k)h ⟩ + (1 - 2 θ) |f(X_k)|^2 h^2 
         + |g(X_k)Δ B_k|^2 + 2/θ⟨ X_k, g(X_k)Δ B_k ⟩
         - 2(1 - θ)/θ⟨ X_k - θ f(X_k)h, g(X_k)Δ B_k ⟩.

Due to the fact that θ∈ [1/2,1] and Condition <ref>, we have

    |X_k+1 - θ f(X_k+1) h |^2
       ≤ |X_k - θ f(X_k) h |^2 + |g(X_k)Δ B_k|^2 + 2(μ |X_k|^2 + a)h 
         + 2/θ⟨ X_k, g(X_k)Δ B_k ⟩ +  2(1 - θ)/θ⟨ X_k - θ f(X_k)h, g(X_k)Δ B_k ⟩.

Summarising both sides yields

    |X_k+1 - θ f(X_k+1) h |^2
       ≤ |X_0 - θ f(X_0) h |^2 + ∑_i=0^k |g(X_i) Δ B_i|^2 + 2ah(k+1) + 2μ h ∑_i=0^k |X_i|^2 
         + 2/θ∑_i=0^k ⟨ X_k, g(X_k)Δ B_k ⟩ +  2(1 - θ)/θ∑_i=0^k ⟨ X_k - θ f(X_k)h, g(X_k)Δ B_k ⟩.

By the elementary inequality, it is not hard to see that

    ( sup_0 ≤ k ≤ n| ∑_i=0^k ⟨ X_i - θ f(X_i)h, g(X_i)Δ B_i ⟩| )
       ≤( ∑_i=0^n |  ⟨ X_i - θ f(X_i)h, g(X_i)Δ B_i ⟩| ) 
       ≤1/2∑_i=0^n | X_i - θ f(X_i)h |^2 + 1/2∑_i=0^n | g(X_i) Δ B_i |^2.

Using (<ref>) and the fact that |Δ B_i|^2 = h, we have

    ( sup_0 ≤ k ≤ n∑_i=0^k |g(X_i) Δ B_i|^2 ) ≤∑_i=0^n  (σ |X_i|^2 + b)h.

Applying the elementary inequality and (<ref>), we have

    ( sup_0 ≤ k ≤ n∑_i=0^k ⟨ X_i, g(X_i)Δ B_i ⟩)
       ≤∑_i=0^n | ⟨ X_i, g(X_i)Δ B_i ⟩| 
       ≤1/2∑_i=0^n  |X_i|^2 + ≤1/2∑_i=0^n  (σ |X_i|^2 + b)h.

Now, taking the expectation and the supreme on both sides of (<ref>) and using (<ref>), (<ref>), and (<ref>) we have

    ( sup_0 ≤ k ≤ n |X_k+1 - θ f(X_k+1) h |^2 )    ≤| X_0 - θ f(X_0) h |^2 + K ∑_i=0^n  |X_k|^2 
         +1/2∑_i=0^n | X_i - θ f(X_i)h |^2 + Kh,

where K is a generic constant.
Due to Condition <ref>, we have

    |X_k - θ f(X_k) h |^2    = |X_k|^2 - 2 θ h ⟨ X_k, f(X_k) ⟩ + θ^2 h^2 |f(X_k)|^2 
       ≥ |X_k|^2  - 2 θ h ( - μ |X_k|^2 - a   ) + θ^2 h^2 |f(X_k)|^2  
       ≥ (1 + 2 θ h μ) |X_k|^2 + 2 a θ h.

Applying the discrete version of the Gronwall inequality and (<ref>) to (<ref>), the assertion holds. 
Combining Lemmas <ref>, <ref> and <ref> and using Chebyshev's inequality, we derive the existence and uniqueness of the stationary distribution of the ST method with θ∈ [1/2, 1] from Theorem <ref>.



 §.§ The Convergence



Given Conditions <ref> to <ref>, the convergence of the numerical stationary distribution to the underlying stationary distribution is discussed in this subsection.
Recall that the probability measure induced by the numerical solution, X_k, is denoted by _k(·,·), similarly we denote the probability measure induced by the underlying solution,x(t), by _t(·,·).



Let Conditions <ref> to <ref> hold and fix any initial value x_0 ∈^d. Then, for any given T_1>0 and >0 there exists a sufficiently small t^* >0 such that

    d_(_k  t(x_0,·),_k(x_0,·)) <

provided that t <  t^* and k  t ≤ T_1.

The result can be derived from the finite time strong convergence of the ST method <cit.>.
Now we are ready to show that the numerical stationary distribution converges to the underlying stationary distribution as time step diminishes.

Given Conditions <ref> to <ref>, then

    lim_ t → 0 d_ (Π_ t(·), π(·)) = 0.



Fix any initial value x_0 ∈^d and set >0 to be an arbitrary real number. Due to the existence and uniqueness of the stationary distribution of the underlying equation, there exists a Θ^* >0 such that for any t > Θ^*

    d_ (_t(x_0,·),π(·)) < /3.

Similarly, by Theorem <ref>, there exists a pair of t^**>0 and Θ^** >0 such that

    d_ (_k(x_0,·),Π_ t(·)) < /3

for all t <  t^** and k  t > Θ^**.
Let Θ = max(Θ^*,Θ^**), from Lemma <ref> there exists a t^* such that for any t <  t^* and k  t < Θ + 1

    d_ (_k  t(x_0,·),_k(x_0,·)) < /3.

Therefore, for any t < min( t^*, t^**), set k = [Θ / t] + 1/ t, we see the assertion holds by the triangle inequality. 



§ SIMULATIONS



We present three numerical results in this section to demonstrate the theoretical results. The first one is a linear scale SDE with the true stationary distribution known. The second one is also a scale SDE but with the super-linear drift coefficient, of which the stationary distribution can be found by solving some ordinary differential equation. The third one is a two dimensional case.




    dx(t)=-α x(t)dt+σ dB(t)   on     t≥0.


Given any initial value X_0=x(0)∈, from (<ref>) we have

    X_k+1=X_k-αθ X_k+1 h -(1-θ)α X_k h+σΔ B_k.

This gives that X_k+1 is normally distributed with mean

    E(X_k+1)=(1-α h/1+αθ h)^k+1 x(0).

And the variance is

    Var(X_k+1)
       =(1+θα h)^-2/(1+θα h)^2(1-α h+θα h)^2Var(X_k)+σ^2h(1+θα h)^-2
       =σ^2h[(1+θα h)^-2+(1+θα h)^-4(1-α h+θα h)^2
       +(1+θα h)^-6(1-α h+θα h)^4+⋯ +(1+θα h)^-2(k+1)(1-α h+θα h)^2k]
       =1-[(1+θα h)^-2(k+1)(1-α h+θα h)^2(k+1)]/(1+αθ h)^2-((1-α h+θα h)^2
       =σ^2/2α-α^2h+2α^2θ h.

So the distribution of the solution generated by the ST method approaches the normal distribution N(0,σ^2/2α-α^2h+2α^2θ h) as k→∞.
Choosing α = σ =2, we draw several pictures. In this setting, the true stationary distribution is the standard normal distribution.




Figure <ref> shows the empirical density function of the numerical solution to (<ref>). Here the step size is chosen to be 0.001, the terminal time is 10 and the initial value is 2. 1000 sample paths with the θ = 1/2 are used to draw the graph. It can be seen that with the time advancing the density function is tending to an stable one, which indicates the existence of the stationary distribution.





Now, we use the Kolmogorov-Smirnov test (K-S test) <cit.> to measure the difference between the numerical stationary distribution and the true stationary distribution. Figure <ref> displays the changes in the p value as the time advances. It can be seen that after roughly t=1.8, the K-S test indicates that one can not reject that the samples generated by the numerical method are from the true distribution with the 95% confidence.





The convergent rate of the numerical stationary distribution with different choices of θ to the true stationary distribution is plotted in Figure <ref>. The step sizes, 2^-1, 2^-2, 2^-3, 2^-4, are used at T=10. It can be seen that the convergent rate is approximately one.



Next, we consider the SDE with the super-linear drift coefficient.



    dx(t)=-0.5(x(t)+x^3(t))dt+dB(t),

with x(0)=x_0.


The corresponding Kolmogorov-Fokker-Planck equation for the theoretical probability density function of the stationary distribution p(x) is

    0.5d^2 p(x)/dx^2 - d/dx (-0.5(x+x^3)p(x)) = 0.

And the exact solution is known to be <cit.>

    p(x) = 1/I_1/4(1/8)+I_-1/4(1/8)exp(1/8-1/2 x^2 - 1/4x^4),

where I_ν(x) is a modified Bessel function of the first kind.
Figure <ref> shows the changes of the  empirical density function with the time advancing. It can be seen that with the time variable increasing the center of the density function rapidly moves from the initial value, 2, to the theoretical centre one. And the density function is quite stable as time goes large.





Figure <ref> shows the convergent rate at T=10 with the step sizes, 2^-1, 2^-2, 2^-3, 2^-4. It can be seen that the rate for the super-linear case is not as good as the linear case, but the plots still show the convergence as the step size getting small.






Figure <ref> shows the change of the p value in the K-S test with the time increasing. It can be seen that after the time t=3 approximately one can not reject the numerical samples are from the true distribution with 95% confidence. This is a little bit worse than the linear case, as one need to wait a bit longer to see the stationary distribution.




We turn to the two dimensional super-linear SDEs.


    d[
    
    x_1(t)
    
    x_2(t)
    ]
    =[
    
    -x_1^3(t)-5x_1(t)+x_2(t)+5
    
    -x_2^3(t)-x_1(t)-5x_2(t)+5
    ]dt+[
    
    x_1(t)-x_2(t)+3
    
    -x_1(t)-x_2(t)+3
    ]dB(t),

with some initial data.


We check Condition <ref>, (<ref>) and (<ref>). Using the fundamental inequality that ab ≤ 1/2(a^2 + b^2), we have

    ⟨ x-y,   f(x) - f(y)⟩
       =(x_1-y_1,x_2-y_2)(
    
    f_1(x)-f_1(y)
    
    f_2(x)-f_2(y)
    )
       =(x_1(t)-y_1(t),x_2(t)-y_2(t))(
    
    -x_1^3(t)-5x_1(t)+x_2(t)+y_1^3(t)+5y_1(t)-y_2(t)
    
    -x_2^3(t)-5x_2(t)-x_1(t)+y_2^3(t)+5y_2(t)+y_1(t)
    )
       =((-x_1^3(t)+y_1^3(t))+5(y_1(t)-x_1(t))+(x_2(t)-y_2(t)))(x_1(t)-y_1(t))
            +((-x_2^3(t)+y_2^3(t))+5(y_2(t)-x_2(t))+(-x_1(t)+y_1(t)))(x_2(t)-y_2(t))
       =(-x_1^3(t)+y_1^3(t))(x_1(t)-y_1(t))+(-x_2^3(t)+y_2^3(t))(x_2(t)-y_2(t))
            -5((x_1(t)-y_1(t))^2+(x_2(t)-y_2(t))^2)
       =-(x_1(t)-y_1(t))^2(y_1^2(t)+x_1(t)y_1(t)+x_1^2(t))-(x_2(t)-y_2(t))^2(y_2^2(t)+x_2(t)y_2(t)+x_2^2(t))
            -5((x_1(t)-y_1(t))^2+(x_2(t)-y_2(t))^2).

Since

    y_1^2(t)+x_1(t)y_1(t)+x_1^2(t)
    ≥y_1^2(t)+2x_1(t)y_1(t)+x_1^2(t)/2≥(x_1(t)+y_1(t))^2/2≥0,

we have

    ⟨(x-y),(f(x) - f(y)) ⟩   ≤ ((x_1(t)-y_1(t))^2+(x_2(t)-y_2(t))^2)-5((x_1(t)-y_1(t))^2+(x_2(t)-y_2(t))^2)
       =-4((x_1(t)-y_1(t))^2+(x_2(t)-y_2(t))^2).

Also, we have

    |g(x)-g(y)|^2   =((x_1(t)-y_1(t))-(x_2(t)-y_2(t)))^2+((y_2(t)-x_2(t))-(x_1(t)-y_1(t)))^2
       =2((x_1(t)-y_1(t))^2+(x_2(t)-y_2(t))^2),

and

    2⟨(x-y),(f(x) - f(y)) ⟩ + |g(x)-g(y)|^2≤ -6((x_1(t)-y_1(t))^2+(x_2(t)-y_2(t))^2).




We plot the two dimensional empirical density function at different time. The initial values are [2,3]^T, the step size is 0.1 and 2 × 10^6 sample points are used to draw the plots. We can see from Figure <ref> that when the time is small the density function changes quite a lot even within a small time interval. But with time goes by, the density function stabilise to some certain shape, which could be regarded as the stationary distribution. Figure <ref> shows almost no difference between the empirical density functions at T=18 and T=20.








§ CONCLUSION
 
In this paper, we study the numerical stationary distributions generated by the stochastic theta methods. Both the drift and diffusion coefficients are required to satisfy the global Lipschitz condition when θ∈ [0, 1/2), but some super-linear terms are allowed to appear in the drift coefficient when θ∈ [1/2, 1]. Three numerical examples are given to show that the convergence and convergent rate of the numerical stationary distributions to their true counterparts. The plots also indicate that the numerical stationary distributions from the numerical solutions to SDEs could be used to approximate some non-linear deterministic differential equations.
















§ REFERENCES



10

abdulle2017optimal
A. Abdulle, I. Almuslimani, and G. Vilmart.
Optimal explicit stabilized integrator of weak order one for stiff
  and ergodic stochastic differential equations.
arXiv preprint arXiv:1708.08145, 2017.

BSY2016
J. Bao, J. Shao, and C. Yuan.
Approximation of invariant measures for regime-switching diffusions.
Potential Anal., 44(4):707–727, 2016.

BBKR2012a
G. Berkolaiko, E. Buckwar, C. Kelly, and A. Rodkina.
Almost sure asymptotic stability analysis of the
  θ-Maruyama method applied to a test system with stabilising and
  destabilising stochastic perturbations.
LMS J. Comput. Math., 15:71–83, 2012.

CW2012
L. Chen and F. Wu.
Almost sure exponential stability of the θ-method for
  stochastic differential equations.
Statist. Probab. Lett., 82(9):1669–1676, 2012.

Hig2000a
D. J. Higham.
Mean-square and asymptotic stability of the stochastic theta method.
SIAM J. Numer. Anal., 38(3):753–769 (electronic), 2000.

HJK2011a
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(2130):1563–1576, 2011.

IW1981a
N. Ikeda and S. Watanabe.
Stochastic differential equations and diffusion processes,
  volume 24 of North-Holland Mathematical Library.
North-Holland Publishing Co., Amsterdam, 1981.

LM2015
W. Liu and X. Mao.
Numerical stationary distribution and its convergence for nonlinear
  stochastic differential equations.
J. Comput. Appl. Math., 276:16–29, 2015.

M1991a
X. Mao.
Stability of stochastic differential equations with respect to
  semimartingales, volume 251 of Pitman Research Notes in Mathematics
  Series.
Longman Scientific & Technical, Harlow, 1991.

M2008a
X. Mao.
Stochastic differential equations and applications.
Horwood Publishing Limited, Chichester, second edition, 2008.

MS2013a
X. Mao and L. Szpruch.
Strong convergence rates for backward Euler-Maruyama method for
  non-linear dissipative-type stochastic differential equations with
  super-linear diffusion coefficients.
Stochastics, 85(1):144–171, 2013.

MYY2005a
X. Mao, C. Yuan, and G. Yin.
Numerical method for stationary distribution of stochastic
  differential equations with Markovian switching.
J. Comput. Appl. Math., 174(1):1–27, 2005.

M1951a
F. J. Massey J.
The kolmogorov-smirnov test for goodness of fit.
Journal of the American statistical Association,
  46(253):68–78, 1951.

QLH2014
Q. Qiu, W. Liu, and L. Hu.
Asymptotic moment boundedness of the stochastic theta method and its
  application for stochastic differential equations.
Adv. Difference Equ., pages 2014:310, 14, 2014.

Rod05
A. Rodkina and H. Schurz.
Almost sure asymptotic stability of drift-implicit θ-methods
  for bilinear ordinary stochastic differential equations in R^1.
J. Comput. Appl. Math., 180(1):13–31, 2005.

SOO1973
T. T. Soong.
Random differential equations in science and engineering.
Academic Press [Harcourt Brace Jovanovich Publishers], New York,
  1973.
Mathematics in Science and Engineering, Vol. 103.

talay1990second
D. Talay.
Second-order discretization schemes of stochastic differential
  systems for the computation of the invariant law.
Stochastics: An International Journal of Probability and
  Stochastic Processes, 29(1):13–36, 1990.

WWL2008
W. Wang, L. Wen, and S. Li.
Nonlinear stability of θ-methods for neutral differential
  equations in Banach space.
Appl. Math. Comput., 2008.

YM2003a
C. Yuan and X. Mao.
Asymptotic stability in distribution of stochastic differential
  equations with Markovian switching.
Stochastic Process. Appl., 103(2):277–291, 2003.

YM2004a
C. Yuan and X. Mao.
Stability in distribution of numerical solutions for stochastic
  differential equations.
Stochastic Anal. Appl., 22(5):1133–1150, 2004.

YM2005a
C. Yuan and X. Mao.
Stationary distributions of Euler-Maruyama-type stochastic
  difference equations with Markovian switching and their convergence.
J. Difference Equ. Appl., 11(1):29–48, 2005.

ZW2014
X. Zong and F. Wu.
Choice of θ and mean-square exponential stability in the
  stochastic theta method of stochastic differential equations.
J. Comput. Appl. Math., 255:837–847, 2014.

ZWH2015
X. Zong, F. Wu, and C. Huang.
Theta schemes for SDDEs with non-globally Lipschitz continuous
  coefficients.
J. Comput. Appl. Math., 278:258–277, 2015.




