A brief review on stability investigations of numerical methods for systems of stochastic differential equations

1.   Introduction

Many authors have already investigated asymptotic stability of numerical methods of stochastic differential equations (SDEs). Numerical methods for approximating systems of SDEs experience a wide range of applications since most SDEs are not explicitly solvable. Thus, one has to be able to trust the qualitative behavior of those methods, which is adequate to that of underlying continuous time counterparts of SDEs. Convergence investigations may give some quantitative results as maximum step sizes $h$ of partitions tend to zero. However, those results do not answer which step size guarantees us adequate replication of qualitative behavior of exact solutions of SDEs. Moreover, stability forms one of the key pillars of modern axiomatic approaches of numerical analysis (for stochastic settings, cf. ^[51] to get to a more transparent understanding of this fact). This is why qualitative studies of long-term behavior of their numerical methods come into play and is of great importance. This paper is supposed to be a brief review on stability investigations of numerical methods for SDEs driven by standard Wiener processes, in terms of one-point motions.

In the present literature, the main focus is on the two criteria of asymptotic almost sure stability and mean square stability. Weaker stability criteria, such as stability in the mean sense, does not give us new insight compared to deterministic studies. Those two criteria are the most promising ones, which may give some new insights and which are also the most feasible ones. Thus, they have experienced a lot of attraction in the recent 40 years of research. Petersen ^[34] was one of the first who experimented with stochastic-numerical stability on computers. Nowadays, the qualitative attempts have gone far beyond those initial stages. A 3rd, but also important topic of stability investigations, is that of correct replication of stationary invariant probability distributions for systems with additive, state-independent noises. We shall only consider here the case of replication of distributions for finite-dimensional stationary Ornstein-Uhlenbeck processes under discretization (more general, the replication of distributions for stationary Langevin equations in infinite dimensions is another very valuable aim to pursue; see e.g. Cui et al. ^[7]). In this respect, our aim is to report on some of the most essential results and techniques of stability investigations for stochastic methods of discretization, which are applicable in every dimension $d\!\ge\! 1$ and not just characteristics of a specific numerical technique.

Without loss of generality, we shall stick to methods confined to Itô-calculus with the major assumptions of predictable integrands. Thus, we cannot easily make use of different implicit diffusion terms converging to other calculi solutions. However, under natural smoothness assumptions (such as coefficients in $C^1(\mathbb{R}^d)$), one can transform all standard stochastic calculi into each other (e.g., see Arnold ^[2]). Therefore, the restriction to Itô diffusion parts in this paper is not so confining as one thinks at the first glance. The big advantage of Itô calculus is that one can directly apply martingale theory and isometry relations---a situation which is totally different in locally anticipative calculi, such as the commonly used Stratonovich calculus in physics. Different stochastic calculi force us to exploit only changes of diffusion parts under numerical integration, which are adequate to the defined calculus parameters (note that all stochastic calculi may coincide in some special cases, such as for systems with state-independent noises). In general, that gives the user of the toolbox of numerical methods an additional restriction and complexity that is not apparent in deterministic settings. Thus, our presentation consistently follows the Itô interpretation in integrating the noise parts by explicit schemes (in contrast to implicit discretization of drift parts).

The following numerical methods for Itô SDEs (2.1) with drift $a$ and diffusion vectors $b^j$ ($j = 1, 2, ..., m$) are especially under inspection by this paper. The classic forward Euler-Maruyama method (FEM) has the scheme

along partitions $(t_n)_{n \in \mathbb{N}}$ with current step sizes $h_n = t_{n+1}-t_n$ and current noise increments $\Delta W^j_n$, which are independent and symmetrically distributed (in what follows, the noise increments $\Delta W^j_n$ need to be independent and must satisfy the finite moment identities

(at least). They can be Gaussian $\mathcal{N}(0, h_n)$, they could be symmetrically uniformly distributed $\mathcal{U}(-\sqrt{3h_n}, +\sqrt{3h_n})$; or, in case of path-wise approximations, they have to be $\Delta W^j_n = W^j(t_{n+1}) - W^j(t_n)$, or their symmetric truncation for Wiener processes $W = (W^1, W^2, ..., W^m)$). The backward drift-implicit Euler method (BEM) is of the form

The drift-implicit midpoint method (MM) possesses the representation

The drift-implicit trapezoidal methods (TM) follow

More generally, the class of (outer) drift-implicit theta methods (OTM) are given by

with implicitness $\theta \in \mathbb{R}^1$. The TMs (1.4) trivially represent a special case of OTMs (1.5) with $\theta = 0.5$. Inner drift-implicit theta methods (ITM) are governed by the scheme

with constant matrix parameter $\Theta \in \mathbb{R}^{d \times d}$ – a family introduced in ^[53] in a more general context. The midpoint methods (1.3) are a special case of this general class with $\Theta = \theta I_{d \times d}$ and $\theta = 0.5$. Note that OTMs (1.5) and ITMs (1.6) coincide for linear systems (that fact is obviously not the case for nonlinear systems). Finally, the (outer) drift-implicit theta Milstein methods (TMM) obey the scheme

with multiple Itô integrals

and the inner drift-implicit theta Milstein methods (ITMM) follow

with unit matrix $I_{d \times d} \in \mathbb{R}^{d \times d}\!.$ In passing, note that one can rewrite drift-implicit theta Milstein schemes (1.7) and (1.8) by using the stochastic concept of Levy areas, but that fact does not help in stability analysis much (just to improve the understanding of impact of commutative noise on simplification of scheme structure and when reduced complexity occurs, cf. ^[23]). Anyway, other higher order numerical methods involve the generation of multiple integrals $I_n^{j_1, j_2}$ and a large deal of complexity. Moreover, one can show that, for the asymptotic stability of them, one needs the asymptotic stability of equilibria of lower order parts in stochastic settings (see ^[43]). Thus, we put our main focus on stability of lower order numerical methods here. Another promising class for treating the problem of numerical stability is presented by balanced implicit methods (BIMs)

with uniformly bounded weights $c^j \ge 0$ (i.e., all $d \times d$-matrix-valued functions $c^j\!: \mathbb{R}^d \to \mathbb{R}^{d \times d}$) and $\Delta W^0_n: = h_n$ (see ^[29,44]). BIMs (1.9) are especially appropriate to treat the problem of a.s. asymptotic stability and form an easy way to implement implicit schemes due their linear-implicit character (note that a.s. = almost sure). Nonzero $c^0$ matrix-weights are needed for stabilization of moments and nonzero $c^j$ with $j = 1, 2, ..., m$ for a.s. stabilization, a.s. nonnegativity and path-wise boundedness on bounded domains (illustrated in ^[44] in detail). Another recent trend of stochastic numerical analysis is to combine drift-implicit theta methods (1.5) and (1.6) with balanced techniques (1.9)—a kind of double implicit method (first ideas in this direction are presented in ^[53]). This can extend the stability domains, i.e., the domains of parameters of SDEs for which the numerical method is asymptotically stable (in whatever probabilistic sense). We do believe that adequateness and asymptotic exactness in the sense of stability in every dimension $d \ge 1$ are more important topics than quantitative extensions illustrated on 1D complex domains with almost no practical relevance to multidimensional settings. For an overview on numerical methods for SDEs, see ^{[1,4,5,17,18,19,20,23,28,32,35,49,53,57,58]}. For general stochastic approximation principles and why stability is an important issue to control, see ^[51].

