# Accuracy of Explicit Euler method (finite difference) decreases as Δx decreases, shouldn't it increase?

The price of a commodity can be described by the Schwartz mean reverting SDE $$dS = \alpha(\mu-\log S)Sdt + \sigma S dW, \qquad \begin{array}.W = \text{ Standard Brownian motion} \\ \alpha = \text{ strength of mean reversion}\end{array}$$

From it is possible to derive the PDE for the price of the forward contract having the commodity as underlying asset $$\tag1\frac{\partial F}{\partial t} + \alpha\Big(\mu-\frac{\mu-r}\alpha -\log S\Big)S\frac{\partial F}{\partial S}+\frac12\sigma^2S^2\frac{\partial^2F}{\partial S^2} = 0$$

whose analytical solution is

$$F(S,\tau)=\exp\bigg(e^{-\alpha\tau}\log S +\Big(\mu-\frac{\sigma^2}{2\alpha}-\frac{\mu-r}{\alpha}\Big)(1-e^{-\alpha\tau})+\frac{\sigma^2}{4\alpha}(1-e^{-2\alpha\tau})\bigg)$$

where $$\tau=T-t$$ is the time to expiry ($$T$$ is the time of delivery/expiry).

Using Euler explicit method, i.e. forward difference on $$\dfrac{\partial F}{\partial t}$$ and central difference on $$\dfrac{\partial F}{\partial S}$$ and $$\dfrac{\partial^2F}{\partial S^2}$$, we can discretize eq (1) as $$F^{n+1}_i = a F^n_{i-1} + b F^n_i + c F^n_{i+1}$$ where $$a = \dfrac{S\Delta t}{2\Delta S}\bigg(\alpha\mu-(\mu-r)-\alpha\log(S)-\dfrac{\sigma^2S}{\Delta S}\bigg)$$

$$b = \bigg(1-\sigma^2S^2\dfrac{\Delta t}{\Delta S^2}\bigg)$$ and $$c = \dfrac{S\Delta t}{2\Delta S}\bigg(-\alpha\mu+(\mu-r)+\alpha\log(S)-\dfrac{\sigma^2S}{\Delta S}\bigg)$$.

To run Explicit Euler we have then to choose the number $$N$$ of time steps, which also set $$\Delta t$$ since $$\Delta t = T/N$$, and the size of $$\Delta S$$. Since the finite difference scheme divides the cartesian plane (time is the X-axis, and spot price is the Y-axis) in a grid, if we take more time steps and/or smaller $$\Delta S$$ the grid will be more dense and the accuracy of the approximation should increase.

However, the code I wrote using the equations above doesn't work in this way, in particular to have big accuracy I have to use a large $$\Delta S$$, and the accuracy decreases when using small values of $$\Delta S$$ to point that by using $$\Delta=0.1$$ the relative error explodes to $$10^{165}$$ as you can see in the image below (dS stands for $$\Delta S$$).

Even if my question involves finance topics, I think the problem is purely numerical or due to a wrong discretization formula, this is why I asked on scicomp.

Here is the matlab code if you would like to inspect it

%% Data and parameters
spot_prices = [ 22.93 15.45 12.61 12.84 15.38 13.43 11.58 15.10 14.87 14.90 15.22 16.11 18.65 17.75 18.30 18.68 19.44 20.07 21.34 20.31 19.53 19.86 18.85 17.27 17.13 16.80 16.20 17.86 17.42 16.53 15.50 15.52 14.54 13.77 14.14 16.38 18.02 17.94 19.48 21.07 20.12 20.05 19.78 18.58 19.59 20.10 19.86 21.10 22.86 22.11 20.39 18.43 18.20 16.70 18.45 27.31 33.51 36.04 32.33 27.28 25.23 20.48 19.90 20.83 21.23 20.19 21.40 21.69 21.89 23.23 22.46 19.50 18.79 19.01 18.92 20.23 20.98 22.38 21.78 21.34 21.88 21.69 20.34 19.41 19.03 20.09 20.32 20.25 19.95 19.09 17.89 18.01 17.50 18.15 16.61 14.51 15.03 14.78 14.68 16.42 17.89 19.06 19.65 18.38 17.45 17.72 18.07 17.16 18.04 18.57 18.54 19.90 19.74 18.45 17.33 18.02 18.23 17.43 17.99 19.03 18.85 19.09 21.33 23.50 21.17 20.42 21.30 21.90 23.97 24.88 23.71 25.23 25.13 22.18 20.97 19.70 20.82 19.26 19.66 19.95 19.80 21.33 20.19 18.33 16.72 16.06 15.12 15.35 14.91 13.72 14.17 13.47 15.03 14.46 13.00 11.35 12.51 12.01 14.68 17.31 17.72 17.92 20.10 21.28 23.80 22.69 25.00 26.10 27.26 29.37 29.84 25.72 28.79 31.82 29.70 31.26 33.88 33.11 34.42 28.44 29.59 29.61 27.24 27.49 28.63 27.60 26.42 27.37 26.20 22.17 19.64 19.39 19.71 20.72 24.53 26.18 27.04 25.52 26.97 28.39 ];
S = spot_prices; % real data

r = 0.1;    % yearly instantaneous interest rate
T = 1/2;   % expiry time

alpha = 0.069217; %
sigma = 0.087598; % values estimated from data
mu = 3.058244;    %

