0
$\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 decreses when using small values of $\Delta S$ to point that by using $\Delta S=0.1$ the relative error explodes to $10^{165}$ as you can see in the image below (dS stands for $\Delta S$).

enter image description here

Here is the matlab code, I think there is an error somewhere, moreover I'm not sure about the boundary conditions F(1) and F(end)

%% 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$
  • $\begingroup$ Without delving much into your question is there a possibility that machine precision is a limiting factor and, beyond that, larger errors are inescapable? $\endgroup$ – Attack68 Sep 24 at 12:27
  • $\begingroup$ @Attack68 I don't think so since matlab uses 16 digits of precision by default $\endgroup$ – sound wave Sep 24 at 12:29
  • $\begingroup$ are you confident on your implementation? Are there any unit tests you could do to test your code in another way? It seems to me that if error increases for smaller grid size above machine precision then the error might be the code.. cannot discount that possibility. $\endgroup$ – Attack68 Sep 24 at 12:35
  • $\begingroup$ @Attack68 The exact solution seems computed correctly (plot is ok), while I'm not sure about the euler method implementation even if the formula should be correct (did the passages many times to verify), but the evidence of some error in the code is the explosion of the relative error when ΔS=0.1, but i cannot find where the error is. What do you mean with unit tests? $\endgroup$ – sound wave Sep 24 at 12:44
  • $\begingroup$ Ok last comment looking at this, it seems that $\Delta S$ is a denominator, so that terms become larger as it becomes smaller. I haven't done the work but my instinct says something is wrong here and it might be discretization formula. $\endgroup$ – Attack68 Sep 24 at 13:02
1
$\begingroup$

I haven't inspect everything in detail but I think I have a useful remark. Since you know the analytical solution, you can compute from there the boundary conditions.

For $t = T$, you have $\tau = 0$. You get $F(t = T; \tau = 0) = S(T)$ when replacing in the analytical solution. Then, F(end) = S(end), as you stated in your code.

On the other hand, when $t = 0$, $\tau = T$. If you replace that in the analytical solution, you get the same expression that you wrote but replacing $S$ with $S(0)$ and $\tau$ with $T$ so, in your code, F(1) = 0 do not match with the analytical solution. Since I also believe that F(1) should be zero, there might be something that we are missing in the analytical solution...

Please, forgive me if I am mistaken, but I thought it would be a useful comment.

PS: comments below address the underlying problem.

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ I see... there is something I am missing. I would recommend reading section 2.2 Finite Difference Discretization from the Andersen and Piterbarg Interest Rate Modeling book. It seems that there you will find all the answers you are looking. For me, there is something wrong in the fact that you are solving only in the time grid and not in the rectangular domain, i.e., including the S discretization. Where are those spot prices coming from? Analytical solution of its SDE? $\endgroup$ – rvignolo Sep 24 at 13:51
  • 1
    $\begingroup$ thank you i will check the book asap! the spot prices are real data, is it possible to discretize also S when using real data or is possible only when using computer generated one? thanks again for support1 $\endgroup$ – sound wave Sep 24 at 13:57
  • 1
    $\begingroup$ Okay, so there is the problem. You have to solve the PDE by applying a discretization in both directions, $S$ and $t$. Forget about market data (use it for initial condition only). Sorry for not being able to point out everything, I am more familiar with MC methods! Hope it helps. I also edited my answer because I believe that the forward contract start value should also be zero. So there might be something in the analytical solution that is either wrong or misleading. Thank you! $\endgroup$ – rvignolo Sep 24 at 14:08
  • 1
    $\begingroup$ Sorry, I missed your comment! For Monte Carlo simulation, you have to simulate paths for $S$ by applying a discretization scheme to its SDE. In this case, you are going to use the Euler Maruyama scheme. For the second case, i.e. PDE method, you also don't have to use market data. You have to apply a discretization scheme to both axis, $S$ and $t$ and solve the PDE. Take a look at Andersen, believe me it is going to be really useful! $\endgroup$ – rvignolo Sep 24 at 19:41
  • 1
    $\begingroup$ Thank you, in this case since we have the analytical solution, so we simply use it to perform the monte carlo method, we don't need to use a discretization scheme. I'm now reading the book and trying to write the code for the discretization scheme to both axis, but I don't understand how can we use it to approximate the green path in the image above without using the market data? Thanks again $\endgroup$ – sound wave Sep 24 at 19:56

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.