Stochastic Theta Method Analysis

• The stochastic theta method is a one-parameter family of discretizations for SDEs, balancing implicit and explicit drift approaches.
• It accommodates diverse drift and diffusion conditions, ensuring stability and strong convergence in both classical and nonstandard noise settings.
• The method efficiently approximates long-term statistical properties, including stationary distributions and invariant measures, vital for asymptotic analysis.

Searching arXiv for recent and foundational papers on the stochastic theta method to ground the article in the supplied literature. {"query":"stochastic theta method SDE arXiv stochastic theta method invariant measure superlinear time-changed Lévy SDAE", "max_results": 10} The stochastic theta method is a one-parameter family of implicit-explicit one-step discretizations for stochastic differential equations. In the classical Itô setting

$dx(t)=f(x(t))\,dt+g(x(t))\,dB(t), \qquad x(0)=x_0,$

it is written as

$$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$

where $t_k=kh$, $\Delta B_k=B(t_{k+1})-B(t_k)$, and $\theta\in[0,1]$ controls how much of the drift is treated implicitly (Jiang et al., 2018). The method interpolates between explicit Euler–Maruyama at $\theta=0$, backward Euler–Maruyama at $\theta=1$, and semi-implicit schemes for $0<\theta<1$; recent analyses treat stationary distributions, invariant measures, mean-square contractivity, random periodic solutions, time-changed Brownian and Lévy noise, free stochastic differential equations, stochastic differential-algebraic equations, neutral stochastic delay equations, and stochastic hybrid systems (Jiang et al., 2018, Li et al., 18 Jun 2026, Niu et al., 2024, Chen et al., 28 Mar 2026, Chen et al., 15 Sep 2025).

1. Classical formulation and parameterization

The defining feature of the stochastic theta method is the split treatment of the drift. The parameter $\theta\in[0,1]$ controls the balance between an implicit term $\theta h\,f(X_{k+1})$ and an explicit term $$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$0, while the diffusion is treated explicitly as $$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$1 (Jiang et al., 2018). The same structural pattern recurs in later variants for time-changed equations, free stochastic differential equations, and stochastic differential-algebraic equations, even when the state space, filtration, or noise model changes (Chen et al., 27 Mar 2025, Niu et al., 2024, Chen et al., 15 Sep 2025).

For the classical Itô case, the step can be rewritten as

$$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$2

which isolates the implicitness in an inverse map (Jiang et al., 2018). In the time-changed setting, an analogous device is used with

$$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$3

so that the update is defined by solving $$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$4 on the operational-time grid induced by the subordinator (Chen et al., 27 Mar 2025). In free probability, the same pattern appears as

$$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$5

and $$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$6 is identified as the free backward Euler method (Niu et al., 2024).

