4
$\begingroup$

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.

enter image description here

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')
$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – davidhigh Sep 24 at 17:05
  • $\begingroup$ @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 $\endgroup$ – sound wave Sep 24 at 17:41
2
$\begingroup$

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

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ 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$? $\endgroup$ – sound wave Sep 25 at 8:57
  • 1
    $\begingroup$ 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. $\endgroup$ – Peter Frolkovič Sep 25 at 11:38
  • $\begingroup$ 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? $\endgroup$ – sound wave Sep 29 at 17:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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