A naïve approach to solving stochastic differential equations (SDEs) would be:

  • take a regular multi-step Runge–Kutta method,
  • use a sufficiently fine discretisation of the underlying Wiener process,
  • make each step of the Runge–Kutta method analogous to an Euler–Maruyama.

Now, this fails on multiple levels and I understand why. However, now I am tasked to convince people of this fact who have little knowledge of Runge–Kutta methods and stochastic differential equations to begin with. All the arguments I am aware of are nothing I can communicate well in the given context. Hence, I am looking for an easily understandable argument that the above approach is doomed.

  • $\begingroup$ @BiswajitBanerjee: I am aware of this and I do indeed not claim that I have understood this to the deepest possible extent. Still I do not think that providing all the arguments here will improve the answer as those who can provide an answer are aware of them. Moreover, this case is somewhat special as it is about explaining why something does not work, for which there naturally are many answers, starting with “we tested it and it failed”. $\endgroup$
    – Wrzlprmft
    Commented Jul 28, 2017 at 11:23
  • $\begingroup$ I wasn't talking about experts on stochastic ODEs but the average reader who understands random variables and RK when I said "us". However, I won't bother you further if you don't want to provide an example of your thinking. $\endgroup$ Commented Jul 28, 2017 at 21:30

1 Answer 1


Let's take a stochastic differential equation:

$$ X_t = f(t,X_t)dt + g(t,X_t)dW_t $$

Here's a few different arguments which lead to intuitive understandings of why the mathematics behind the higher order methods is necessary. I will be discussing in terms of strong order, which is the same as saying "for a given Brownian motion $W(t)$, how well does the numerical integral solve that trajectory?"

Regularity of the Equation

First of all, your proposed method fails to take into account the fact that $X_t$ is not continuously differentiable. Actually, you can use Rossler's results to show that extending normal RK methods like you suggested will result in convergent methods, but they will only have strong order 0.5. The reason is because they were derived using calculus with $X_t$ being differentiable and having a Taylor series. Brownian motion is not differentiable, and instead has Holder continuity of $\alpha < 0.5$ a.s.

However, like in perturbation theory, processes which are not regular enough are not expandable in terms of a Taylor series, but with Holder regularity $\alpha$ they can be expanded in terms of a Puiseux series with terms of $\alpha$, i.e. for Brownian motion there is an extension to the Taylor series notion which is expanded in terms of something like $\frac{1}{2}$ derivatives. Like in regular calculus, the first term is the "linear term", i.e. change $dt$ to $\Delta t$ and $dW_t$ to $N(0,dt)$ and you get something about right. This is why the methods, including things like Euler-Maruyama, converge with strong order 0.5: they get the first term in the Taylor series correct. However, the higher order terms need to have the corrections for the fact that $X_t$ is not continuously differentiable, which is why normal methods fail to do so.

Instantaneous Correlations and Iterated Integrals

That's a quick heuristic explanation, but there's a bit more to it. Let's look at a few other details. A Taylor series is not just the expansion in terms of derivatives, but it can also be thought of as the number higher order terms to integrate. $X_t = X_0 + \Delta t f(t,X_t)$ is integrating once. But if you add the $dt^2$ term, to get the units right you need to do double integrals. $dt^2$ is easy to integrate twice, but what is $dW_t^i dW_t^j$? These are the instantaneous correlations between the Brownian motions. You need to know this to compute the double integral. If you're only looking at averages, you can lop this off. But in any trajectory there are correlations between the different Brownian motions of a system of differential equations. Assuming there are no correlations between the Brownian motions is another way of characterizing the Maruyama extension of deterministic methods, but to get the next term in the series (the 1.0 term) you have to get this right. The Milstein correction is precisely adding these correlation terms. When noise is diagonal, this is equivalent to understanding that there is no correlation except with itself, but correlation with ones self is just the variance which is $dt$, and so there must be a correction of $dW_t^2$ vs $dt$, i.e. $dW^2 - dt$. When there is non-diagonal noise, these double integrals have to be approximated so that the instantaneous correlations of the Brownian motions are taken into account, and the common approximation here is the Wiktorsson approximation which is what then makes non-diagonal noise simulations so complicated (since there is no analytical solution even to the double integrals).

Average Effect of Diffusion

But this leads us to another way of thinking about the problem. Thinking of expanding in terms of moments, in some heuristic sense the first order term, the strong order 1.0 or $\mathcal{O}(\Delta t)$ term, must get the average movements correct, right? Here's a question: what's the derivative of $g$ in time? The easiest answer would be to define the derivative the normal way:

but this isn't actually correct when putting $g$ into the context of the SDE. If we think about the derivative of $g$ in terms of how much it changes $X_t$, it's not always on average pointing in the same direction since it's always multiplied by this random factor $dW_t$. The question is: what is the average size of this $dW_t$? Diffusion has changes on average on the scale of $\sqrt{\Delta t}$, so in reality the affect that $g(t,X_t)$ has is more like

$$ \frac{g(t+\Delta t,X_{t+\Delta t}) - g(t,X_t)}{\sqrt{\Delta t}} $$

You can show more rigorously that the numerical derivative should be this with $X_{t + \Delta t} = X_t + g(t,X_t)\sqrt{\Delta t}$ as the "predictor forward in time".

But intuitively, this is just understanding the average effect that $g$ has on the trajectory of $X_t$: about $g(t,X_t)\sqrt{\Delta t}$. In a Runge-Kutta method, an internal step at time $c_i$ is supposed to be an approximation of the value of $X_{t + c_i\Delta t}$, but even from this quick physical heuristic argument about diffusion we see that the easy extension of a Runge-Kutta method is already wrong on average: it's wrong by about $g(t,X_t)\sqrt{c_i \Delta t}$ which is another way to explain why it's at most strong order 0.5 (it's surprising that methods still work! But you can attribute this to the fact that the sum of the stages in an RK method must be 1, and so this error is somewhat canceled out). Interestingly, this heuristic argument goes pretty deep, since higher order stochastic Runge-Kutta methods like those due to Rossler have corrections which are precisely related to $g(t,X_t)\sqrt{\Delta t}$.


Those are 3 different heuristic ways to understand why higher orders must involve stochastic calculus. Higher orders must take into account the fact that Holder regularity is 1/2 and thus there are extra terms in the Taylor series, they must take into account the instantaneous correlations, and they must at least take into account the average effects of the diffusion term. Otherwise they are doomed to not be correct to $\mathcal{O}(\Delta t)$, and instead only satisfy the "linear approximation" of the first term and receive $\mathcal{O}(\sqrt{\Delta t})$.

Of course, in some circumstances there are ways to find appropriate generalizations that do give higher order methods, but I'll leave this as a dangling thread because that's one point of a paper I'll be submitting soon. Hope this helps.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.