Milstein Numerical Scheme

• The Milstein scheme is a higher-order discretization method for stochastic differential equations that incorporates derivative corrections to improve simulation accuracy.
• It extends to complex systems including SPDEs, rough SDEs, delay equations, and mean-field models using tailored updates and spectral techniques.
• Adaptations like tamed, truncated, and derivative-free variants ensure strong order convergence even for non-globally Lipschitz or irregular drift conditions.

The Milstein numerical scheme is a higher-order strong discretization method for stochastic differential equations (SDEs) and their extensions, designed to improve simulation accuracy over the basic Euler–Maruyama approach by incorporating derivative information from the diffusion coefficients. The scheme has been extensively developed both for SDEs and for more complex systems such as stochastic partial differential equations (SPDEs), SDEs with rough drivers (e.g., fractional Brownian motion), McKean–Vlasov equations, delay equations, and other mean-field or path-dependent processes. The technical core of Milstein methods is the inclusion of a so-called "Milstein correction," which typically involves terms reflecting the curvature or nonlinear interaction of the driving noise.

1. Core Principles and Mathematical Formulation

The classical Milstein scheme for an Itô SDE

$dX_t = b(X_t)\,dt + \sigma(X_t)\,dB_t$

updates a discretized solution via

$$X_{k+1} = X_k + h\,b(X_k) + \sigma(X_k)\,\Delta B_k + \frac{1}{2}\sigma(X_k)\sigma'(X_k)\left[(\Delta B_k)^2 - h\right],$$

where $\Delta B_k = B_{t_{k+1}} - B_{t_k}$ and $h$ is the time step. This correction term, involving the derivative of $\sigma$, captures interactions between increments of the Brownian motion within each interval, making the scheme strongly convergent of order 1, in contrast to the order $1/2$ achieved by Euler–Maruyama.

The Milstein correction extends to multidimensional and infinite-dimensional settings (e.g., SPDEs), where its derivative structure is generalized through Fréchet or Lions derivatives and, in mean-field models, derivatives with respect to measures. Implementations for rough SDEs (e.g., driven by fBm) and SPDEs invoke spectral projections, commutativity assumptions, or rough path theory to maintain accuracy without prohibitive computational complexity.

2. Extensions to Infinite Dimensions and SPDEs

Milstein-type discretizations for SPDEs rely on a mild formulation (semigroup treatment) and require analogues of the commutativity condition that allows the simplification of double stochastic integrals (Jentzen et al., 2010). For example, the SPDE

$$X_t = e^{A t}\xi + \int_0^t e^{A(t-s)}F(X_s)\,ds + \int_0^t e^{A(t-s)}B(X_s)\,dW_s$$

leads to the Milstein update: $$Y_{n+1}^N = P_N e^{A h}\left[Y_n^N + hF(Y_n^N) + B(Y_n^N)\Delta W_n^{(N)} + \frac{1}{2}B'(Y_n^N)(B(Y_n^N)\Delta W_n^{(N)}) \Delta W_n^{(N)} - \frac{h}{2}\sum_j \eta_j B'(Y_n^N)(B(Y_n^N)g_j)g_j\right],$$ where $P_N$ is the spectral projection and $\Delta W_n^{(N)}$ is the noise increment in the finite-dimensional noise subspace. The commutativity condition, expressed as symmetry in the bilinear form $(v, u) \mapsto B'(v)[B(v)u]$, is essential to avoid intractable cross terms and to achieve computational efficiency and higher-order convergence.

Derivative-free Milstein schemes for SPDEs (Leonhard et al., 2015) replace analytic derivatives of the diffusion coefficient by finite-difference approximations to lower computation costs, achieving an error-versus-cost rate ("effective order of convergence") that surpasses both Euler and classical Milstein, especially as the problem dimension increases.

3. Adaptations for Irregularities and Non-Globally Lipschitz Drift

Several variants of Milstein schemes address SDEs with super-linearly growing or locally bounded coefficients:

• Tamed Milstein schemes (Wang et al., 2011, Fan et al., 2022) regularize the drift by nonlinear damping, e.g.,

$$\widetilde{\mu}(x) = \frac{\mu(x)}{1 + h\|\mu(x)\|},$$

embedded into the standard Milstein update to maintain boundedness and strong order 1 convergence even if only a one-sided Lipschitz condition holds.

