# Problem with Richardson extrapolation method for weak convergence in SDE

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.

• 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 – Kirill Mar 2 '16 at 13:52
• 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! – nabla Mar 2 '16 at 14:58
• @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 :) – Christian Clason Mar 2 '16 at 15:03
• Sure thing (and some more random words so I can answer). – nabla Mar 2 '16 at 15:10

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$.
• 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. – nabla Mar 2 '16 at 15:26