Our review is organized in 10 sections as follows. After introducing and presenting most commonly known stochastic numerical methods, Section 2 points out the open problem of relevance of classes of test equations. Section 3 discusses the notion of moment stability and stability functions, and its limitation to simultaneously diagonalizable SDEs. Section 4 is devoted to almost sure stability while exploiting classic dynamic and semi-martingale approaches. Asymptotic mean square stability of numerical methods for linear systems of SDEs is investigated by the help of stability operators in Sections 5 and 6. Section 7 treats the case of mean square stability while discretizing monotone, nonlinear test SDEs, which admit moment-dissipative dynamics (cf. ^[44,47,48]). Section 8 supplements the investigation with some comments on attempts to study stability of other numerical methods in dimensions 2D, 3D, etc. Section 9 delivers another approach to discuss stability and its preservation by discretizations of bilinear test systems of SDEs with state-independent noises (i.e., additive noise case). Eventually, Section 10 closes our contribution with some remarks, conclusions, and a brief summary.

2.   A general note on relevance of test equations

The problem of the relevant test equation for which criteria and class of SDEs has not been completely solved in large detail. Most of the authors follow directly the concept of linear test equations in one dimension. That cannot be considered as most relevant for general linear or nonlinear systems of SDEs in multi-dimensions. In passing, for theory of SDEs, see ^{[1,2,10,21,22,33,36]}, among many others. However, for $d$-dimensional systems of Itô SDEs

with smooth coefficients $a, b^j \in C^1(\mathbb{R}^d)$ and driven by independent, standard Wiener processes $W^j$, one could follow the following "localizing" argumentation: Suppose that SDEs (2.1) have a non-random equilibrium $x^* \in \mathbb{R}^d$ of the drift term such that

Then

where $J$ is the Jacobian matrix operator and $\xi^j$ are intermediate random values between $X(t)$ and $x^*$. Thus, by defining $Z(t) = X(t)-x^*$ and separating Eq (2.2) into systems with additive (i.e., state-independent) and multiplicative (i.e., state-dependent) noises, we gain the linearized systems (it has been our suggestion to separate into both for the sake of more transparency, cf. ^[44])

where $A = Ja(x^*)$, $B^j = Jb^j(x^*)$, and $b_0^j = b^j(x^*)$ (as $\xi^0(t) \approx x^*$ locally). The first case (2.3) is that for the situation $a(x^*) = b^j(x^*) = 0$, whereas the second case (2.4) is when $a$ and $b^j$ have no common equilibria with $b^j_0 = b^j(x^*)\not = 0$. Unfortunately, there are not many systems (2.3) that can be simultaneously diagonalized by (non-random) matrix-multiplication (for example, in the class of normal matrices $A, B^{(j)}$, we could continue with diagonalization). That would be a too restrictive assumption in general. At this point, the stability analysis bifurcates into the fairly simple 1D case (i.e., the simultaneously diagonalizable case of $A, B^{(j)}$ by a single matrix $L$) and the non-diagonalizable, but much more complex and richer case (i.e., the bilinear system (2.3) needs to serve as a relevant test equation). Most of the literature follows the diagonalizable case without giving any justification for such a 1D simplification (maybe by mimicking what is just common in the deterministic case). As we all know nowadays, stochastic integration is much more complex and offers a much greater variety. However, there is a first book ^[44] and publication ^[43] addressing the non-diagonalizable case and related numerical stability investigations in every dimension $d \ge 1$. The very restricted, diagonalizable case leads to the two types of test equations:

with complex constants $\lambda, \sigma \in \mathbb{C}$. However, these 2 cases need to be treated in truly different ways. In the first case, we have the trivial solution $x^* = 0$ as a stationary, deterministic equilibrium, whereas in the second case, the stationary solution is a random variable with all their moment characteristics and own nondegenerate, probability distribution (indeed Gaussian if $Y_0$ is deterministic or Gaussian). It turns out that the class of drift-implicit theta methods (1.5) is already able to provide adequate numerical solutions with correct asymptotic stability behavior for linear systems (2.3) and (2.4) in any dimension $d\ge 1$. This has been proven for $\theta = 0.5$ (i.e., for the midpoint and trapezoidal sub-case) in ^[45,46]. In fact, for asymptotic exactness, drift-implicit MMs (1.3) (i.e., theta methods with $\theta = 0.5$ or $\Theta = 0.5 I_{d \times d}$) with any non-random step size $h$ provide numerical solutions with adequate asymptotic behavior for linear SDEs (2.3) and (2.4) (i.e., no bias in its asymptotic moment characteristics and its distribution for any step size $h$ as integration time tends to infinity, hence equivalent to that of underlying SDE). They have indeed no asymptotic, probabilistic bias of any nature for any constant step size $h$. That is one of the reasons why we prefer to start with midpoint-type methods for numerical approximation of systems of SDEs and modify them to improve their qualitative behavior in stochastic settings. Other reasons are their deterministic contraction-exactness ^[54], deterministic energy-exactness ^[53], and deterministic symplectic character (the latter 2 topics are not subjects of this paper). Improved versions of drift-implicit midpoint methods can also be energy-exact for oscillators with additive noise (cf. ^[53]). Unfortunately, they themselves are not necessarily exact for all functionals under random perturbations. So, some care is needed, but they qualitatively discretize the deterministic parts of SDEs (i.e., under absence of noise terms) fairly well.

However, as we have already noted, linear stability analysis for 1D SDEs (2.5) and (2.6) has not been proven to be relevant for the general multidimensional case. The mathematical reason is that two or more matrices for the drift and the diffusion terms do not necessarily commute with each other and they cannot be simultaneously diagonalized in order to reduce the problem to 1D test equations. That is the main reason why one has to study fully multidimensional versions including the "non–commuting case". For this purpose, to bypass the problem of relevance of 1D test equations for linear systems, we introduce the concept of "stability operators" rather than the more traditional, but much easier concept of "stability functions". Indeed, this approach is feasible and delivers very assertive results on asymptotic stability, exploiting the positivity and monotonicity of related stability operators.

There is also an attempt to establish nonlinear test equations for numerical methods of SDEs. So far, there is not much known w.r.t. numerical stability in truly nonlinear, multidimensional situation, except for the treatment of moment-dissipative systems. A few exceptions in this respect are, for example ^{[15,16,44,48,50,54,59]}. However, the class of $p$-th moment dissipative, nonlinear SDEs forms one of the most promising and tractable classes. Uniformly $p$-th moment dissipative SDEs with drift $a$ and diffusion terms $b^j$ satisfy the one-sided boundedness condition (OSB)

with positive constant $c_{OSB} > 0$, where $\|\cdot \|_d$ denotes the Euclidean norm and $\langle\cdot, \cdot\rangle_d$ is the Euclidean scalar product of $\mathbb{R}^d$. Under this condition (2.7), one can show asymptotic $p$-th moment stability of trivial solution $x^* = 0$ for continuous time SDEs (2.1); see ^[47]. Most numerical methods for nonlinear SDEs have not been systematically investigated w.r.t. the preservation of asymptotic $p$-th moment stability of trivial equilibrium $x^* = 0 \in \mathbb{R}^d$ for all possible step sizes $h_n$ and any dimension $d\ge 1$. However, in ^[44], it is shown that the BEMs (1.2) guarantee that the (global) asymptotically $p$-th moment stable trivial equilibrium is preserved under their discretization for all admissible, nonrandom step sizes $h_n$ (i.e., for $p = 2$, mean square $A$-dissipativity, see Theorem 8.4.1, and for nonlinear mean square $A$-stability of BEMs (1.2) in any dimension $d \ge 1$, see Corollary 8.5.3 in chapter 8 of ^[44]). For equidistant step sizes $h$ and drift-implicit backward Euler method (1.2), see also ^[15,59], taylored to Hamiltonian systems. For exponential mean square stability of a split-step theta method with $\theta \in (0.5, 1]$ applied to nonlinear dissipative test equations with any equidistant, nonrandom step size $h$, see ^[16].