• Projected and split-step Milstein methods (Beyn et al., 2015) apply projection or implicit steps to enforce stability and handle superlinear coefficients under global monotonicity.
• Truncated and modified-truncated Milstein schemes (Jiang et al., 2022) truncate coefficients and their derivatives outside a growing bound $h(\Delta)$, only requiring local boundedness rather than global Lipschitzness. With careful selection of truncation functions, strong convergence can approach order 1 and exponential stability of the discretized system is preserved.

These mechanisms ensure robustness in practice for a broader class of SDEs, including neutral or delay equations (Fan et al., 2022, Griggs et al., 21 Aug 2025).

4. Milstein Schemes for SDEs Driven by Rough Signals and Fractional Brownian Motion

For SDEs driven by fractional Brownian motion (fBm) with $H \in (1/4, 1/2)$, Milstein-type techniques are designed to avoid the direct simulation of iterated integrals ("Lévy areas") that are not available in closed form for fBm:

• Modified Milstein schemes (Deya et al., 2010, Huang, 2021) approximate Lévy area terms using products of increments:

$$(1/2) \sum_{i,j} \mathcal{D}^{(i)}\sigma^{(j)}(Y_n) \Delta B_n^{(i)} \Delta B_n^{(j)},$$

rather than the true double integrals. A stochastic backward error analysis justifies that such schemes retain a nearly optimal convergence rate $(2H - 1/2)^-$ for sufficiently smooth coefficients (Huang, 2021). Remarkably, for $H \in (1/4, 1/3]$—where it might be expected that a second-order expansion would fail—these modified Milstein schemes still converge, matching the best achievable rate under practical implementability.

5. Milstein Methods in Mean-Field, Delay, and Path-Dependent Settings

Milstein-type approaches have been successfully generalized to:

• McKean–Vlasov equations and interacting particle systems: Implementations require correction terms related to Lions derivatives to account for dependencies on the law of the process (Bao et al., 2020, Biswas et al., 2022). When the drift is not differentiable, randomization of the drift evaluation over each step (using uniformly drawn random variables) regularizes the scheme and achieves strong order 1 under minimal regularity.
• Stochastic delay-differential equations (SDDEs): Accurate simulation, particularly with indivisible or multiple delays, requires careful treatment of delayed arguments and delayed iterated stochastic integrals (Griggs et al., 21 Aug 2025). Strong order 1 can be preserved using augmented time meshes that contain all necessary evaluation points, as well as refined techniques for simulating double integrals on such meshes.
• Stochastic Volterra equations: For SVIEs with weakly singular kernels and non-differentiable drift, randomized Milstein schemes employing stochastic quadrature replace Taylor expansions and achieve theoretically sharp convergence rates (Wang et al., 2023).

6. Analysis of Invariant Measures and Limit Theorems

Recent work has examined the long-time behavior of Milstein schemes, especially their ability to approximate invariant measures of ergodic SDEs:

• Central limit theorems (CLT) and moderate deviations: For irreducible SDEs with suitable coefficients, both the original system and its Milstein discretization possess unique invariant measures; the empirical measure from the Milstein chain satisfies a CLT with explicit variance determined by solutions to associated Stein equations. Moderate deviation principles have also been established, demonstrating precise exponential approximations for the tail probabilities of the error normalized by estimated variances, under Bernstein-type martingale difference controls (Chen et al., 3 Oct 2025).

7. Implementation, Computational Complexity, and Performance

Computational efficiency is a critical consideration for Milstein-type schemes. In high dimensions, especially for SPDEs, derivative-free variants and cost models based on real function evaluations are essential: the complexity for classical Milstein may scale as $O(N^2 K M)$ (for $N$ spatial modes, $K$ noise modes, $M$ steps), while the derivative-free approaches reduce this to $O(N K M)$, leading to much faster decay of error for a given cost (Leonhard et al., 2015). For delay equations and systems with irregular sampling (e.g., SDDEs with arbitrary delays), mesh augmentation and refined trapezoidal schemes are crucial for retaining high order of convergence without incurring prohibitive interpolation or error accumulation penalties (Griggs et al., 21 Aug 2025).

In summary, the Milstein numerical scheme—across its classical, tamed, truncated, derivative-free, randomized, and high-order extensions—constitutes a central methodology in the strong simulation of SDEs, SPDEs, and related systems. It balances high accuracy (order 1 or more) with practical implementability, scalability, and, when paired with advanced techniques in rough path theory, randomness injection, or mesh design, remains robust even in irregular or non-standard settings.