% explicitfd_Euro_call.m
% Explicit finite difference pricing of a European call option, see
% page 60 of Clewlow and Strickland "Implementing derivatives models".
clear
K=100
S=100 % stock price at time 0
sig=0.2 % volatility
T=1 % end-time (in years)
r=0.06 % risk-free interest rate
div=0.03 % dividend yield
N=3000 % number of time steps
dx=0.025 % space step
Nj=200 % number of space steps

dt=T/N
if dt/dx^2>1/(3*sig^2)  disp('stability warning');  end
nu=r-div-0.5*sig^2
edx=exp(dx)
pu=0.5*dt*((sig/dx)^2+nu/dx) % note changed from tree_Euro_call.m
pm=1.0-dt*(sig/dx)^2-r*dt % note changed from tree_Euro_call.m
pd=0.5*dt*((sig/dx)^2-nu/dx) % note changed from tree_Euro_call.m
% note no explicit discount factor needed
tic
% initialize asset prices at maturity:
joffset=Nj+1; % used to translate from tree j coordinate to MATLAB vector index via index=j+joffset
St(-Nj+joffset) = S*exp(-Nj*dx);
for j=-Nj+1:Nj
    St(j+joffset)=St(j-1+joffset)*edx;
end 

% initialize option values at maturity:
C=zeros(N,2*Nj+1);
for j=-Nj:Nj
    C(N,j+joffset) = max( 0 , St(j+joffset)-K );
end

% step back through lattice:
for i=N-1:-1:1
    for j=-Nj+1:Nj-1
        C(i,j+joffset) = pu*C(i+1,j+1+joffset) + pm*C(i+1,j+joffset) + pd*C(i+1,j-1+joffset);
    end   % for j
    % boundary conditions:
    C(i,-Nj+joffset) = C(i,-Nj+1+joffset);
    C(i,Nj+joffset) = C(i,Nj-1+joffset) + (St(Nj+joffset)-St(Nj-1+joffset));
end % for i
% call value equals C with i=0 and j=0:
European_call=pu*C(1,0+1+joffset) + pm*C(1,0+joffset) + pd*C(1,0-1+joffset)
toc