I have implemented the Richardson extrapolation of the Euler-Maruyama method to 4th order, to estimate the moments of SDE.

The Euler-Maruyama works, and I would expect the Richardson extrapolation to work too, because when I suppress the noise contribution it gives the correct solution in the resulting ordinary differential equations.

However, when I add the noise term, I don't get the moments (the variance, in particular) I would expect, as I do with the regular Euler method, so apparently I am doing something wrong.

This is my MATLAB code for the extrapolation:

function [t, y] = Richardson_weak_4th_order(a, b, t_interval, y0, h)

[t1, y1] = euler_maruyama(a, b, t_interval, y0, h);
[~, y2]  = euler_maruyama(a, b, t_interval, y0, h/2);
[~, y3]  = euler_maruyama(a, b, t_interval, y0, h/4);
[~, y4]  = euler_maruyama(a, b, t_interval, y0, h/8);

% Extrapolation
y2 = y2(1:2:end,:);
y3 = y3(1:4:end,:);
y4 = y4(1:8:end,:);

t = t1;
y = 1/21*(64*y4 - 56*y3 + 14*y2 - y1);

For simplicity, I am testing with an Ornstein-Uhlenbeck process $$ \mathbb{d} X_t = (m - X_t)\mathbb{d}t + \sigma \mathbb{d}W_t $$ for which $\mathbb{E}(X) \simeq m$ and $\operatorname{Var}(X) \simeq \frac{\sigma^2}{2}$. The results are OK with the Euler method, but the variance is quite bigger with the extrapolation. What am I missing?

Edit: The problem is Richardson extrapolation was not implemented on the estimated moments (as it should), but to realizations of the process.

  • 1
    $\begingroup$ Are you applying extrapolation to the paths of the stochastic process themselves, or to the moments? Do you keep the same Brownian motion path for each of the EM paths? IMO, it's not terribly clear if Richardson extrapolation should work at all. Related question: scicomp.stackexchange.com/q/396/713 $\endgroup$
    – Kirill
    Mar 2, 2016 at 13:52
  • $\begingroup$ After writing the post I took a look at "Numerical solution of SDE through computer experiments" and found that, indeed, I was applying Richardson to the process realizations and not to the moments. Thanks! $\endgroup$
    – nabla
    Mar 2, 2016 at 14:58
  • $\begingroup$ @Kirill Would you mind writing up your comment as an answer (maybe in a bit more detail)? Otherwise the Community bot will keep poking the question :) $\endgroup$ Mar 2, 2016 at 15:03
  • $\begingroup$ Sure thing (and some more random words so I can answer). $\endgroup$
    – nabla
    Mar 2, 2016 at 15:10

1 Answer 1


You are applying Richardson extrapolation to increasingly-accurate independent sample paths of an OU process.

How would this work even for independent identically normally distributed random variables $y_1,y_2,y_3,y_4\sim\mathrm{N}(\mu,\sigma^2)$? These have already "converged" to the true distribution $\mathrm{N}(\mu,\sigma^2)$, and you might expect the extrapolation to just reproduce the variables, but $$ \tfrac1{21}(64y_4-56y_3+14y_2-y_1) \sim \mathrm{N}(\mu, 16.8\sigma^2). $$ If they are independent, that destroys any chance that this linear combination would do the right thing.

According to a cursory search on Google Scholar, Richardson extrapolation can indeed be helpful. But I think it's necessary to be careful about how you pick your samples for it, and how independent they are. The above argument with $y_i$ would apply as well when $y_i$ are estimates of a sample moment of the generated process paths: the variance of the sample moments will be amplified by the process of Richardson extrapolation, so this is something you need to handle carefully.

So Richardson extrapolation can only be applied, with a little care, to sample moments, but not at all not sample paths. For sample moments specifically, it would indeed help: independent $y_k$ would be distributed as $\mathrm{N}(\mu,2^{-2(k-1)}\sigma^2)$, and then $$ \tfrac1{21}(64y_4-56y_3+14y_2-y_1) \sim \mathrm{N}(\mu, 0.7\sigma^2). $$ If you were to do the same thing with $y_k$'s coming from the same sample, then by my calculations the variance would (for sample means) be $ 0.15\sigma^2$.

  • $\begingroup$ If you add that the method has to be applied to the moments (or any $\mathbb{E}(g(X_t))$ in general), I will accept your answer. $\endgroup$
    – nabla
    Mar 2, 2016 at 15:26
  • $\begingroup$ @gerd Please feel free to suggest edits to this answer. $\endgroup$
    – Kirill
    Mar 2, 2016 at 15:28
  • 1
    $\begingroup$ @gerd I looked it up in Glasserman's Monte Carlo Methods, it's section 6.2.4, page 360. $\endgroup$
    – Kirill
    Mar 2, 2016 at 16:05

Your Answer

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

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