Problem with Richardson extrapolation method for weak convergence in SDE

Question

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.

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$.