{'fraction_non_alphanumeric': 0.045661506439443215, 'fraction_numerical': 0.020424092623910498, 'mean_word_length': 4.039331366764995, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stationary distribution of the stochastic theta method for nonlinear stochastic differential equations

The existence and uniqueness of the stationary distribution of the numerical solution generated by the stochastic theta method is studied. When the parameter theta 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.


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 [12,20,21], the authors studied the approximates to stationary distributions of SDEs and SDEs with Markovian switching by using the Euler-Maruyama method. In [2], 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 [6], 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 [8] for those SDEs with the super-linear drift coefficient. Higher order methods were also discussed for sampling the invariant measures [1,17].
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 [5], 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 [15]. The ST method was applied to a test system with stabilising and destabilising stochastic perturbations and almost sure asymptotic stability was analysed in [3]. The abilities to preserve the almost sure and the mean square exponential stabilities were discussed for different choices of the theta in [4] and [22], respectively. The asymptotic boundedness of the ST method was studied in [14]. 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 2. Section 3 contains the main results. Some numerical examples are used to illustrate the theoretical results in Section 4. We summarize the paper by Section 5.

Mathematical Preliminaries
In this paper, let (Ω, F , P) be a complete probability space with a filtration {F t } t≥0 satisfying the usual conditions that it is right continuous and increasing while F 0 contains all P-null sets. Let | · | denote the Euclidean norm in R 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 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 Now, we present the conditions on the coefficients.
Condition 2.2 Assume there exists a constant K 2 < 0 such that for any x, y ∈ R d In addition, we require that The next two conditions can be derived from Conditions 2.1 and 2.2 but with a little bit complicated coefficients. For the simplicity, we give two new conditions as follows.

Condition 2.3
There exist constants µ < 0 and a > 0 such that for any x ∈ R d x, f (x) ≤ µ|x| 2 + a.

Condition 2.4
There exist positive constants σ, κ, b and c such that for any In addition, we require that 2µ + σ < 0. (2.5) The existence and uniqueness of the underlying SDE (2.1) has been broadly studied. We refer the readers to Theorem 3.1 in [19] 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 2.11.
The stochastic theta method to SDE (2.1) is defined by is the Brownian motion increment and t k = kh, for The proof of the next lemma is similar to those in [8,11]. Proof. It is useful to write (2.6) as Define a function G : for θhK 2 < 1, we know that G has the inverse function G −1 : monotone. The ST method (2.6) can be written as Proof. If X k = x and X 0 = x, by (2.6) we see and 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 2.5, 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 2.8, we cite the following classical result (see, for example, Lemma 9.2 on page 87 of [10]).

Lemma 2.7
Let h(x, ω) be a scalar bounded measurable random function of x, independent of F s . Let ζ be an F s -measurable random variable. Then where H(x) = Eh(x, ω).
For any x ∈ R d and any Borel set B ⊂ R d , define Proof. The homogeneous property follows Lemma 2.6, so we only need to show the Markov property. Define The proof is complete.
Therefore, we see that P(·, ·) is the one-step transition probability and P 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, [9]). Denote the family of all probability measures on R d by P(R d ). Define by L the family of for any x, y ∈ R d . For P 1 , P 2 ∈ P(R d ), define metric d L by The weak convergence of probability measures can be illustrated in terms of metric d L [7]. That is, a sequence of probability measures {P k } k≥1 in P(R d ) converge weakly to a probability measure P ∈ P(R d ) if and only if lim k→∞ d L (P k , P) = 0.
Then we define the stationary distribution for {X k } k≥0 by using the concept of weak convergence.
Definition 2.10 For any initial value x ∈ R d and a given step size In [21], 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.
Theorem 2.11 Assume that the following three requirements are fulfilled.
• For any ε > 0 and • For any ε > 0 and any compact subset K of R d , there exists a positive integer • For any ε > 0, n ≥ 1 and any compact subset Then the numerical solution generated by the stochastic theta method {X k } k≥0 has a unique stationary distribution Π ∆t .

