🔢 mathematics

Implementation of Milstein Schemes for Stochastic Delay-Differential Equations with Arbitrary Fixed Delays

This contribution presents numerical implementation techniques for Euler–Maruyama and Milstein schemes applied to stochastic delay differential equations with arbitrary fixed delays, employing linear interpolation to achieve a convergence order of 1/2 and an extended time grid with varying step sizes to achieve a convergence order of 1, thereby overcoming the limitations of earlier methods restricted to commensurate delays.

Original authors: Mitchell T. Griggs, Kevin Burrage, Pamela M. Burrage

Published 2026-05-05

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Mitchell T. Griggs, Kevin Burrage, Pamela M. Burrage

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to predict the course of a boat drifting down a river. The river is not calm; it is churned up by random waves (this is the "stochastic" part). But here is the twist: the boat's steering does not react only to the waves hitting it right now. It also reacts to waves that hit it five minutes ago or even ten minutes ago.

Scientists call this a Stochastic Delay Differential Equation (SDDE). It is a mathematical problem used to model things like stock prices, biological populations, or chemical reactions where the past influences the present in a random way.

The article you are asking about deals with a very specific, tricky problem: What happens when these "past times" do not align neatly with your clock?

The Problem: The "Misaligned Clock"

To solve these equations on a computer, scientists usually divide time into tiny, equal segments, like a ruler with markings every millimeter. This is called a uniform time grid.

The Simple Case (Divisible Delays): Imagine the delay is exactly 5 seconds and your clock ticks every second. If the computer needs to know what happened 5 seconds ago, it simply looks at the mark exactly 5 ticks back. Easy as ABC.
The Difficult Case (Indivisible Delays): Now imagine the delay is $\pi$ seconds (about 3.14159...), but your clock still ticks every second. If the computer needs to know what happened 3.14 seconds ago, there is no tick mark there! It is stuck between the 3-second mark and the 4-second mark.

For a long time, computers could only handle the "simple case." If the delays were messy numbers (like $\pi$ or $\sqrt{2}$ ) that did not fit the clock, scientists had to either:

Guess: Estimate the missing value by drawing a straight line between the 3-second and 4-second marks (Linear Interpolation). This works, but it is slow and inaccurate.
Give Up: Only study problems where the delays happened to be "nice" numbers.

The Solution: Building a Custom Ruler

The authors, Mitchell Griggs, Kevin Burrage, and Pamela Burrage, say: "Why force the clock to fit the delay? Let us build a custom ruler that fits the delay."

They developed a method to create an extended time grid. Think of it as taking your standard ruler and adding extra, custom-made tick marks exactly where the "past" moments lie, no matter how messy the numbers are.

Here is how their solution works in simple terms:

1. The "Smart" Ruler (Extended Grid)

Instead of forcing the computer to check every second, the computer creates a list of every single moment it needs to check.

It checks the present.
It checks the present minus the delay.
It checks the present minus two delays (because sometimes the past depends on the past of the past).
It sorts all these times into a single, ordered list.

Now, even if the delay is $\pi$ , the computer has a specific tick mark for exactly $\pi$ seconds back. It no longer needs to guess.

2. The "High-Definition" Camera (Milstein Scheme)

The article focuses on a specific mathematical tool called the Milstein Scheme.

The Low-Res Camera (Euler-Maruyama): This is like a blurry photo of the boat. It is fast, but not very accurate. It also works well with the "guessing" method, but it cannot capture the fine details of the boat's movement.
The High-Res Camera (Milstein): This is a super-sharp camera that perfectly captures the boat's movement. However, it requires very precise data. If you feed it a "guessed" value (from the blurry photo method), the image is ruined, and the high-res camera loses its advantage.

The authors' major breakthrough is showing that by using their custom ruler (extended grid), they can feed the High-Res Camera (Milstein Scheme) the exact data it needs, even if the delays are messy numbers. This allows the computer to obtain the "High-Definition" result for problems that were previously too inaccurate to solve precisely.

The Trade-off: More Marks, More Work

There is a catch. Building this custom ruler requires more work.

If you have a simple delay, your ruler has 100 marks.
If you have a messy delay, your ruler might need 1,000 marks to cover all the "in-between" moments.

The article shows that while this requires more computer memory and time, the result is worth it. The "High-Res" method (Order 1) becomes significantly more accurate than the "Low-Res" method (Order 1/2), especially when you need precise answers.

Real Example in the Article

The authors tested this on a model for Financial Options (specifically a "digital call option").

Imagine you are betting on whether a stock price will be above a certain number in the future.
The stock price depends on its volatility (how wildly it jumps) from the past.
If that past time is a messy number (like $\pi$ seconds ago), the old methods would deliver a blurry, inaccurate bet.
The new method builds the custom ruler, calculates the exact past volatility, and delivers a much sharper, more reliable prediction of the bet's value.

Summary

The article is essentially a guide on how to stop forcing messy time delays into clean, uniform time windows. Instead, it teaches us how to build a flexible, custom timeline that contains every single moment the mathematics requires. This allows us to apply our most powerful, high-precision calculation tools to problems that were previously too messy to solve accurately.

Technical Summary: Implementation of Milstein Schemes for Stochastic Delay Differential Equations with Arbitrary Fixed Delays

1. Problem Statement

The paper addresses the numerical simulation of Itô stochastic delay differential equations (SDDEs) with multiple fixed delays ( $\tau_1, \dots, \tau_K$ ) that do not align with a uniform time grid. While theoretical frameworks exist for strong convergence orders of $1/2$ (Euler–Maruyama) and $1$ (Milstein) for such equations, the practical implementation of order-1 schemes is largely restricted to "divisible delays"—that is, cases where all delays are integer multiples of a common step size $h$ .

