pith. sign in

arxiv: 2508.15365 · v2 · submitted 2025-08-21 · 🧮 math.NA · cs.NA

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

Mitchell T. Griggs , Kevin Burrage , Pamela M. Burrage This is my paper

Pith reviewed 2026-05-18 22:20 UTC · model grok-4.3

classification 🧮 math.NA cs.NA
keywords stochastic delay-differential equationsMilstein schemeEuler-Maruyama schemenumerical simulationconvergence orderaugmented time meshlinear interpolationiterated stochastic integrals
0
0 comments X

The pith

Numerical schemes achieve strong orders 1/2 and 1 for SDDEs with arbitrary fixed delays

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

This paper develops implementation techniques for solving stochastic delay-differential equations when fixed delays do not coincide with the simulation time mesh. For schemes with strong order 1/2, it uses a fixed step size and linear interpolation to approximate values at delayed times. For strong order 1 schemes like Milstein, it builds an augmented mesh including all necessary points and extends a method to simulate the required delayed stochastic integrals. These approaches allow general delays without the previous restriction to divisible cases. Computational examples verify that the expected convergence orders are attained.

Core claim

The central claim is that simulation techniques using linear interpolation on fixed meshes and augmented time meshes with extended delayed integral simulation allow Euler-Maruyama and Milstein schemes to be applied to SDDEs with multiple fixed delays that are not restricted to the time mesh, achieving their theoretical strong convergence orders of 1/2 and 1 respectively, as confirmed by numerical experiments.

What carries the argument

The augmented time mesh construction that incorporates all delayed time points together with the extension of the established method for simulating delayed iterated stochastic integrals on that mesh.

Load-bearing premise

Linear interpolation on the fixed mesh and the augmented mesh with extended integral simulation do not degrade the strong convergence orders below their theoretical values.

What would settle it

Numerical tests with specific indivisible delays showing measured convergence rates matching 0.5 for Euler-Maruyama and 1.0 for Milstein would support the claim, while significant deviation would falsify it.

Figures

Figures reproduced from arXiv: 2508.15365 by Kevin Burrage, Mitchell T. Griggs, Pamela M. Burrage.