The method is therefore best viewed as a one-parameter family rather than a single algorithm. This family contains explicit, semi-implicit, and fully implicit drift treatments, and the available theory shows that the choice of $$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$7 materially changes solvability requirements, stability thresholds, and the coefficient classes that can be handled (Jiang et al., 2018, D'Ambrosio et al., 2020).

2. Assumptions, solvability, and admissible coefficient growth

A central theme in the literature is that the assumptions required on $$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$8 and $$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$9 depend on $t_k=kh$0. For $t_k=kh$1, the classical stationary-distribution analysis requires stronger global Lipschitz-type conditions on both drift and diffusion. For $t_k=kh$2, the method can tolerate super-linear drift, requiring only a Lipschitz condition on the diffusion together with a one-sided dissipativity condition on the drift (Jiang et al., 2018).

In the classical monotonicity framework, well-posedness follows from

$t_k=kh$3

together with the step-size condition

$t_k=kh$4

Under this condition, $t_k=kh$5 is monotone and invertible, so $t_k=kh$6 is uniquely solvable (Jiang et al., 2018). For $t_k=kh$7, the same paper replaces global Lipschitz drift assumptions by

$t_k=kh$8

$t_k=kh$9

with

$\Delta B_k=B(t_{k+1})-B(t_k)$0

thereby shifting the theory from global Lipschitz continuity to dissipativity and linear-growth control (Jiang et al., 2018).

A more recent invariant-measure analysis allows both drift and diffusion to grow super-linearly. There the assumptions are

$\Delta B_k=B(t_{k+1})-B(t_k)$1

$\Delta B_k=B(t_{k+1})-B(t_k)$2

and

$\Delta B_k=B(t_{k+1})-B(t_k)$3

with $\Delta B_k=B(t_{k+1})-B(t_k)$4 (Li et al., 18 Jun 2026). In this regime, the method is well defined for any $\Delta B_k=B(t_{k+1})-B(t_k)$5 under the monotonicity assumption, and the analysis yields uniform second-moment bounds and exponential contraction between two numerical solutions (Li et al., 18 Jun 2026).

The same pattern persists in nonstandard settings. For time-changed stochastic differential equations with local Lipschitz coefficients, the implicit equation $\Delta B_k=B(t_{k+1})-B(t_k)$6 is uniquely solvable for sufficiently small $\Delta B_k=B(t_{k+1})-B(t_k)$7, and the threshold

$\Delta B_k=B(t_{k+1})-B(t_k)$8

appears in the moment analysis for $\Delta B_k=B(t_{k+1})-B(t_k)$9 (Chen et al., 27 Mar 2025). For index-1 stochastic differential-algebraic equations with time-dependent singular matrices, well-posedness of the stochastic theta step requires

$\theta\in[0,1]$0

and the discrete solution remains on the constraint manifold at all time levels (Chen et al., 15 Sep 2025, Zhu et al., 4 Jun 2026).

3. Stationary distributions, invariant measures, and random periodic solutions

One major line of work studies the long-time statistical objects generated by the numerical chain rather than only finite-time pathwise error. In the nonlinear Itô setting, the stochastic theta method generates a Markov chain with a unique stationary distribution $\theta\in[0,1]$1 once tightness, asymptotic contractivity of two trajectories, and uniform finite-time boundedness are verified. The numerical stationary distribution converges weakly to the true stationary distribution $\theta\in[0,1]$2 in the bounded-Lipschitz metric: $\theta\in[0,1]$3 (Jiang et al., 2018).

For super-linearly growing drift and diffusion, the invariant-measure theory is formulated in $\theta\in[0,1]$4-Wasserstein distance. Under Assumptions 2.1, 2.2, and 2.3, and for

$\theta\in[0,1]$5

the numerical method converges in $\theta\in[0,1]$6-Wasserstein distance to a unique invariant measure $\theta\in[0,1]$7, and the approximation error satisfies

$\theta\in[0,1]$8

for any $\theta\in[0,1]$9 (Li et al., 18 Jun 2026).

A different long-time object appears in periodic nonautonomous equations. For semilinear stochastic differential equations with periodic coefficients and non-globally Lipschitz nonlinearities, stochastic theta methods with $\theta=0$0 admit a unique random periodic solution $\theta=0$1, and this numerical random periodic solution converges strongly in mean square to the exact random periodic solution $\theta=0$2. The convergence order is $\theta=0$3 for multiplicative noise and $\theta=0$4 for additive noise under the additional smoothness assumption (Chen et al., 2024).

Long-time object	Representative conclusion	Source
Stationary distribution	unique $\theta=0$5, with $\theta=0$6	(Jiang et al., 2018)
Invariant measure	unique $\theta=0$7, with $\theta=0$8	(Li et al., 18 Jun 2026)
Random periodic solution	unique $\theta=0$9, with mean-square order $\theta=1$0 or $\theta=1$1	(Chen et al., 2024)

These results show that the stochastic theta method is used not only for trajectory approximation but also for approximating asymptotic probabilistic structures. The literature distinguishes clearly between stationary distributions for time-homogeneous Markov chains, invariant measures in Wasserstein distance, and random periodic solutions in periodic nonautonomous systems (Jiang et al., 2018, Li et al., 18 Jun 2026, Chen et al., 2024).

4. Strong and weak convergence across problem classes

For classical and generalized Itô equations, the dominant strong order is often $\theta=1$2, but the precise statement depends on the noise model and coefficient regularity. In free stochastic differential equations, the free stochastic theta methods strongly converge in $\theta=1$3 with order $\theta=1$4: $\theta=1$5 (Niu et al., 2024). For index-1 stochastic differential-algebraic equations under global monotonicity conditions and non-global Lipschitz coefficients, each stochastic theta method with $\theta=1$6 achieves the mean square convergence rate

$\theta=1$7

that is, RMS strong order $\theta=1$8 (Chen et al., 15 Sep 2025).

In time-changed settings, the convergence rate may depend explicitly on the random clock. For stochastic differential equations driven by a time-changed Brownian motion and locally Lipschitz time-space-dependent coefficients, the main theorem states that for $\theta=1$9, $0<\theta<1$0, and any $0<\theta<1$1,

$0<\theta<1$2

The final theorem is a convergence-to-zero result rather than an explicit optimal rate statement, but the derivation indicates that the error is controlled by terms of the form $0<\theta<1$3 and $0<\theta<1$4, together with stochastic-clock-dependent factors (Chen et al., 27 Mar 2025).

For time-changed Lévy noise beyond Lipschitz continuity, the strong convergence rate is stated explicitly: $0<\theta<1$5 under Assumptions 2.1–2.5, $0<\theta<1$6, $0<\theta<1$7, and $0<\theta<1$8 (Chen, 18 Aug 2025). The rate is the minimum of the time regularities of the drift, diffusion, and jump coefficients and the intrinsic $0<\theta<1$9 penalty induced by the inverse $\theta\in[0,1]$0-stable subordinator (Chen, 18 Aug 2025).

Weak convergence results are also available. For stochastic differential equations driven by time-changed Lévy noise under global Lipschitz and linear-growth conditions, the stochastic theta method with $\theta\in[0,1]$1 achieves weak order one: $\theta\in[0,1]$2 for $\theta\in[0,1]$3 (Chen et al., 28 Mar 2026). The proof combines a global weak convergence estimate of order one for the corresponding non-time-changed Lévy equation, the Kolmogorov backward partial integro-differential equation, an approximation of the inverse subordinator, and the duality principle (Chen et al., 28 Mar 2026). A parallel weak-order-one theory has been established for structure-preserving stochastic theta methods for index-1 stochastic differential-algebraic equations with time-dependent singular matrices: $\theta\in[0,1]$4 for $\theta\in[0,1]$5 and $\theta\in[0,1]$6 (Zhu et al., 4 Jun 2026).

5. Stability, contractivity, and the role of implicitness

The stability theory of stochastic theta methods is highly sensitive to $\theta\in[0,1]$7. For classical Itô stochastic differential equations satisfying a one-sided Lipschitz condition for the drift and a global Lipschitz condition for the diffusion, the exact dynamics are exponentially mean-square contractive when

$\theta\in[0,1]$8

For the stochastic $\theta\in[0,1]$9-Maruyama method,

$\theta h\,f(X_{k+1})$0

where

$\theta h\,f(X_{k+1})$1

The admissible contractivity region is

$\theta h\,f(X_{k+1})$2

so the implicit Euler–Maruyama method is unconditionally mean-square contractive (D'Ambrosio et al., 2020).

The advantage of larger $\theta h\,f(X_{k+1})$3 is even more visible in nonclassical noise models. For free stochastic differential equations, if $\theta h\,f(X_{k+1})$4, then the free stochastic theta method is exponentially stable in mean square for

$\theta h\,f(X_{k+1})$5

whereas for $\theta h\,f(X_{k+1})$6 the free backward Euler method is exponentially stable in mean square for any $\theta h\,f(X_{k+1})$7 (Niu et al., 2024). For time-changed stochastic differential equations driven by Brownian motion, the main numerical stability theorem states that the stochastic theta method is asymptotically mean-square stable for all $\theta h\,f(X_{k+1})$8 whenever $\theta h\,f(X_{k+1})$9; for $$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$00, stability requires

$$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$01

under the additional bound $$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$02 (Chen et al., 27 Mar 2025).

For stochastic differential equations driven by fractional Brownian motion with Hurst parameter $$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$03, the mean-square stability picture is more delicate because the increments are correlated and there is no martingale structure. For the linear test problem, the stochastic theta method preserves mean-square stability if

$$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$04

and also if

$$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$05

By contrast, if $$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$06, the method does not preserve stability unconditionally (Li et al., 2021). The same paper identifies an open regime,

$$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$07

for which numerical evidence suggests stability may still hold (Li et al., 2021).

Across these analyses, a common conclusion is explicit in the source material: larger $$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$08 means more implicitness, and more implicitness improves stability (Chen et al., 27 Mar 2025, Niu et al., 2024, Li et al., 2021).

6. Extensions, related schemes, and terminological ambiguities

The stochastic theta method has been extended well beyond standard Brownian-motion-driven Itô equations. For neutral stochastic differential delay equations with non-globally Lipschitz coefficients, the theta-EM scheme

$$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$09

has strong order $$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$10 in the Brownian case,

$$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$11

and almost sure order $$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$12; a corresponding jump version also converges strongly, with a weaker rate in the sense described in the paper (Tan et al., 2017).

For stochastic hybrid systems written as Random Time Change equations, the $$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$13-Maruyama family is analyzed as a semi-implicit Maruyama-type one-step method: $$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$14 In this setting, the principal theorem states that strong consistency implies strong convergence in mean, and the local error analysis is carried out via Itô–Taylor expansions of the exact solution and the approximation process (Riedler et al., 2013).

A common misconception concerns the appearance of the word “theta” in model names rather than in numerical discretizations. The paper on a “mean-reverting theta-rho model” for asset prices does not introduce or analyze a stochastic theta method. It studies a highly nonlinear mean-reverting asset price model with stochastic volatility and develops a truncated Euler–Maruyama method for the coupled system; the paper itself states that, although earlier Euler–Maruyama work on theta-type models is discussed in the literature review, the new method is explicitly a truncated EM method rather than a stochastic theta method (Coffie, 2022).

This terminological distinction matters because “theta” may refer either to a model parameter, as in mean-reverting $$X_{k+1}=X_k+\theta h\,f(X_{k+1})+(1-\theta)h\,f(X_k)+g(X_k)\Delta B_k,$$15-type processes, or to the implicitness parameter in the stochastic theta discretization. The modern numerical-analysis literature uses “stochastic theta method” in the latter sense: a one-parameter implicit-explicit scheme whose drift splitting determines solvability, stability, and approximation properties across a broad range of stochastic dynamical systems (Jiang et al., 2018, D'Ambrosio et al., 2020, Li et al., 18 Jun 2026).