All in all, to the best of our knowledge, the problem of stochastic test equations for stochastic numerical methods has not been solved in any comprehensive manner so far. The best that we can do right now is to test stability of stochastic numerical methods with respect to nonlinear moment-dissipative systems of SDEs satisfying Eq (2.7). This suggestion was initiated by the book ^[44].

3.   Notion of moment stability and stability functions of some known numerical methods

As a mathematical-based paper, we make use of standard notations using quantifiers as common in mathematics and computer science in what follows. Below, $N(y)$ refers to the topologically relevant, but standard notion of (open) neighborhood of involved vector $y$.

Definition 3.1 (Local $p$-th Moment and Mean Square Stability). Let $y^* = x^*$ be an equilibrium of both numerical method $Y$ and SDE (2.1). A numerical method with values $Y = (Y_n(y_0))_{n \in \mathbb{N}}$ started with initial value $y_0$ and applied to SDEs (2.1) is said to be asymptotically $p$-th moment stable at $x^*$ or to have an asymptotically $p$-th moment stable equilibrium $y^* = x^*$ iff

where $N(y^*)$ denotes an open neighborhood of equilibrium $y^*$ and $Y_0(y_0) = y_0$. In the case of $p = 2$, we also speak of asymptotic mean square stability of $y^*$. If this property is true for all $y_0 \in \mathbb{D} \subseteq \mathbb{R}^d$, with $\mathbb{D}$ as the domain of definition of underlying SDE, then we make use of the phrase global asymptotically $p$-th moment stable equilibrium $y^*$.

For SDEs (2.5), one knows the mean square stability functions $R: \mathbb{C} \times \mathbb{C} \to \mathbb{C}$ of several numerical methods. The concept of those stability functions has been introduced by Saito and Mitsui ^[30,41,42] and studied by Hernandez and Spigler ^[12], Artemiev ^[3], Higham ^[13], and Schurz ^[43] at the very early stages. The explicit Euler-Maruyama method (1.1) for SDEs (2.5) possesses the mean square stability function $R_2^E$ satisfying

The (outer and inner) drift-implicit theta methods (1.5) and (1.6) for SDEs (2.5) have the mean square stability function $R_2^{(\theta)}$

The mean square behavior of (outer and inner) drift-implicit theta Milstein methods (1.7) and (1.8) for SDEs (2.5) is governed by the mean square stability function $R_2^{(\theta M)}$ satisfying

The advantage of analyzing stability functions $R: \mathbb{C} \times \mathbb{C} \to \mathbb{C}$ is only seen for illustration purposes with simple stability criteria (being relevant only for simultaneously diagonalizable SDEs) by involving just the modulus of $R$ with $u = \lambda h$, $v = \sigma \sqrt{h} \in \mathbb{C}$, such that

Additionally, this simplified approach is purely restricted to the case of equidistant step sizes $h$. Surely, it does not explain the case of variable step sizes $h_n$ (which would lead to the more complex, infinite Weierstrass products of complex analysis). Moreover, those functions are not proven to be relevant for the fully multidimensional, non-diagonalizable case! For multidimensional case, they have to be replaced by the concept of stability operators and the criteria governed by their spectral radii, as originally introduced in ^[43,44]. Furthermore, one can illustrate a problem with drift-implicit theta Milstein methods (1.7). Its higher order Milstein term leads to destabilization of mean square dynamics since the difference between stability functions of drift-implicit theta methods (1.5) and theta Milstein methods (1.7) consists of exclusively positive terms; hence, it gives a nondecreasing monotonicity relation for the drift-implicit theta Milstein methods (see Comparison Theorem 3.2 in ^[43]). This is the reason why we may suspect that the Dahlquist order barrier for mean square $A$-stable numerical methods could earlier be taken on with $r = 1.0$ or $r = 1.5$ in stochastics. Therefore, it gives us a good reason why one should prefer lower order numerical methods controlling both numerical stability and convergence in an adequate fashion, instead of higher order numerical methods with very restricted step sizes. Here, order just stands for order of local convergence, and not adequate long-term behavior in the sense of asymptotic stability.

The fact that we have to stay within Itô-calculus makes the adequate construction of correctly replicating numerical methods even harder. Recall that Itô calculus requires the predictability of stochastic integrands, and implicit numerical methods (except for linear-implicit or certain partial-implicit ones) do not maintain that predictability requirement for nonlinear SDEs. Thus, the discretization of diffusion terms "enjoys" much less degree of freedom for introducing implicitness. The same is true for exploitation of random variable step sizes $h_n$, which need to be predictable as well. Recall that predictability represents a requirement that the integrands $b^j(X(t))$ and their discrete replacements depend only on its past. So, these facts have limited all initiatives to introduce implicitness through stochastic diffusion terms (except for linear systems of SDEs discretized by balanced implicit methods (1.9) or other linear-implicit variants) while maintaining higher order of convergence.

4.   Almost sure stability of linear systems

Almost sure stability (also called T-stability by Saito and Mitsui ^[40]) is a truly probabilistic concept. We follow the standards, which are set by probability theory here. As commonly used, the abbreviation "a.s." stands for "almost surely".

Definition 4.1 (Local A.S. Asymptotic Stability or T-Stability). A numerical method with values $Y = (Y_n(y_0))_{n \in \mathbb{N}}$ with its equilibrium $y^* = x^*$ for SDEs (2.1) is said to be a.s. asymptotically stable at $x^*$ or have an a.s. asymptotically stable equilibrium $y^* = x^*$ iff

where $N(y^*)$ denotes an open neighborhood of equilibrium $y^*$ and $Y_0(y_0) = y_0$. The method $Y = (Y_n(y_0))_{n \in \mathbb{N}}$ is called a.s. A-stable iff its equilibrium $y^* = x^*$ is a.s. asymptotically stable whenever the SDEs (2.1) have an a.s. asymptotically stable equilibrium $x^*$.

Almost sure stability of linear systems of numerical methods in $\mathbb{R}^d$ can be studied by the following theorem, based on Kolmogorov's SLLN (Strong Law of Large Numbers) ^[55].

Theorem 4.2 (General Theorem on A.S. Asymptotic Stability of Numerical Methods ^[50]). Let $V = (V(n))_{n \in \mathbb{N}}$ be a sequence of nonnegative random variables $V(n) : (\Omega, {\mathcal F}_n, \mathbb{P}) \to \mathbb{R}^1_+$ with $V(0) > 0$ satisfying the recursive scheme

where $G(n) : (\Omega, {\mathcal F}_n, \mathbb{P}) \to \mathbb{R}^1_+$ are independent random variables such that there exists a nonrandom sequence $b = (b_n)_{n \in \mathbb{N}}$ with $b_n \to + \infty$ as $n \to + \infty$

Then, $V = (V(n))_{n \in \mathbb{N}}$ forms an a.s. (globally) asymptotically stable sequence, i.e., we have

Moreover, if

then $V = (V(n))_{n \to + \infty}$ is a.s. (globally) asymptotically unstable sequence, i.e., we have $\lim_{n \to + \infty} V(n) \; = \; +\infty$ (a.s.) for all nonrandom $y_0 \neq 0$.

Proof. See Schurz ^[50] based on Kolmogorov's SLLN ^[55] (i.e., Strong Law of Large Numbers). □

For application of Theorem 4.2, we may set $V(n) = \|Y_n\|_d$ or $V(n) = \|Y_n\|_d^2$ for dimension $d \ge 1$ and check conditions (4.2) or (4.3), resp. To illustrate the applicability of Theorem 4.2 in a fairly simple manner, the following linear test equation for almost sure stability is suggested by Itô SDEs

Its strong solution $X = (X(t))_{t \ge 0}$ forms a naturally induced martingale, which a.s. converges to the trivial equilibrium $x^* = 0$ whenever $\sigma \not = 0$. Then, the following result provides a mathematical evidence that the numerical experiments for BIMs (1.9) in ^[29] lead to the correct observation of numerical a.s. asymptotic stability of $x^* = 0$. It extends results that are found in ^[14,44,49].