Figure 1
Figure 1. Figure 1: Integration regions of Iij (tn, tn+1) and I τk ij (tn, tn+1) shown together on the same plane. To achieve strong order of convergence 1, the integrals are simulated using approximations (18) based on the values of the Wiener paths evaluated on the finest time mesh. If the finest mesh uses a fixed step size hR that divides τk then the delayed and region also aligns with this mesh. In the case i = j, simulat… view at source ↗
Figure 2
Figure 2. Figure 2: Mean-square error graphs for the numerical solutions of equation ( [PITH_FULL_IMAGE:figures/full_fig_p009_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: Visual depiction of the (past and present) scheme values required to simulate a subsequent value, using the Milstein [PITH_FULL_IMAGE:figures/full_fig_p011_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: Integration regions of Iij (tn, tn+1) and I τk ij (tn, tn+1), shown over a refined time mesh of fixed step size h R, where the scheme step size h = Fnh R (for natural number Fn) does not divide τk. In this case, the approximation (18) is not suited for I τk ij (tn, tn+1), since the refined mesh does not coincide with the corners of the delayed integration region. 3.4. Discussion of solutions for these issu… view at source ↗
Figure 5
Figure 5. Figure 5: Simulating the delayed scheme value Y (tn − τk) using linear interpolation between the scheme values at the nearest mesh points. The nearest mesh points are those immediately before and after the delayed time tn − τk. Definition 5. Suppose Y = (Yn)n=−p,...,N is a numerical approximation of the solution X = (X(tn))t∈[−τ,T] for (1), with Yn = ϕ(tn) for tn ⩽ 0. Further to this, suppose Y is a scheme of fixed … view at source ↗
Figure 6
Figure 6. Figure 6: Mean-square error graphs for the fixed-step LI numerical solutions of equation ( [PITH_FULL_IMAGE:figures/full_fig_p014_6.png] view at source ↗
Figure 7
Figure 7. Figure 7: Augmented mesh constructions used to simulate numerical schemes with OoC [PITH_FULL_IMAGE:figures/full_fig_p017_7.png] view at source ↗
Figure 8
Figure 8. Figure 8: Mean-square error graphs for the numerical solutions of equation ( [PITH_FULL_IMAGE:figures/full_fig_p018_8.png] view at source ↗
Figure 9
Figure 9. Figure 9: Mean-square error graphs for the numerical solutions of equation ( [PITH_FULL_IMAGE:figures/full_fig_p019_9.png] view at source ↗
read the original abstract

This paper develops methods for numerically solving stochastic delay-differential equations (SDDEs) with multiple fixed delays that do not align with a uniform time mesh. We focus on numerical schemes of strong convergence orders $1/2$ and $1$, such as the Euler--Maruyama and Milstein schemes, respectively. Although numerical schemes for SDDEs with delays $\tau_1,\ldots,\tau_K$ are theoretically established, their implementations require evaluations at both present times such as $t_n$, and also at delayed times such as $t_n-\tau_k$ and $t_n-\tau_l-\tau_k$. As a result, previous simulations of these schemes have been largely restricted to the case of divisible delays. We develop simulation techniques for the general case of indivisible delays where delayed times such as $t_n-\tau_k$ are not restricted to a uniform time mesh. To achieve order of convergence (OoC) $1/2$, we implement the schemes with a fixed step size while using linear interpolation to approximate delayed scheme values. To achieve OoC $1$, we construct an augmented time mesh that includes all time points required to evaluate the schemes, which necessitates using a varying step size. We also introduce a technique to simulate delayed iterated stochastic integrals on the augmented time mesh, by extending an established method from the divisible-delays setting. We then confirm that the numerical schemes achieve their theoretical convergence orders with computational examples.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

2 major / 1 minor

Summary. The manuscript develops practical implementation techniques for Euler--Maruyama (strong order 1/2) and Milstein (strong order 1) schemes applied to SDDEs with multiple fixed delays that are not necessarily divisible by the time step. For order 1/2 it employs a fixed-step discretization together with linear interpolation to evaluate the scheme at delayed times. For order 1 it constructs an augmented mesh that includes all required delay points (necessitating variable step sizes) and extends an existing method to simulate the additional delayed iterated stochastic integrals. The authors state that computational examples confirm that both implementations attain their respective theoretical strong convergence orders.

Significance. If the proposed implementations preserve the known theoretical orders without introducing order-reducing interpolation or mesh-variation errors, the work would remove a long-standing practical restriction on the use of Milstein-type schemes for SDDEs, thereby enabling higher-accuracy simulation for the general case of incommensurate delays that arise in many applications. The extension of the delayed-integral simulation technique to the augmented-mesh setting is a concrete technical contribution that builds directly on prior divisible-delay results.

major comments (2)
  1. [Abstract / OoC 1/2 implementation] Abstract, paragraph on OoC 1/2 implementation: the assertion that linear interpolation on a fixed mesh preserves strong order 1/2 is load-bearing for the central claim, yet the manuscript supplies no explicit bound on the interpolation remainder inside the Itô-Taylor expansion for the delayed terms, nor does it quantify how this remainder interacts with the stochastic integrals.
  2. [Abstract / OoC 1 implementation] Abstract, paragraph on OoC 1 implementation: the augmented-mesh construction with variable step sizes and the extended simulation of delayed iterated integrals are presented as achieving strong order 1, but the text does not verify that the variable steps satisfy the Lipschitz or moment conditions under which the Milstein scheme retains order 1 when delays are indivisible; this is a load-bearing gap because the schemes were previously established only for the divisible-delay case.
minor comments (1)
  1. The manuscript should include explicit error tables, comparison against known exact solutions (or high-accuracy reference solutions), and a brief discussion of how the observed convergence rates behave under different noise structures or delay ratios.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for the careful reading and for identifying the key theoretical points that require clarification in our presentation. We address each major comment below and outline the revisions we will make to strengthen the manuscript.

read point-by-point responses

  1. Referee: [Abstract / OoC 1/2 implementation] Abstract, paragraph on OoC 1/2 implementation: the assertion that linear interpolation on a fixed mesh preserves strong order 1/2 is load-bearing for the central claim, yet the manuscript supplies no explicit bound on the interpolation remainder inside the Itô-Taylor expansion for the delayed terms, nor does it quantify how this remainder interacts with the stochastic integrals.

    Authors: We agree that an explicit treatment of the interpolation remainder would improve the rigor of the argument. In the revised manuscript we will add a short paragraph (or appendix remark) deriving a mean-square bound on the linear interpolation error for the delayed arguments. Under the standard Lipschitz and linear-growth assumptions, this error is O(h) uniformly in the mean-square sense and, when inserted into the Itô-Taylor expansion of the delayed terms, produces a local truncation contribution that remains compatible with the global strong order 1/2 of the Euler–Maruyama scheme. The non-anticipative nature of the interpolation error allows it to be absorbed into the existing Gronwall-type estimates used for SDDEs; we will cite the relevant moment bounds to make this interaction explicit. revision: yes

  2. Referee: [Abstract / OoC 1 implementation] Abstract, paragraph on OoC 1 implementation: the augmented-mesh construction with variable step sizes and the extended simulation of delayed iterated integrals are presented as achieving strong order 1, but the text does not verify that the variable steps satisfy the Lipschitz or moment conditions under which the Milstein scheme retains order 1 when delays are indivisible; this is a load-bearing gap because the schemes were previously established only for the divisible-delay case.

    Authors: The referee is correct that the original strong-order-1 analysis for SDDEs assumed commensurate delays permitting a uniform mesh. For the augmented mesh the step sizes are generated by the fixed delays together with a base step h; consequently there are only finitely many distinct step lengths, each bounded above and below by positive multiples of h. Under the global Lipschitz and linear-growth conditions already stated in the paper, these variable steps satisfy the uniform moment bounds and local Lipschitz requirements needed for the Milstein local truncation error to remain O(h^{3/2}) in the mean-square sense. The extension of the delayed iterated-integral simulation technique preserves the necessary adaptedness and moment properties. We will insert a concise verification subsection that records these facts and thereby justifies retention of strong order 1. revision: yes

Circularity Check

0 steps flagged

No circularity: implementation extends external theory with numerical verification

full rationale

The paper appeals to pre-existing theoretical results establishing strong convergence orders for Milstein-type schemes on SDDEs (with divisible delays) and supplies concrete implementation techniques—fixed-step linear interpolation for order 1/2 and augmented-mesh construction plus extended delayed-integral simulation for order 1—together with computational confirmation on examples. No step reduces a claimed result to a fitted parameter, self-defined quantity, or load-bearing self-citation chain inside the work; the central claims rest on external theory plus direct numerical testing rather than internal tautology.

Axiom & Free-Parameter Ledger

0 free parameters · 1 axioms · 0 invented entities

The work rests on standard stochastic calculus for SDDEs and on the existence of theoretical convergence proofs for the schemes when delays are commensurate. No new free parameters, axioms, or invented entities are introduced in the abstract.

axioms (1)
  • domain assumption Theoretical strong convergence orders of Euler-Maruyama (1/2) and Milstein (1) schemes hold for SDDEs when all required delayed values lie on the mesh.
    Invoked when the authors state they achieve the theoretical orders after implementing the general-delay case.

pith-pipeline@v0.9.0 · 5798 in / 1471 out tokens · 32592 ms · 2026-05-18T22:20:30.640836+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

32 extracted references · 32 canonical work pages

  1. [1]

    Lee, J.-H

    M.-K. Lee, J.-H. Kim, J. Kim, A delay financial model with stochastic volatility; martingale method, Physica A: Statistical Mechanics and its Applications 390 (2011) 2909–2919.doi:10.1016/j.physa.2011.03.032

    work page doi:10.1016/j.physa.2011.03.032 2011

  2. [2]

    Hodgson, C

    R. Hodgson, C. Reisinger, Numerical methods for mckean–vlasov equations via the mean-field and master equations approach, Applied Numerical Mathematics 179 (2022) 1640–1684.doi:10.1214/18-AAP1429

    work page doi:10.1214/18-aap1429 2022

  3. [3]

    Boulet, R

    J. Boulet, R. Balasubramaniam, A. Daffertshofer, A. Longtin, Stochastic two-delay differential model of delayed visual feedback effects on postural dynamics, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 368 (2010) 423–438.doi:10.1098/rsta.2009.0214

    work page doi:10.1098/rsta.2009.0214 2010

  4. [4]

    Hartung, J

    F. Hartung, J. Turi, Parameter identification in a respiratory control system model with delay, in: J. Batzel, M. Bachar, F. Kappel (Eds.), Mathematical Modeling and Validation in Physiology, Vol. 2064 of Lecture Notes in Mathematics, Springer, Berlin, Germany, 2013, pp. 105–118.doi:10.1007/978-3-642-32882-4

    work page doi:10.1007/978-3-642-32882-4 2064

  5. [5]

    W. Wang, L. Wang, W. Chen, Stochastic nicholson’s blowflies delayed differential equations, Applied Mathematics Letters 87 (2019) 20–26.doi:10.1016/j.aml.2018.07.020

    work page doi:10.1016/j.aml.2018.07.020 2019

  6. [6]

    Maruyama, Continuous markov processes and stochastic equations, Rendiconti del Circolo Matematico di Palermo 4 (1955) 48–90.doi:10.1007/BF02846028

    G. Maruyama, Continuous markov processes and stochastic equations, Rendiconti del Circolo Matematico di Palermo 4 (1955) 48–90.doi:10.1007/BF02846028

    work page doi:10.1007/bf02846028 1955

  7. [7]

    Milstein, Approximate integration of stochastic differential equations, Theory of Probability & Its Applications 19 (1975) 557–562.doi:10.1137/1119062

    G. Milstein, Approximate integration of stochastic differential equations, Theory of Probability & Its Applications 19 (1975) 557–562.doi:10.1137/1119062

    work page doi:10.1137/1119062 1975

  8. [8]

    Küchler, E

    U. Küchler, E. Platen, Strong discrete time approximation of stochastic differential equations with time delay, Math- ematics and Computers in Simulation 54 (2000) 189–205.doi:10.1016/S0378-4754(00)00224-X

    work page doi:10.1016/s0378-4754(00)00224-x 2000

  9. [9]

    Hofmann, T

    N. Hofmann, T. Müller-Gronbach, A modified Milstein scheme for approximation of stochastic delay differential equations with constant time lag, Journal of Computational and Applied Mathematics 197 (2006) 89–121.doi: 10.1016/j.cam.2005.10.027

    work page doi:10.1016/j.cam.2005.10.027 2006

  10. [10]

    Milstein, Numerical Integration of Stochastic Differential Equations, Springer, Berlin, Germany, 1995

    G. Milstein, Numerical Integration of Stochastic Differential Equations, Springer, Berlin, Germany, 1995. doi: 10.1007/978-94-015-8455-5

    work page doi:10.1007/978-94-015-8455-5 1995

  11. [11]

    W. Cao, Z. Zhang, G. Karniadakis, Numerical methods for stochastic delay differential equations via the wong–zakai approximation, SIAM Journal on Scientific Computing 37 (1) (2015) A295–A318.doi:10.1137/130942024

    work page doi:10.1137/130942024 2015

  12. [12]

    and Manton, J

    P. Kloeden, T. Shardlow, The Milstein scheme for stochastic delay differential equations without using anticipative calculus, Stochastic Analysis and Applications 30 (2) (2012) 181–202.doi:10.1080/07362994.2012.628907

    work page doi:10.1080/07362994.2012.628907 2012

  13. [13]

    Griggs, K

    M. Griggs, K. Burrage, P. Burrage, Magnus methods for stochastic delay-differential equations, arXiv preprint (2025). doi:10.48550/arXiv.2506.16908

    work page doi:10.48550/arxiv.2506.16908 2025

  14. [14]

    Bellman, K

    R. Bellman, K. Cooke, Differential-Difference Equations, RAND Corporation, Santa Monica, CA, USA, 1963

    work page 1963

  15. [15]

    Kloeden, E

    P. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, 3rd Edition, Springer, Berlin, Ger- many, 1999. doi:10.1007/978-3-662-12616-5

    work page doi:10.1007/978-3-662-12616-5 1999

  16. [16]

    Magnus, On the exponential solution of differential equations for a linear operator, Communications on Pure and Applied Mathematics 7 (1954) 649–673.doi:10.1002/cpa.3160070404

    W. Magnus, On the exponential solution of differential equations for a linear operator, Communications on Pure and Applied Mathematics 7 (1954) 649–673.doi:10.1002/cpa.3160070404

    work page doi:10.1002/cpa.3160070404 1954

  17. [17]

    Burrage, P

    K. Burrage, P. Burrage, High strong order methods for non-commutative stochastic ordinary differential equation systems and the Magnus formula, Physica D: Nonlinear Phenomena 133 (1999) 34–48.doi:10.1016/S0167-2789(99) 00097-4

    work page doi:10.1016/s0167-2789(99 1999

  18. [18]

    K. Kamm, S. Pagliarani, A. Pascucci, On the stochastic Magnus expansion and its application to spdes, Journal of Scientific Computing 89 (2021) 1–34.doi:10.1007/s10915-021-01633-6

    work page doi:10.1007/s10915-021-01633-6 2021

  19. [19]

    G. Yang, K. Burrage, Y. Komori, P. Burrage, X. Ding, A class of new Magnus-type methods for semi-linear non-commutative Itô stochastic differential equations, Numerical Algorithms 88 (2021) 1641–1665.doi:10.1007/ s11075-021-01089-7

    work page 2021

  20. [20]

    Q. Guo, X. Mao, R. Yue, The truncated Euler–Maruyama method for stochastic differential delay equations, Numer- ical Algorithms 255 (2018) 599–624.doi:10.1007/s11075-017-0391-0. 21

    work page doi:10.1007/s11075-017-0391-0 2018

  21. [21]

    Milstein, M

    G. Milstein, M. Tretyakov, Stochastic Numerics for Mathematical Physics, 2nd Edition, Springer, Berlin, Germany,

    work page

  22. [22]

    doi:10.1007/978-3-030-82040-4

    work page doi:10.1007/978-3-030-82040-4

  23. [23]

    Kloeden, E

    P. Kloeden, E. Platen, I. Wright, The approximation of multiple stochastic integrals, Stochastic Analysis and Appli- cations 40 (4) (1992) 431–441.doi:10.1080/07362999208809281

    work page doi:10.1080/07362999208809281 1992

  24. [24]

    M. Wiktorsson, Joint characteristic function and simultaneous simulation of iterated Itô integrals for multiple inde- pendent Brownian motions, Annals of Applied Probability 11 (2) (2001) 470–487.doi:10.1214/aoap/1015345301

    work page doi:10.1214/aoap/1015345301 2001

  25. [25]

    Mrongowius, A

    J. Mrongowius, A. Rößler, On the approximation and simulation of iterated stochastic integrals and the corresponding lévy areas in terms of a multidimensional brownian motion, Stochastic Analysis and Applications 40 (1) (2022) 397–

    work page 2022

  26. [26]

    doi:10.1080/07362994.2021.1922291

    work page doi:10.1080/07362994.2021.1922291 2021

  27. [27]

    D. Kuznetsov, Development and application of the Fourier method for the numerical solution of Itô stochas- tic differential equations, Computational Mathematics and Mathematical Physics 58 (2018) 1058–1070. doi: 10.1134/S0965542518070096

    work page doi:10.1134/s0965542518070096 2018

  28. [28]

    D. Kuznetsov, A comparative analysis of efficiency of using the Legendre polynomials and trigonometric functions for the numerical solution of Itô stochastic differential equations, Computational Mathematics and Mathematical Physics 59 (2019) 1236–1250.doi:10.1134/S0965542519080116

    work page doi:10.1134/s0965542519080116 2019

  29. [29]

    Hu, S.-E

    Y. Hu, S.-E. Mohammed, F. Yan, Discrete-time approximations of stochastic delay equations: the Milstein scheme, The Annals of Probability 32 (1A) (2004) 265–314.doi:10.1214/aop/1078415836

    work page doi:10.1214/aop/1078415836 2004

  30. [30]

    J. Ma, J. Yong, Forward-Backward Stochastic Differential Equations and Their Applications, Springer, Berlin, Ger- many, 2007. doi:10.1007/978-3-540-48831-6

    work page doi:10.1007/978-3-540-48831-6 2007

  31. [31]

    12949 of Lecture Notes in Computer Science, Springer, Cham, 2021, pp

    D.Conte, R.D’Ambrosio, G.Giordano, B.Paternoster, ContinuousextensionofEuler–Maruyamamethodforstochas- tic differential equations, in: Computational Science and Its Applications – ICCSA 2021, Vol. 12949 of Lecture Notes in Computer Science, Springer, Cham, 2021, pp. 135–145.doi:10.1007/978-3-030-86653-2_10

    work page doi:10.1007/978-3-030-86653-2_10 2021

  32. [32]

    Fodor, H

    G. Fodor, H. T. Sykora, D. Bachrathy, A collocation method for stochastic delay differential equations, Probabilistic Engineering Mechanics 74 (2023) 103515.doi:10.1016/j.probengmech.2023.103515. 22

    work page doi:10.1016/j.probengmech.2023.103515 2023