%% Exact solution
t = linspace(0,T,numel(S));
tau = T-t; % needed in order to get the analytical solution (can be seen as changing the direction of time)
F = exp( exp(-alpha*tau).*log(S) + (mu-sigma^2/2/alpha-(mu-r)/alpha)*(1-exp(-alpha*tau)) + sigma^2/4/alpha*(1-exp(-2*alpha*tau)) ); % analytical solution
F(1) = 0; % I think since there is no cost in entering a forward contract
plot(t,S)
hold on
plot(t,F,'g')
Exact_solution = F;

%% Explicit Euler approximation of the solution
S1 = S(2:end-1);  % all but endpoints
N = 3000; % number of time steps
dt = T/N; % delta t
dS = 1e1; % delta S, by decreasing dt and/or dS the approximation should improve
for m = 1:N
F(2:end-1) = S1*dt/2/dS.*( alpha*mu-(mu-r)-alpha*log(S1)-sigma^2*S1/dS).*F(1:end-2) ...
+                                  (1+sigma^2*S1.^2*dt/dS^2).*F(2:end-1) ...
+ S1*dt/2/dS.*(-alpha*mu+(mu-r)+alpha*log(S1)-sigma^2*S1/dS).*F(3:end);
F(1) = 0; % correct?
F(end) = S(end); % correct?
end
plot(t,F,'r.')
legend('Spot prices','Forward prices from exact solution','Forward prices from Explicit Euler')
title("dS = " + dS + ", relative error = " + norm( F-Exact_solution,2 ) / norm( Exact_solution,2 ))
xlabel('time')
ylabel('price')

• Have a look at the Courant–Friedrichs–Lewy condition. It basically says that you can't simply decrease $\Delta S$ to get a better accuracy for explicit Euler, but you have to decrease $\Delta T$ accordingly. If that bothers you, take a look at other solvers which are unconditionally stable, e.g.implicit Euler. – davidhigh Sep 24 at 17:05
• @davidhigh Thank you for the reply, I don't understand what is in this case the magnitude of velocity appearing in the formula for 1 dimensional CFL, anyway since $\Delta t=0.5/3000$ and $\Delta S=0.1$, their ratio is about 0.0016 way smaller than 1, so I think that CLF is satisfied. On this paper page 28 web.maths.unsw.edu.au/~qlegia/thesisfinal.pdf it is written that the explicit finite difference scheme is stable if and only if $0<\dfrac{\Delta t}{(\Delta S)^2}\le\dfrac12$, which is also true with the previous values – sound wave Sep 24 at 17:41

Looking briefly to your implementation, it seems to me that you may not discretize correctly the (independent) variable $$S$$. The coefficient of your PDE are depending on $$S$$, so you should evaluate them using correct (discretizied) values. Choosing smaller values od $$\Delta S$$ requires more discrete values of $$S$$, you can not choose a fixed number of "prices" $$S$$.

Al least I would define your coefficients $$a, b, c$$ as follows:

$$a = \dfrac{S_i\Delta t}{2\Delta S}\bigg(\alpha\mu-(\mu-r)-\alpha\log(S_i)-\dfrac{\sigma^2 S_i}{\Delta S}\bigg)$$

$$b = \bigg(1-\sigma^2S^2_i\dfrac{\Delta t}{\Delta S^2}\bigg)$$ and $$c = \dfrac{S_i\Delta t}{2\Delta S}\bigg(-\alpha\mu+(\mu-r)+\alpha\log(S_i)-\dfrac{\sigma^2S_i}{\Delta S}\bigg)$$,

where

$$S_{i+1}=S_i + \Delta S$$

and you have to define properly $$S_0$$ and the $$\Delta S$$, so that the interval of your interest for $$S$$, say, $$(S_{\min},S_{\max})$$ is covered.

If you are interested in values of $$F$$ for particular values of $$S$$, you can interpolate from two neighbor values, say $$S_i$$ and $$S_{i+1}$$ for properly chosen $$i$$.

• Thank you for the help, yes you are right about the fact that I did not discretize the variable $S$. I did not do that because I'm dealing with market real data (200 prices) which are indeed already discrete, so I don't know how to discretize them. In my mind when I think about discretization I think of an interval of real numbers, say [0,T], which has to be split into many subintervals. Since market data is a finite set, to discretize it we have to manually add other values by taking the mean of two consecutive values? For example with $S_1$ and $S_2$ we can find a middle value $(S_1-S_2)/2$? – sound wave Sep 25 at 8:57
• You need not to discretize your market data for S, in fact, they should not be used at all in the numerical solution of your PDE!! As S is independent variable, you should consider it in some real interval (for which your market prices are inside) and then discretize this interval depending on $\Delta S$. Then you compute the numerical solution $F_i^n$ that approximate the exact values $F(S_i,t^n)$. Of course, your discrete values S_i will be different to your market values, but to evaluate numerical solution everywhere, use an interpolation. – Peter Frolkovič Sep 25 at 11:38
• I did what you said and it works now thanks! But wouldn't be better to construct the vector of prices to be $(..., S_0-2\Delta t, S_0-\Delta t, S_0, S_0+\Delta t, S_0+2\Delta t)$ where $S_0=22.93$ is the first spot price of my real data? In this way there would be no need to interpolate and the result will be better I guess. What do you think? – sound wave Sep 29 at 17:48