Theorem 4.3 (A.s. A-Stability of a Subclass of BIMs (1.9) ^[50]). The BIMs $Y = (Y_n)_{n \in \mathbb{N}}$ governed by Eq (1.9) with scalar weights $c^0 = 0$ and $c^1 = |\sigma|$ applied to martingale test Eq (4.4) with any constant $\sigma \in \mathbb{R}^1\setminus \{0\}$ possess (globally) asymptotically a.s. stable trivial equilibrium $x^* = 0$ for any equidistant step size $h = \Delta t_n > 0$.

Proof. See Schurz ^[50]. It is based on an application of Theorem 4.2 with $V(n) = |Y_n|$ obeying the recursion $V(n+1) = G(n) V(n)$ (see Eq (4.1)) with independent, identically distributed random variables $G(n)$ satisfying

Thus, Theorem 4.2 with condition (4.2) and $b_n = n$ yields the claim of Theorem 4.3. □

Remark 4.4 (Existence of a.s. A-Stable Methods). Theorem 4.3 explains that the construction of a.s. $A$-stable (i.e., $T$-stable) numerical methods is indeed possible for linear SDEs, established constructively through an appropriate subclass of BIMs (1.9). Remarkably, there are not many known numerical methods which preserve a.s. asymptotic stability for SDEs under their discretization with any step size $h$!

There is another alternative to verify a.s. asymptotic stability of numerical methods while using semi-martingale convergence theorems (MCTs) (for MCTs, e.g., see Shiryaev ^[55] or Liptser and Shiryaev ^[56], as applied in Rodkina and Schurz ^[39,52]). To illustrate this idea, we consider drift-implicit theta methods (1.5) applied to linear systems of SDEs (2.3) and follow ideas borrowed from ^[39,52]. First, note that, for bilinear SDEs (2.3), their scheme with equidistant step sizes $h$, reduces to

provided that either matrix $A$ is negative semi-definite or the step sizes $h$ is sufficiently small such that inverses $(I-\theta h A)^{-1}$ exist. Let $O_{d \times d}$ be the matrix of $\mathbb{R}^{d \times d}$ with all zero entries and $\mathcal{S}_+^{d \times d}$ be the set of all positive semi-definite matrices of $\mathbb{R}^{d \times d}\!\!.$

Theorem 4.5 (A.s. Stability of Drift-implicit Theta Methods for Bilinear SDEs ^[52]). Assume that $Y_0 \in L^2(\Omega, \mathcal{F}, \mathbb{P})$ is independent from the $\sigma$-algebra of all underlying Wiener processes $W^j$, and there is a matrix $P \in \; \stackrel{\circ}{ \mathcal{S}_+^{d \times d}}$ such that $Q < O_{d\times d}$ where $I = I_{d\times d}$ and

in the sense of negative definite matrices (i.e., all eigenvalues of $Q$ are strictly negative). Then, the drift-implicit theta method (1.5) has a (global) a.s. asymptotically stable equilibrium $y^* = 0$.

Proof. Adapt the proof of Theorem 3.10 of ^[52] using the vector norm $V(n) = Y_n^T P Y_n$ instead of $Y_n^T Y_n$. For this purpose, use well-known semi-martingale convergence theorems (see MCTs in ^[55] or ^[56]) applied to the decomposition of

with $M_0 = Y_0$, $A_0^{(1)} = 0 = A_0^{(2)}$, where $M = (M_n)_{n \in \mathbb{N}}$ is the martingale part of $V = (V(n))_{n \in \mathbb{N}}$ and $A^{(k)} = (A_n^{(k)})_{n \in \mathbb{N}}$ are nondecreasing stochastic processes of finite variation. □

Remarkably, Theorem 4.5 supports the commonly believed intuition that, for sufficiently small step sizes, drift-implicit theta methods produce almost surely stable dynamics whenever the underlying bilinear SDE has an a.s. asymptotically stable equilibrium $x^* = 0$. This is definitely the case once the SDE additionally has an asymptotically mean square stable equilibrium $x^*$. However, for bilinear systems with pure diffusions (i.e., under absence of drift terms), those numerical methods may produce a.s. unstable dynamics. Then balanced terms as in BIMs (1.9) are needed, as we have seen above with 1D test Eq (4.4).

In passing, note that a significant simplification of investigations w.r.t. almost sure stability restricted to dimension $1$ has been suggested by the concept of T-stability for numerical solutions by Saito and Mitsui ^[40] at first.

5.   Mean square stability of linear systems by positive operators

In order to bypass the unsolved problem of a relevant reduction of multidimensional test systems of SDEs to 1D test equations, we suggest to follow the concept of positive stability operators (as first done in detail in ^[43,44]). In passing, we note that positive operators and their usefulness are explained in Krasnosel'skij et al. ^[24]. Let $\mathcal{S}_+^{d\times d}$ denote the space of positive-definite $d \times d$ matrices with nonrandom, real-valued entries. Recall the relation $S_{\!1} < S_{\!2}$ on the space $\mathcal{S}_+^{d \times d}$ of positive-definite matrices means that

Definition 5.1 (Mean Square Stability Operators ^[43,44]). The family of positive operators $\mathcal{L}_n: \mathcal{S}_+^{d\times d} \longrightarrow \mathcal{S}_+^{d\times d}$ such that

is called the family of (mean square) stability operators $(\mathcal{L}_n)_{n \in \mathbb{N}}$ of numerical method $Y = (Y_n)_{n \in \mathbb{N}}$.

Remark 5.2 (Meaning of Stability Operators). It actually suffices to restrict stability operators $\mathcal{L}_n$ to be defined on the class of rank-one matrices $S \in \mathcal{S}_+^{d\times d}$ since the rank of positive-definite matrices is kept under congruent matrix transforms. Moreover, under the operation of conditional expectations, this definition of stability operators $\mathcal{L}_n$ can be also rewritten as the requirement

Successive iterations of mean square stability operators $\mathcal{L}$ of numerical methods $Y = (Y_n)_{n \in \mathbb{N}}$ equivalently describe the temporal evolution of all second moments of random process $Y = (Y_n)_{n \in \mathbb{N}}$. That means that, for equidistant step sizes $h$, we find that the powers $\mathcal{L}^n$ of a single operator $\mathcal{L}$ obey the relationship

Hence, the actions and convergence of positive stability operators control the qualitative long-term behavior of the entire network of numerical dynamics in the mean square sense. Notice that, by the concept of trace of matrices, the connection to real-valued functionals is obvious (the "usual standard")

Theorem 5.3 (Mean Square Stability Operators for Drift-Implicit Theta Methods ^[43]). Assume that the initial data $Y_0 \in L^2(\Omega, \mathcal{F}, \mathbb{P})$ are independent of the $\sigma$-algebra generated by all Wiener processes $W^j$ and all step sizes $h_n > 0$ are nonrandom such that

Then, the mean square stability operators $\mathcal{L}_n$ associated to the drift-implicit theta methods (1.5) and (1.6) with $\Theta = \theta I_{d \times d}$ applied to linear test SDEs (2.3) with any nonrandom, current step sizes $h_n$ are determined by the mappings

where $I_{d\times d} \in \mathcal{S}_+^{d\times d}$ represents the $d \times d$ unit matrix, and, for $n \in \mathbb{N}$ and nonrandom step sizes $h_n$, we have

Proof. Simple exercise using basic properties of conditional expectations and linear algebra ^[43]. □

The mean square stability of numerical methods of bilinear SDEs (2.3) with equidistant step sizes $h$ can be decided by an application of the following simple theorem.