The core difficulty arises because high-order schemes (particularly the Milstein scheme) require evaluations of the solution at delayed time points ( $t_n - \tau_k$ ) and delayed iterated stochastic integrals ( $I^{\tau_k}_{ij}$ ). When delays are "indivisible" (i.e., no common divisor exists for the delays and the step size), these required time points do not coincide with the simulation grid. Consequently, standard fixed-step implementations fail to evaluate these terms exactly, forcing the use of approximations that degrade the convergence order to $1/2$ .

2. Methodology

The authors propose a systematic approach to achieve a strong convergence order (OoC) of 1 for SDDEs with arbitrary fixed delays. The methodology is divided into two primary strategies:

A. Linear Interpolation (LI) for OoC 1/2

To provide a computationally efficient benchmark, the authors implement schemes with a fixed step size $h$ . When delayed time points ( $t_n - \tau_k$ ) lie between grid points, the delayed solution values are approximated via linear interpolation (LI) between the nearest simulated values.

Schemes: Euler–Maruyama (EM), Milstein, Magnus–Euler–Maruyama (MEM), and Magnus–Milstein (MM) are adapted using LI.
Result: This approach is simple and efficient but limits the convergence order to $1/2$ , even when using refined integral approximations, because the interpolation error dominates the global error.

B. Extended Time Grid for OoC 1

To achieve the theoretical OoC of 1, the authors construct an extended time grid that explicitly includes all time points required for evaluating the Milstein scheme.

Construction: The grid is built by starting with an initial set of observation time points (including the final time $T$ and multiples of the delays) and recursively adding all necessary delayed points ( $t_n - \tau_k$ , $t_n - \tau_l - \tau_k$ , etc.) until the set is closed under delay operations within the interval $[0, T]$ .
Variable Step Size: This results in a non-uniform grid with variable step sizes.
Integral Simulation: The authors extend the trapezoidal approximation method for iterated stochastic integrals (originally developed for divisible delays) to this extended grid. They define an Extended Refined Time Grid (ARTM) by applying the same construction logic with a finer initial step size. This ensures that for every interval $[t_n, t_{n+1}]$ in the scheme grid, the corresponding delayed intervals $[t_n - \tau_k, t_{n+1} - \tau_k]$ are fully resolved by the ARTM, enabling the exact simulation of delayed iterated integrals $I^{\tau_k}_{ij}$ using the trapezoidal rule.

3. Main Contributions

Framework for Extended Grids: The paper introduces a constructive algorithm (Algorithm 1) to generate a time grid that accommodates arbitrary fixed delays, enabling the exact evaluation of delayed scheme values without interpolation.
Extension of Trapezoidal Integration: The authors generalize the trapezoidal approximation for iterated stochastic integrals to the setting of indivisible delays by utilizing the ARTM, thereby ensuring that delayed integration intervals align with grid points.
Convergence Verification: The paper provides a theoretical proof (Theorem 4) and numerical evidence that the Milstein and Magnus–Milstein schemes, when implemented on the extended grid with refined integrals, achieve their theoretical strong convergence order of 1, even for indivisible delays.
Comparative Analysis: The work systematically compares the performance of fixed-step schemes with LI (OoC 1/2) against those on extended grids (OoC 1), demonstrating that while the extended grid incurs higher computational costs (due to grid growth), it is necessary to achieve order-1 accuracy.

4. Results

Convergence Orders: Numerical experiments confirm that both fixed-step schemes and those on extended grids behave as expected for divisible delays. However, for indivisible delays, fixed-step schemes using linear interpolation degrade to OoC 1/2. In contrast, schemes on extended grids maintain OoC 1.
Grid Size vs. Delay Magnitude: The size of the extended grid grows quadratically with the number of delay multiples. The authors note that the method is most practical when delays are large relative to the final time $T$ . For small delays, the grid size can become prohibitively large.
Computational Efficiency: Table 3 shows that while the extended grid requires more computation time than fixed-step methods with LI, it can be more efficient than LI methods when high accuracy is required. For instance, the Milstein scheme on the extended grid proved faster than the EM-LI scheme in specific configurations with indivisible delays to reach a certain error level.
Application to Financial Models: The method was applied to a Delayed Local Volatility Model (DLVM) for pricing digital call options. The Milstein scheme on the extended grid provided stable and accurate estimates of option values at coarser time resolutions compared to lower-order methods.

5. Significance and Claims

The paper claims to bridge a critical gap between the theoretical existence of high-order SDDE solvers and their practical implementation for general delays.

Practical Implementation: The authors emphasize that while the continuous Milstein scheme was theoretically established by Kloeden and Shardlow [14], no standard implementation strategy existed for indivisible delays. This work provides that strategy.
Robustness: The method proves robust across various delay configurations, including cases with multiple delays and irrational delay values (e.g., $\pi/4$ ).
Limitations and Future Work: The authors modestly acknowledge that the size of the extended grid can become a bottleneck for small delays. They suggest that for such cases, the LI method (OoC 1/2) remains a viable alternative. Mentioned future extensions include handling time-dependent or distributed delays, as well as exploring adaptive refinement or parallel implementations to mitigate computational costs.

In summary, the paper concludes that by abandoning the restriction to a uniform time grid in favor of a systematically constructed extended grid, strong convergence of order 1 for SDDEs with arbitrary fixed delays becomes achievable, enabling more precise simulations for complex stochastic models in fields such as finance and ecology.