Rabelais


Total Posts: 12 
Joined: Sep 2020 


The price of a commodity can be described by the Schwartz mean reverting SDE
where W is the standard Brownian motion and alpha is the strength of mean reversion. From it is possible to derive the PDE for the price of the forward contract having the commodity as underlying asset
(1) whose analytical solution is
where is the time to expiry (T is the time of delivery/expiry). Using Euler explicit method, i.e. forward difference on and central difference on and , we can discretize eq (1) as
where , and .
To run Explicit Euler we have then to choose the number N of time steps, which also set Δt since Δt = T/N, and the size of ΔS. Since the finite difference scheme divides the cartesian plane (time is the Xaxis, and spot price is the Yaxis) in a grid, if we take more time steps and/or smaller ΔS the grid will be more dense and the accuracy of the approximation should increase (provided that ).
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 ΔS, and the accuracy decreases when using small values of ΔS to point that by using ΔS=0.1 the relative error explodes to 10^165 as you can see in the image below (dS stands for ΔS).
If you want you can inspect the matlab code (see attachment), I'm sure there is an error somewhere but cannot find it. I think one problem is the fact that I'm using real data (vector spot_prices in the code) and so I cannot discretize along the yaxis, right?
Attached File: accuracy_test.rar




bullero


Total Posts: 75 
Joined: Feb 2018 


Does your implementation check if the condition is satisfied or not?
For any positive dt letting dS to go to zero (from right) will cause this (stability?) condition to explode above 0.5. 


Rabelais


Total Posts: 12 
Joined: Sep 2020 


Yes the condition is satisfied, when Δt = 0.5/3000 and ΔS = 0.1 then Δt/(ΔS)^2 = 0.0166... < 1/2, however the relative error explodes to 10^156 as shown in the image 



nikol


Total Posts: 1345 
Joined: Jun 2005 


I would plot a, b, c as function of time to be 100% sure and to spot the source of the problem.
log(S) around zero is huge. If S is thrown into negative region, then log(S) gets complex, but Matlab will continue delivering Re(value) without crash. 
... What is a man
If his chief good and market of his time
Be but to sleep and feed? (c) 

Jim


Total Posts: 166 
Joined: Jun 2004 


Since you are using reals instead of double precision, you are bumping up into machine precision errors.
In a nutshell, with infinite precision taking a smaller and smaller dS when computing a finite difference derivative gets you closer to the true numerical derivative. But computers only have finite precision. A floating point representation of X can only be represented in the computer as (X + e1), where e1 is the error. So instead of computing f' = (f(x + h)  f(x))/h, you really are getting f' =((f(x + h) + e1)  (f(x) + e2))/h. So as h>0, you get a true f' + (e1  e2)/h. Since e1 and e2 don't disappear as h>0, your finite difference derivative gets better for a while then it gets worse due to machine precision.
And you notice the error gets worse when S is large. That is because to represent a larger S there are fewer bits to represent the fraction to the right of the decimal point.