Theorem 5.4 (Sufficient Criterium for Mean Square Stability ^[43,44]). Assume that there is a matrix

for the stability operator belonging to numerical method $Y = (Y_n)_{n \in \mathbb{N}}$ applied to SDEs (2.3) with equidistant step sizes $h$ and equilibrium $y^* = 0$. Then, the numerical method $Y = (Y_n)_{n \in \mathbb{N}}$ has a (global) asymptotically mean square stable equilibrium $y^* = 0$.

Proof. Obviously, the associated stability operator $\mathcal{L}$ has the trivial fixed point $O_{d \times d}$ and the bilinear SDEs (2.3) have the trivial equilibrium $0$. Moreover, $\mathcal{L}$ is positive-definite and linear; hence, it is monotone on positive cone $\mathcal{S}_+^{d \times d}$. Under the hypothesis of Theorem 5.4, we know that

Actually, the Perron-Frobenius theorem and the theory of positive operators by Krasnosel'skij, Lifshits and Sobolev ^[24] initiated by Krein tells us that $\rho$ is connected to the spectral radius of operator $\mathcal{L}$, which is equal to the maximum of all eigenvalues of $\mathcal{L}$. Hence, thanks to the linearity of $\mathcal{L}$, we find that

since $0 < \rho < 1$. Obviously, we may conclude that

Therefore, by the relation to traces of matrices, we may conclude that

as $trace(vv^T) = \|v\|_d^2 = v^Tv = \langle v, v\rangle_d$ for all vectors $v \in \mathbb{R}^d$. This fact is equivalent to the asymptotic mean square stability of trivial equilibrium $y^* = 0$ for the associated numerical method. □

Theorem 5.4 can be applied to decide about mean square stability of several numerical methods when applied to linear test systems (2.3). We illustrate this fact by the following theorem on mean square stability for drift-implicit theta methods (1.5) and (1.6) with $\Theta = \theta I_{d \times d}$. For this purpose, from Arnold ^[2] (Corollary 11.4.14) and Khas'minskij ^[22], we recall that linear SDEs (2.3) have an asymptotically mean square stable equilibrium (and also exponentially mean square stable equilibrium) iff the Lyapunov identity

is satisfied.

Theorem 5.5 (Mean Square Stability of Drift-Implicit Theta Methods (1.5) ^[43]). Assume that mean square Lyapunov identity (5.1) is fulfilled, and the nonrandom and equidistant step size $h > 0$ is such that

Then, the drift-implicit theta methods (1.5) and (1.6) with $\Theta = \theta I_{d \times d}$, both with equidistant step size $h$ and implicitness $\theta \ge 0.5$ possess a (global) asymptotically mean square stable equilibrium $y^* = 0$ for Itô SDEs (2.3) in every dimension $d\ge 1$.

Proof. We borrow the ideas from ^[43,44]. Recall the form of associated mean square stability operator $\mathcal{L}$ related to drift-implicit theta methods (1.5) and (1.6) with $\Theta = \theta I_{d \times d}$ (both applied to linear systems (2.3)), which is determined by the claim of Theorem 5.3. Now, we begin with analyzing $\mathcal{L}$ by simultaneous multiplication of $\mathcal{L}(S) < S$ with matrices $(I_{d\times d}-\theta h A)$ from the left and $(I_{d\times d}-\theta h A)^T$ from the right (Recall that simultaneous linear transformations do not change the eigenvalues and monotonicity relations). This procedure gives

Now, take $S = P$ where $P \in \, \stackrel{\circ}{ \mathcal{S}_+^{d \times d}}$ is the matrix belonging to $Q \in \, \stackrel{\circ}{ \mathcal{S}_+^{d \times d}}$ in mean square Lyapunov identity (5.1). Consequently, we obtain that

since $\theta \ge 0.5$ and

Finally, apply Theorem 5.4 with $S^* = P$ to mean square stability operator $\mathcal{L}$ of drift-implicit theta methods (1.5) and (1.6) with $\Theta = \theta I_{d \times d}$ when applied to linear systems (2.3). This gives the conclusion on asymptotic mean square stability of equilibrium $y^* = 0$ in every dimension $d$, completing the proof of Theorem 5.5. □

Actually, the drift-implicit theta methods (1.5) and (1.6) with $\theta = 0.5$, $\Theta = 0.5 I_{d \times d}$ are the only methods of all drift-implicit theta methods with $\theta \in \mathbb{R}$, $\Theta = \theta I_{d \times d}$, which admit the equivalence of

This remarkable property of preservation of mean square stability is shown for bilinear systems (2.3) and (2.4) and every dimension $d \ge 1$ in ^[43,44,45].

Theorem 5.6 (Preservation of Mean Square Stability By Midpoint Methods (1.3) ^[45]). For all bilinear SDEs (2.3), the drift-implicit midpoint method (1.3) ($\theta = 0.5$, hence also trapezoidal) is the only drift-implicit theta method (1.5) that has an asymptotically mean square stable equilibrium $x^* = 0 = y^*$ iff Lyapunov identity (5.1) holds (i.e., the necessary and sufficient condition for exponential stability of trivial equilibrium $x^* = 0$ for underlying bilinear system (2.3)).

Proof. For details using fixed point theory of operators, see ^[45] or ^[44]. □

Remark 5.7 (Equivalence of Asymptotic and Exponential Stability for Bilinear Systems). In theory of dynamical systems, it is a commonly known fact that the notions of asymptotic stability and exponential stability are equivalent for autonomous, bilinear SDEs. Warning: This is not the case for nonlinear systems of SDEs in general! However, the proof techniques of those properties may significantly differ (even more in discrete time cases).

6.   Mean square stability of drift-implicit theta Milstein methods

The mean square stability of drift-implicit theta Milstein methods (1.7) and (1.8) applied to systems of linear SDEs (2.3) with equidistant, nonrandom step sizes $h > 0$ is understood by analyzing the associated mean square stability operator $\mathcal{L}^{(\theta M)}$ given as follows.

Theorem 6.1 (Mean Square Stability Operator for Drift-Implicit theta Milstein Methods ^[43]). Assume that the initial data $Y_0 \in L^2(\Omega, \mathcal{F}, \mathbb{P})$ are independent of the $\sigma$-algebra generated by all Wiener processes $W^j$ and all step sizes $h_n > 0$ are nonrandom such that

Then, the mean square stability operator $\mathcal{L}$ associated to the drift-implicit theta Milstein methods (1.7) and (1.8) with $\Theta = \theta I_{d \times d}$ applied to linear test SDEs (2.3) with any nonrandom, current step sizes $h_n$ is determined by the mappings $S \!\!\in\! \mathcal{S}_+^{d \times d} \!\!\mapsto\!\! \mathcal{L}_n^{(\theta M)}(S) =$

where $I_{d\times d} \in \mathcal{S}_+^{d\times d}$ represents the $d \times d$ unit matrix, and, for $n \in \mathbb{N}$ and nonrandom step sizes $h_n$, we have

Proof. For details, see ^[44]. □

Remark 6.2 (Comparison of Drift-Implicit $\theta$-Methods). For the same $h_n$ and same $\theta$, it is fairly easy to reckon that

Consequently, for mean square stability of trivial equilibrium of drift-implicit theta Milstein methods (1.7) (and methods (1.8)), the mean square stability of trivial equilibrium of drift-implicit theta methods (1.5) (and methods (1.6), resp.) is needed. Thus, theta Milstein methods (1.7) and (1.8) do not improve the stability behavior of numerical solutions, compared to their lower order counterparts. This crucial fact has been proved in ^[43] first. Hence, one needs to find more efficient, higher order, asymptotically stable discretization methods for all step sizes. This observation is true in all dimensions $d\ge 1$! Do higher order, asymptotically stable stochastic numerical methods really exist unconditionally for all step sizes and in all dimensions $d \ge 1$ (with asymptotic stability behavior equivalent to underlying SDE)? This represents an unanswered question on stochastic Dahlquist-type barriers and the existence of higher order, mean square $A$-stable, mean square preserving stochastic numerical methods (i.e., higher order than drift-implicit midpoint methods).

