Heston model's PDE - stability concerns

[Yani Boshev]

I am trying to implement implicit Euler scheme for solving the PDE of Heston stochastic volatility model for European options. I use operator splitting scheme discretizing the spatial variables (stock price and volatility) using implicit Euler and explicit Euler for cross derivative. Unfortunately, when I reduce the length of the time step I receive meaningless results around the upper stock price boundary (for example lower price for call option when the volatility is higher). If the time step is larger the results seem logical (not sure whether they are close to the exact solution but it seems so). Using the discrete maximum principle I obtained the following equation regarding stability of the scheme:
((1+e1*dt)/(1-e2*dt))
e1 - eigenvalues of the matrix containing the explicit terms (cross derivative coefficients)
e2 - eigenvalues of the matrices containing the implicit terms (coefficients of the derivatives of the prices and volatilities)
dt - time step

The problem is that eigenvalues of the matrix with the volatility terms are less in absolute value than the matrix containing the explicit terms making the equation bigger than 1 -> instable scheme. I tried to use exponential fitting factor for general convection-diffusion equation as proposed by Dr. Duffy ( FDM in Financial Engineering) but the results didn't change. I suppose that the cross derivative alters the exponential fitting factor but I don't know how to derive the new fitting factor in case I am right.

As guides I use the aforementioned book of Dr. Duffy and Quantitative Methods in Derivatives Pricing - Domingo Tavella. If you have any suggestions regarding the improvement of the stability of the scheme or have ideas about my mistakes, please comment. If you have more information about the derivation of the exponential fitting factors please tell me (I didn't understand clearly this chapter of the book of Dr. Duffy).

 

[Daniel Duffy]

Please see thesis by Roelof Sheppard (does it all in detail)

http://www.datasimfinancial.com/forum/viewtopic.php?t=101

Here an article on exponential fitting

http://www.datasimfinancial.com/forum/files/fittingidoc_106.pdf

 

[Daniel Duffy]

The FD scheme in Sheppard's thesis is unconditionally stable. Have you checked and debugged your code?

I use operator splitting scheme discretizing the spatial variables (stock price and volatility) using implicit Euler and explicit Euler for cross derivative

Yanenko scheme for X derivs is the standard.

 

[Yani Boshev]

Well, as you can suppose it was my fault. For boundary conditions when \(S -> inf\) I use \(S-strike\) for call option. My mistake was that I discounted the resulting value at every time step to the present value at the \(t-1\) interval. When I removed the discount factor the results seemed to be fine.

However, something else take my attention. I am not sure if it is possible near the upper price boundary when the current volatility is lower than the assumed long term volatility the price of the option to be higher than the price of the option when the current volatility is at the same rate as the assumed long term volatility. The results obtained by my code exhibit such pattern only near the upper price boundary so is this meaningful or I have to search for another mistake?

 

[Daniel Duffy]

Well, as you can suppose it was my fault. For boundary conditions when \(S -> inf\) I use \(S-strike\) for call option. My mistake was that I discounted the resulting value at every time step to the present value at the \(t-1\) interval. When I removed the discount factor the results seemed to be fine.

However, something else take my attention. I am not sure if it is possible near the upper price boundary when the current volatility is lower than the assumed long term volatility the price of the option to be higher than the price of the option when the current volatility is at the same rate as the assumed long term volatility. The results obtained by my code exhibit such pattern only near the upper price boundary so is this meaningful or I have to search for another mistake?

Difficult to say until you give more _precise_ details of FD scheme that you are using.

One can never rule out code errors as well.

 

[Yani Boshev]

Here is my code. Written in Matlab. I am using implicit Euler with exponential fitting factors and explicit Euler for the cross derivative. Also I am using Yanenko scheme. Hope it's readable.

clear all; close all; clc;
str=50;    %strike
n=100;      %number of spatial points
time=1;   
tspan=0:0.01:time;      %vector with the time steps
dt=tspan(2)-tspan(1);  %delta t
dt=dt/2;                %according to the Yanenko scheme solution of the first...
                        %equation is in period t+0.5
 
r=0.05;                %risk-free rate
vega=0.05;              %volatility of volatility
y=0.1;                  %long-term volatility level
om=0.018;              %rate of mean-reversion
rho=0.5;                %correlation   
mpr=0.0125;            %market price of risk
mprV=0.01;              %market price volatility risk
   
S=linspace(0, 100, n+2);   
S1(1:n,1)=S(2:(n+1));      %vector with length "n" containing the price values
ds=S(2)-S(1);              %delta price
   
Vol=linspace(0, 0.5, n+2);
Vol1(1:n,1)=Vol(2:n+1);    %vector with length "n" containing the volatility values
dv=Vol(2)-Vol(1);          %delta volatility
 
[s,vol]=meshgrid(S1,Vol1); 
s=reshape(s,n^2,1);        %vector with size n^2 containing price values
vol=reshape(vol,n^2,1);    %vector with size n^2 containing volatility values
 
 
effv=((om*(y-vol.^2)-vega*vol*(rho*mpr+sqrt(1-rho^2)*mprV)).*dv)./(2)...
    .*coth(((om*(y-vol.^2)-vega*vol*(rho*mpr+sqrt(1-rho^2)*mprV)).*dv)./(2*0.5*vega^2*vol.^2));  %exponential fitting factor volatility
 
effs=((r*s)*ds/2).*coth((r*s).*ds./(2*0.5*vol.^2.*s.^2));      %exponential fitting factor price
 
a=(effs./ds^2+0.5*r*s/ds)*dt;  %upper diagonal matrix A (point S+1,vol)
b=1+(2*effs./ds^2+0.5*r)*dt;    %main diagonal matrix A (point S,vol)
c=(effs./ds^2-0.5*r*s/ds)*dt;  %lower diagonal matrix A (point S-1,vol)
 
f=(0.5*rho*vega*vol.^2.*s/(4*ds*dv))*dt; %coef. of the cross der. Explicit Euler
 
d=(effv./dv^2+0.5*(om*(y-vol.^2)-vega*vol*(rho*mpr+sqrt(1-rho^2)*mprV))/dv)*dt;    %upper diagonal matrix C (point S,vol+1)
e=(1+(2*effv./dv^2+0.5*r)*dt);      %main diagonal matrix C (point S, vol)
j=(effv./dv^2-0.5*(om*(y-vol.^2)-vega*vol*(rho*mpr+sqrt(1-rho^2)*mprV))/dv)*dt;  %lower diagonal matrix C (point S, vol-1)
 
A=(spdiags([-c b -a], [-n 0 n], n^2, n^2));    %building matrix A
 
f1=f;      %diagonal matrix B (cross der. Explicit Euler (point S-1, vol-1))
f2=-f;      %diagonal matrix B (cross der. Explicit Euler (point S-1, vol+1))
f3=-f;      %diagonal matrix B (cross der. Explicit Euler (point S+1, vol-1))
f4=f;      %diagonal matrix B (cross der. Explicit Euler (point S+1, vol+1))
 
%loop which will make the required positions in the vectors equal to zero
%these points will be implemented as Boundary conditions
for l=1:n
f1(l*n,1)=0;
f2((l-1)*n+1,1)=0;
f3(l*n,1)=0;
f4((l-1)*n+1,1)=0;
j(l*n,1)=0;
d((l-1)*n+1,1)=0;
end
 
mdiagB(1:n^2,1)=ones;  %main diagonal matrix B
B=spdiags([f1 f2 mdiagB f3 f4], [-(n+1) -(n-1) 0 n-1 n+1], n^2, n^2); %matrix B (cross der. Explicit Euler)
C=spdiags([-j e -d], [-1 0 1], n^2, n^2);        %matrix C (contains volatility part of the PDE)
 
 
 
%defining BC mat A B C
BCA(1:n^2,1)=zeros;    %vector boundary conditions matrix A
BCB(1:n^2,1)=zeros;    %vector boundary conditions matrix B
BCC(1:n^2,1)=zeros;    %vector boundary conditions matrix C
 
%loop required to fulfill the required positions in the BC vectors
for l=1:n
 
    BCA((n-1)*n+l,1)=-(((r*(S1(end,1)+ds)*ds/2)*coth(r*(S1(end,1)+ds)*ds/(2*0.5*Vol1(l,1)^2*(S1(end,1)+ds)^2)))/ds^2....
    +0.5*r*(S1(end,1)+ds)/ds)*dt*((S1(end,1)+ds)-str);      %boundary conditions points S (max)
 
    BCC(l*n,1)=-(((om*(y-(Vol1(end,1)+dv)^2)-vega*(Vol1(end,1)+dv)*(rho*mpr+sqrt(1-rho^2)*mprV))*dv/(2)...
    *coth(((om*(y-(Vol1(end,1)+dv)^2)-vega*(Vol1(end,1)+dv)*(rho*mpr+sqrt(1-rho^2)*mprV))*dv)/(2*0.5*vega^2*(Vol1(end,1)+dv)^2)))/dv^2....
    +(om*(y-(Vol1(end,1)+dv)^2)-vega*(Vol1(end,1)+dv)*(rho*mpr+sqrt(1-rho^2)*mprV))/dv)*dt*S1(l,1);      %boundary conditions points vol (max)
 
    if l==1
        BCB(l*n,1)=(0.5*rho*vega*(Vol1(end,1)+dv)^2*(S1(l,1)+ds)/(4*ds*dv))*dt*(S1(l,1)+ds); %BC cross derivative - S(2), vol(max)
   
    else if l==n
   
        for l2=1:n
            if l2==1
                BCB((l-1)*n+l2,1)=...
                (0.5*rho*vega*Vol1(l2+1,1)^2*(S1(end,1)+ds)/...
                (4*ds*dv))*dt*(S1(end,1)+ds-str);      %BC cross derivative - S(max), vol(2)
           
            else if l2==n
                  BCB((l-1)*n+l2,1)=...
                  (0.5*rho*vega*(Vol1(end,1)+dv)^2*(S1(end,1)+ds)/...
                  (4*ds*dv))*dt*(S1(end,1)+ds)-(0.5*rho*vega*Vol1(l2-1,1)^2*(S1(end,1)+ds)/...
                  (4*ds*dv))*dt*(S1(end,1)+ds-str)-(0.5*rho*vega*(Vol1(end,1)+dv)^2*(S1(l2,1)-ds)/...
                  (4*ds*dv))*dt*(S1(end,1)-ds);      %BC cross derivative - S(max), vol(max)
               
                else
                  BCB((l-1)*n+l2,1)=...
                  (0.5*rho*vega*Vol1(l2+1,1)^2*(S1(end,1)+ds)/...
                  (4*ds*dv)-(0.5*rho*vega*Vol1(l2-1,1)^2*(S1(end,1)+ds)/...
                  (4*ds*dv)))*dt*(S1(end,1)+ds-str);      %BC cross derivative - S(max) , vol(points between 1 and "n")
                end
            end
      end
       
        else
            BCB(l*n,1)=(0.5*rho*vega*(Vol1(end,1)+dv)^2*(S1(l+1,1))/...
            (4*ds*dv)*(S1(l+1,1))-(0.5*rho*vega*(Vol1(end,1)+dv)^2*(S1(l-1,1))/...
            (4*ds*dv))*(S1(l-1,1)))*dt;        %BC cross derivative - S(points between 1 and "n"), vol(max)
        end
    end
end
 
%defining init conditions CALL option
v0(1:n^2,1)=max(0,s-str);
 
%solution
for t=1:length(tspan)
 
v=A\(B*v0+BCB-BCA);
v0=v;
v=C\(B*v0+BCB-BCC);
v0=v;
end
 
v=reshape(v,n,n)

Now if I decide to add jump component should I generate 3-D mesh or I could generate mesh with log(S) and log(jump size) only. I think that generating 2-D mesh is more appropriate because jump process and stochastic volatility process are not correlated. I think to use the output at every step of the Heston model, Fourier transform it, multiply it with the FT of the jump distribution and then Inverse Fourier Transform. I should use the negative of the generated log(jump size). I suppose I should use log(S) in the non-integral part of the PIDE. Am I right?