Remark 2.12
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 re-

Main Results
In this section, we present the main results of this paper. Since different choices of the parameter θ in (2.6) 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 3.1 and the situation when θ ∈ [1/2, 1] is presented in Section 3.2. The convergence of the numerical stationary distribution to the underlying counterpart is discussed in Section 3.3.
where C 1 is a constant that does not rely on k.

Proof. Applying Conditions 2.3 and 2.4, we have
where and any two initial values x, y ∈ R d with x = y the solutions generated by the ST method Applying Conditions 2.1 and 2.2, we have Then, we have where .
Since E∆B k = 0, we have EQ 1 = 0. Then, we obtain By iteration, we have where C 2 is a constant that can rely on k.
Proof. From (2.6), we have Applying Conditions 2.3 and 2.4, we get That is where

Now we have
Taking the supreme and expectation on both sides gives Then we have where E|∆B k | 2 = h is used. Using the discrete version of the Gronwall inequality, Lemma The proof is complete. Combining for any x, y ∈ R d .
To prove Lemmas 3.6 and 3.7, let us present the following two lemmas and we refer the readers to [18] for the proof. holds.
Now we are ready to present the three main lemmas in this subsection.
where c 1 is a constant that does not rely on k.

It is clear that EG k = 0. By Condition 2.3 and (2.3), we have
.
Proof. Denote If θ ∈ [1/2, θ * ], by the definition of λ and Lemma 3.5, we have It is not hard to see that EM k = 0. By Condition 2.2 and (3.1), we have 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. where c 2 is a constant that can rely on k.
Proof. From (2.6), we have Rewriting the right hand side, we have Due to the fact that θ ∈ [1/2, 1] and Condition 2.3, we have

Summarising both sides yields
By the elementary inequality, it is not hard to see that

Using (2.3) and the fact that E|∆B
Applying the elementary inequality and (3.4), we have Now, taking the expectation and the supreme on both sides of (3.2) and using (3.3), (3.4), where K is a generic constant. Due to Condition 2.3, we have Applying the discrete version of the Gronwall inequality and (3.7) to (3.6), the assertion holds.
Combining Lemmas 3.6, 3.7 and 3.8 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 2.11.

The Convergence
Given Conditions 2.1 to 2.4, 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 P k (·, ·), similarly we denote the probability measure induced by the underlying solution,x(t), byP t (·, ·). Lemma 3.9 Let Conditions 2.1 to 2.4 hold and fix any initial value x 0 ∈ R d . Then, for any given T 1 > 0 and ε > 0 there exists a sufficiently small ∆t * > 0 such that d L (P k∆t (x 0 , ·), P 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 [23].
Now we are ready to show that the numerical stationary distribution converges to the underlying stationary distribution as time step diminishes. Proof. Fix any initial value x 0 ∈ R 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 L (P t (x 0 , ·), π(·)) < ε/3.

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.
Given any initial value X 0 = x(0) ∈ R, from (2.6) we have This gives that X k+1 is normally distributed with mean

And the variance is
So the distribution of the solution generated by the ST method approaches the normal distribution N(0, σ 2 2α−α 2 h+2α 2 θh ) as k → ∞. Choosing α = σ = 2, we draw several pictures. In this setting, the true stationary distribution is the standard normal distribution. Figure 1 shows the empirical density function of the numerical solution to (4.1). 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) [13] to measure the difference between the numerical stationary distribution and the true stationary distribution. Figure  Figure 1: The empirical density function along the number of iterations 2 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 3. 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.
The corresponding Kolmogorov-Fokker-Planck equation for the theoretical probability density function of the stationary distribution p(x) is And the exact solution is known to be [16] p(x) = 1 I 1 where I ν (x) is a modified Bessel function of the first kind. Figure 4 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 5 shows the convergent rate at T = 10 with the step sizes, 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 6 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.
with some initial data.
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

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.