Theorem 6.3 (Conditional Mean Square Stability of Drift-Implicit Theta Milstein Methods ^[43]). Assume that the initial data $Y_0 \in L^2(\Omega, \mathcal{F}, \mathbb{P})$ are independent of the $\sigma$-algebra generated by all Wiener processes $W^j$ and all step sizes $h_n > 0$ are nonrandom such that

and there is a matrix $S^* \in \, \stackrel{\circ}{ \mathcal{S}}_+^{d \times d}$ such that the mean square stability operator $\mathcal{L}^{(\theta M)}$ of theta Milstein methods (1.7) with step size $h > 0$ satisfies

Then, drift-implicit theta Milstein methods (1.7) and (1.8) with $\Theta = \theta I_{d \times d}$ applied to bilinear SDE (2.3) with that nonrandom, equidistant step size $h$ are asymptotically mean square stable (hence, they are also exponentially mean square stable).

Proof. Apply Theorem 5.4 to stability operator $\mathcal{L}^{(\theta M)}$ (similarly as done in ^[43,44]). □

Remark 6.4. Theorem 6.3 is applicable for mean square stable, bilinear systems of SDEs (2.3) for only sufficiently small step sizes $h > 0$. This is due to the destabilizing effect of positive higher order Milstein terms in their mean square stability operator (cf. Theorem 6.1) and apparent for systems with very large noise terms. However, it is indeed applicable whenever Lyapunov identity (5.1) is valid in any dimension $d \ge 1$. It remains the task to find new variants of implementations of drift-implicit Milstein methods, which guarantee mean square stability of them for all situations when Eq (5.1) is satisfied and all possible step sizes $h > 0$ (i.e., to construct mean square A-stable, higher order methods preserving mean square stability for all step sizes iff relation (5.1) holds (i.e., a truly higher order, mean square stability preservation)).

7.   Mean square stability for some drift-monotone, nonlinear test systems

Consider the test class of Itô SDEs (2.1) with monotonically decreasing drift $a$ and sublinearly bounded diffusion terms $b^j$ with trivial equilibrium $x^* = 0$. Suppose that

Using Dynkin's formula with $(i)$–$(iii)$, it is fairly easy to show that all strong solutions $X$ of SDEs (2.1) have mean square exponentially stable trivial equilibria and they satisfy

Such systems belong to the more general class of mean square dissipative SDEs with $-2K_a+K_b \le 0$.

Theorem 7.1 (MS-Stability of Drift-Implicit Theta Methods (1.5) for MS-Dissipative SDEs ^[53]). Assume that $a(0) = 0$, $a \in C^0(N(0))$ on a neighborhood $N(0)$ of $x^* = 0$, (i)–(iii), $h > 0$ is nonrandom. Then, the drift-implicit theta methods (1.5) with $\theta \ge 0.5$ have a (local) asymptotically mean square stable equilibrium $y^* = 0$ when applied to SDEs (2.1) (MS = Mean Square).

Proof. This is part of a more general result of Theorem 1.6.3 with non-equidistant step sizes in ^[53], while verifying the monotonicity of expected Lyapunov functional

along the numerical dynamics associated to Eq (1.5), which are applied to mean square dissipative SDEs (2.1) with implicitness $\theta \ge 0.5$ and nonrandom step sizes $h > 0$. □

In fact, for drift-implicit backward Euler methods (1.2), we may establish exponential mean square stability of $y^*\! = \!0$ for mean square dissipative SDEs (2.1). This is content of the following theorem.

Theorem 7.2 (Exponential MS-Stability of Drift-Implicit Backward Euler Methods for MS-Dissipative SDEs ^[43,48,53]). Assume that $a(0)\! = \!0$, $a \!\in\! C^0(\mathbb{R}^d)$ on a neighborhood $N(0)$ of $x^*\!\! = \!0$, (i)–(iii), and $h\!\! > \!\!0$ is nonrandom. Then, drift-implicit backward Euler methods (1.2) have a (local) exponentially mean square stable equilibrium $y^*\!\! = \!0$ when applied to mean square dissipative SDEs (2.1).

Proof. This is part of a more general result of Theorem 1.3.5 referring to non-equidistant step sizes in ^[53], while verifying the inequality

along the numerical dynamics associated to Eq (1.5), which are applied to mean square dissipative SDEs (2.1) with implicitness $\theta \ge 0.5$ and nonrandom step sizes $h > 0$ in any dimension $d \ge 1$. Thus, exponential mean square stability of trivial equilibrium $y^* = 0$ of drift-implicit backward Euler methods is confirmed by complete induction on $n\in \mathbb{N}$. □

Remark 7.3 (On Mean Square Dissipative Systems). For mean square dissipative systems obeying (ⅱ), the diffusion functions $b^j$ also vanish at $x^*\! = \!0$, i.e., $b^j(x^*)\! = \!0$. So, we do not need to impose an extra condition on equilibria of $b^j$. Mean square dissipative numerical methods are classified in larger detail in ^[48], where one finds a proof of exponential mean square stability for Eq (1.2), including an estimation of exponential growth rates with variable, but nonrandom step sizes $h_n$ (cf. Corollary 5.8 in ^[48]).

Remark 7.4 (The Lack of Stability Investigations for Higher Order Methods). So far, higher order numerical methods (except for drift-implicit theta Milstein methods (1.7)) have not been investigated w.r.t. mean square stability systematically. Classic drift-implicit theta Milstein methods (1.7) do not improve the numerical stability behavior compared to their drift-implicit theta Euler counterparts (1.5) as shown in ^[44] (same for inner Milstein methods (1.8) and lower order inner theta methods (1.6), resp.). However, recent new versions of them such as very recently developed split-step (^[9], tamed (^[26,61,64]), improved (^[63]), or truncated balanced implicit Milstein methods (^[37]) may give promising results w.r.t. nonlinear stability under additional conditions

to control their higher order terms, where $b^{j_1}_i$ is the $i$-th component of vector $b^{j_1}$. In fact, this condition is implied by the additional hypothesis that

which is needed for proving strong convergence of higher order numerical methods.

8.   Some comments on stability analysis for other methods in 2D, 3D, ...

There are also attempts to investigate stability of other stochastic numerical methods. For example, mean square stability of two step methods are studied in D'Ambrosio et al. ^[8]. They focus on preserving long-term statistics related to the dynamics of a linear stochastic damped oscillator whose velocity is Gaussian in the stationary regime. Burrage and Lythe ^[6] carry out an analysis of the stationary density in the case of constant-coefficient 2D linear SDEs, imposing exact stationary statistics in the position variable and the absence of correlation between position and velocity of stochastic oscillations. They study explicit partitioned and two-stage "leapfrog" numerical methods (reverse and Runge-Kutta), which have good properties in the position variable and are symplectic in the limit of zero damping, sharing the midpoint property

Saito and Mitsui ^[42] investigate numerical mean square stability for 2D systems beyond single steps as well. Huang ^[16] claims that drift-implicit theta methods (1.5) with any implicitness $\theta \in [0, 1]$ preserves exponential stability of underlying linear systems (2.3) for sufficiently small step sizes $h$ (see Theorem 3.2, but that is no surprise since it follows from mean square convergence already while arguing with axiomatic principles ^[51]). Moreover, he shows that a split-step theta method (SST) with implicitness $\theta \in [0.5, 1]$

does preserve exponential stability of trivial equilibrium $x^* = 0$ for SDEs (2.3) with drift function $a(t, x)$ and single diffusion function $b(t, x)$ (i.e., $m = 1$), while using any step size $h > 0$. That split-step theta method (SST) can be also considered as a deterministic implicit predictor-corrector method. Exponential mean square stability of split-step backward Euler methods (i.e., the special case when $\theta = 1.0$) was shown in ^[15] earlier. Mean square stability of drifting split-step theta Milstein methods

is considered for 1D SDEs with $m = 1$ (single noise terms) in Eissa et al. ^[9]. Reshniak et al. ^[38] investigated mean square stability of several split-step composite and Adams-Moulton-type Milstein methods and Rosenbrock improvements. Stability matrices with their spectral radii are computed for multidimensional noise terms. Additionally, based on the well-known Routh–Hurwitz theorem, Tocino and Senosiain ^[60] offer a particular criterion yielding mean square stability of theta methods for two-dimensional systems in terms of their coefficients.

Of course, $L^2$-convergence of numerical methods implies that exponential stability can be preserved for sufficiently small step sizes $h$. However, it is rather computationally difficult to provide very efficient estimates on the maximum threshold $h^*$ for step sizes $h < h^*$ in order to preserve exponential stability while using equidistant $h < h^*$ for every dimension $d > 1$ (i.e., the only simple case is $d = 1$). So, constructive, easily implementable, and unconditional results on numerical stability for all step sizes $h$ are important to be found. In fact, asymptotically mean square exact numerical methods are available for bilinear systems (2.3) and (2.4). Such methods are given by midpoint-type methods and their improved modifications, as proven in ^[45,46,53]. Eventually, moment stability exponents for fairly general, implicit stochastic iterations with variable step sizes $h_n$ are estimated in ^[47] in order to measure the "strength" of attraction in the sense of $p$-th moment evolution. There one finds general estimates on the rate of convergence of drift-implicit split-type numerical methods to stable equilibria.

9.   Numerical preservation of long-term statistical quantities for SDEs with additive noise

The stability analysis is somewhat different for systems of bilinear SDEs (2.4) with state-independent diffusion terms, also called the SDEs of Ornstein-Uhlenbeck processes or additive noise with linear drift. Of course, we could ignore the effect of random noise terms by taking just expectations and conducting the same steps as with deterministic ODEs (Ordinary Differential Equations). Such an approach would not give any new insights into stability of stochastic numerical methods, as conducted in ^[11,23]. Hence, a truly probabilistic approach requires a slight change of definitions of stability with additive noise. Talay ^[59] has already been able to prove convergence of the backward Euler scheme (1.2) to the invariant measure when applied to stochastic Hamiltonian systems. A first classification of attraction to stationary probabilistic laws is found in ^[46]. We have found in ^[46] that MMs (1.3) for linear systems (2.4) (hence, also trapezoidal ones, and both inner and outer midpoint methods with $\Theta = 0.5 I_{d \times d}$ and $\theta = 0.5$) have no asymptotic bias in its invariant measure for any nonrandom, equidistant step size $h > 0$. Hence, this mathematical evidence suggests to measure the distance to related statistical quantities of underlying invariant stationary probability laws in the case of additive noise. This gives the motivation to rather make use of the following definitions on preservation of exact stationary values (which also provides a classification for systems with additive noise in general).

Definition 9.1 (Classification of Asymptotically Exact Stochastic Numerical Methods ^[44,46]). Let $X = \{X(t), t \ge 0\} \subseteq \mathbb{R}^d$ be a stationary, ergodic stochastic process and $\|\cdot\|_d$ be any vector norm of $\mathbb{R}^d$. Then, the numerical method $Y = (Y_n)_{n \in \mathbb{N}}$ discretizing stochastic process $X$ is said to be (asymptotically) $p$-th mean preserving $(p \in \mathbb{R}^1)$ iff

Furthermore, $Y$ is called (asymptotically) mean preserving iff

(asymptotically) mean square preserving iff

(asymptotically) variance preserving iff

and (asymptotically) in-law preserving iff

where $\mathcal{L} (.)$ denotes the probability law of the corresponding random variable.

This definition presents different concepts in general. However, in special cases such as preservation of Gaussian distributions, the two properties of mean and mean square (or variance) preservation of Gaussian distributed, numerical methods implies in-law preservation of Gaussian invariant distributions (Gaussian measures). We think of possibly more general cases than Gaussianity and SDEs (2.4) that do not have to possess Gaussian solutions for non-Gaussian initial data or non-Gaussian random coefficients. Moreover, some stochastic numerical methods may easily preserve moments in first moment sense, but not necessarily other $p$-th moments for any step size $h$. Such methods are widely used FEMs (1.1) and BEMs (1.2) with $p\! = \!2$.

For simplicity of this paper, we shall only refer to the Euclidean vector norm $\|\cdot\|_d$ of $\mathbb{R}^d$. The following proposition ^[45,46] establishes the superiority of MMs (1.3) w.r.t. preservation of asymptotic probabilistic laws of Gaussian-distributed SDEs (2.4), including of all of their statistical quantities as time tends to infinity.

Proposition 9.2 (Preservation of Probability Law by Midpoint Methods ^[45,46]). There are only the drift-implicit theta methods (1.5) and (1.6) with $\Theta = 0.5 I_{d \times d}$ among all drift-implicit theta methods with $\Theta = \theta I_{d \times d}, \theta \in \mathbb{R}^1$, which exactly replicate the asymptotic behavior of stationary Ornstein-Uhlenbeck processes (2.4) for all step sizes. More precisely, $Y = (Y_n)_{n \in \mathbb{N}}$ generated by schemes (1.5) with any equidistant, nonrandom step size $h > 0$ and implicitness $\theta = 0.5$ is asymptotically mean, $p$-th mean, mean square, variance and in-law preserving for stationary SDEs (2.4) with nonrandom initial data or Gaussian initial data (independent from all underlying Wiener processes $W^j$).

Proof. For general proof based on fixed point principles, see ^[45]. For general proof, based on diagonalization of drift matrices $A$, see ^[46]. It is obvious that the stationary law of drift-implicit theta methods (1.5) and (1.6) with $\Theta = \theta I_{d \times d}$, nonrandom implicitness $\theta$, nonrandom step sizes $h$ and Gaussian, centered increments $\Delta W_n$ is again Gaussian. Moreover, all eigenvalues of drift matrices $A$ must have negative real parts in order to have a finite, stationary law for linear SDEs (2.4). For simplicity of this presentation, consider the diagonalized, 1D form Eq (2.6)

of SDEs (2.4) with real parameters $\lambda > 0$ and $\sigma > 0$. Asymptotic mean preservation of theta methods with any nonrandom $\theta\ge 0.5$ and any nonrandom step sizes $h > 0$ is clear since it boils down to the well-known investigations of deterministic case. So, it remains to investigate the case for $p = 2$ since the preservation of all other moments follows from the characterization of Gaussian distributions. The drift-implicit theta methods (1.5) and (1.6) with $\Theta = \theta I_{d \times d}$ applied to Ornstein-Uhlenbeck processes (2.4) is governed by the explicit schemes

Thus, their 2nd moments obey the identity

Therefore, we arrive at the nonrandom difference equation

This difference equation can only have finite limits for all $\lambda > 0$ and all step sizes $h > 0$ iff $\theta \ge 0.5$. In this case, we have $0 < a < 1$ and, as a simple calculus exercise, the limit is computed by

Obviously, this generates a statistical bias for any step size $h$ whenever $\theta \not = 0.5$. Moreover, it is exact for all step sizes $h > 0$ and all parameters $\theta > 0$ iff $\theta = 0.5$. Since Gaussian distributions are uniquely characterized by 1st and 2nd moments, the entire stationary probability law must be the same as that of SDE (2.6) when $\theta = 0.5$, including all other $p$-th moments. This ends the proof of proposition 9.2. □

Remark 9.3 (Asymptotically In-Law Preserving Stochastic Numerical Methods). Up to now, there are not any other known (asymptotically) in-law preserving numerical methods for systems of SDEs (2.1) (i.e., for general problems and all step sizes $h$), except for the drift-implicit midpoint methods applied to bilinear SDEs (2.3) and (2.4) for every dimension $d\ge 1$. Note that BIMs (1.9) with appropriate weights are also asymptotically exact in the almost sure sense for the test class of pure, linear diffusions with multiplicative noise. The latter is just proven in dimension $d = 1$, but it can be extended to pure diffusions with positive matrices $B$ for $d > 1$. Of course, all theta methods (1.5) with $\theta \ge 0.5$ are asymptotically exact in mean sense for all nonrandom step sizes $h > 0$. However, asymptotic mean exactness is a very weak notion of exactness. Hence, it does not get a lot of attention in stochastics.

Remark 9.4 (The Advantage and Coincidence of Drift-Implicit Theta Methods for Bilinear SDEs (2.4) With Additive Gaussian Noise). One is also able to compute the characteristics of probability distributions associated to other drift-implicit methods when applied to bilinear SDEs (2.4), including its limits. However, note that drift-implicit theta methods (1.5) and (1.7) (and methods (1.6) and (1.8)) coincide for SDEs (2.4) with additive noise since higher order Milstein terms vanish for diffusions with additive noise. Therefore, we also know about the preservation of asymptotic laws of ergodic SDEs (2.4) by drift-implicit midpoint Milstein methods (i.e., methods (1.7) and (1.8) with $\theta = 0.5$, $\Theta = 0.5 I_{d \times d}$).

Remark 9.5 (Computation of Path-wise Limits for Additive Noise). The drift-implicit theta methods (1.5)–(1.8) applied to linear SDEs (9.1) with additive noise for any equidistant step size $h > 0$ have an explicit representation:

This reduces to Gaussian distributed random variables whenever $Y_0 = y_0$ is nonrandom or Gaussian itself independent from the $\sigma$-algebra induced by $W$, and the computation of Gaussian statistics is very well-known. In that case, path-wise limits of Eq (9.2) can be computed for nonrandom step sizes $h > 0$.

10.   Summary and conclusions

Higher order convergent numerical methods in general do not guarantee better stochastic stability properties. This could already been seen by the families of drift-implicit theta Euler methods compared to drift-implicit theta Milstein methods (cf. ^[43]). The successful construction of higher order, mean square $A$-stable stochastic numerical methods for systems (2.3) with multiplicative noise is still an open problem (i.e., which are mean square stable iff the underlying SDE has a mean square stable equilibrium). Stability and preservation of qualitative properties are two of the key pillars to understanding the qualitative behavior of numerical methods better. While our main focus is on discretization with equidistant step sizes in this paper, many of the major stability results can also be rewritten for variable, but nonrandom step sizes (as initialized by ^[44]). General stability investigations apply to both weakly and strongly convergent approximations (and $L^p$-convergent ^[17], ^[19] ones).

To bypass the problem of adequate test equations for mean square stability, we suggest to exploit the theory of positive operators and study the family of stability operators related to underlying numerical methods in any dimension $d \ge 1$. Those studies can be carried out on a system level using the property of monotonicity of positive operators and computations of spectral radii. For nonlinear, dissipative test equations with monotonicity conditions, we exploit appropriate Lyapunov functionals.

The problem of adequateness of test equations remains unsolved (i.e., the complete solution), except for SDEs with simultaneously diagonazible linearizations and smooth coefficients. Note that those restrictions are hardly met in practice.

Almost sure stability of systems of numerical methods with both uniform and non-equidistant step sizes can be studied by the martingale-type approach and Kolmogorov's strong law of large numbers. Methods with balanced terms can help to stabilize numerical dynamics in a.s. long-term sense.

We assembled here only major results w.r.t. the so-called one-point motion of stochastic dynamical systems. There is a similar approach available for the two-point motion process, cf. see ^[54] leading to the concept of $B$-stability. Moment contractivity and moment dissipativity with possible variable, but with nonrandom step sizes is studied in ^[44,47,48]. A more detailed classification of moment-dissipative stochastic numerical methods is given in ^[48] (together with that of underlying SDEs).

In fact, for all uniform step sizes $h > 0$, asymptotically exact numerical methods are available for systems of SDEs. The contributions ^[45,46,54] have shown that midpoint-type methods provide exact ones w.r.t. asymptotic stability. Nowadays, split-step and partial-implicit implementations are claimed to provide the most efficient discretization of both bilinear and nonlinear SDEs (in the combined sense of convergence, stability, and other qualitative properties). Balanced terms as met in balanced implicit methods (1.9) are needed to tackle the problem of constructing and verifying adequate a.s. positivity, a.s. boundedness, and adequate a.s. stability properties.

There are recent attempts to carry out a mean square stability analysis of stochastic numerical approximations in infinite dimension; see Lang et al. ^[25]. Moreover, one can also study the replication of the exact stationary invariant distribution by numerical approximations (as in Burrage and Lythe ^[6] or Weng and Liu ^[62], i.e., asymptotic stability w.r.t. invariant measures for nonlinear SDEs). That type of analysis has been initiated by Talay ^[59] w.r.t. the implicit Euler-Maruyama method. There is the class of split-step methods, which is promising in that direction, too. For example, the split-step technique based on the BEMs (1.2) to approximate invariant distributions of finite-dimensional problems was suggested by Mattingly et al. ^[27]. Additionally, Cui et al ^[7] investigated certain splitting schemes for stochastic versions of the Langevin equation in infinite dimensions. General split-step methods for SDEs are explained in Moro and Schurz ^[31]. All these methods need to still be refined in order to get optimal (i.e., exact) replication of invariant distributions at (possibly all) finite step sizes $h > 0$, but they are justified by convergence as $h \downarrow 0$. We have not focused on the stable replication of invariant probability distributions by numerical methods in infinite dimensions in this review. That would have sprinkled the frame of this concise review paper. However, all in all, (stochastic) "dynamic consistency" is the key in such attempts of adequate analysis (for origin and description of that notion, cf. ^[49,51,53]). Stability investigations in all dimensions play an important role in it, beyond any doubt.

Of course, one can boil down the analysis to commonly suggested one-dimensional ( = 1D) test equations in order to simplify the representation and related studies for elementary classroom experiences. However, the problem of adequateness of those 1D test equations to any dimension $d \ge 1$ remains fairly open. For the very, very restricted subclass of systems of linear SDEs with simultaneously diagonalizable drift and diffusion parts, 1D equations might be considered as relevant and may serve as adequate test equations. However, in the general and more practice-relevant cases, such a simplification is not mathematically justified at all.

All in all, stability investigations enjoy the major tools of Lyapunov functionals, monotonicity relations of positive operators, and semi-martingale decompositions and their convergence theorems, and not just mean square evolutions studied by simple stability functions in 1D. It is possible to construct asymptotically exact numerical methods for SDEs (i.e., in the sense that integration time tends to infinity), while exploiting any nonrandom, equidistant step sizes $h$ (i.e., the stochastic preservation of qualitative properties under unconditional discretization). In that respect, stable midpoint-based and balanced implicit methods (or their combinations with Milstein techniques, cf. ^[9,37]) are the most promising in stochastics. They can be easily implemented in the form of their linear-implicit or split-step modifications while still guaranteeing asymptotic stability for many systems of test equations.

Use of AI tools declaration

The author has not used any Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

There is no conflict of interest of any nature. We claim by no means that our modest review attempt exploits or refers to all available literature.