% explicitfd_Amer_put.m
% Explicit finite difference pricing of a American put option, see
% page 62 of Clewlow and Strickland "Implementing derivatives models".
% Here we have not implemented the storage efficiency improvement.
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=3 % number of time steps
dx=0.2 % space step
Nj=3 % 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) 
pm=1.0-dt*(sig/dx)^2-r*dt 
pd=0.5*dt*((sig/dx)^2-nu/dx) 
% 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 , K - St(j+joffset) ); % put payoff
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 for put:
    C(i,-Nj+joffset) = C(i,-Nj+1+joffset)+(St(-Nj+1+joffset)-St(-Nj+joffset));
    C(i,Nj+joffset) = C(i,Nj-1+joffset);
    % early exercise condition:
    for j=-Nj:Nj
        C(i,j+joffset) = max( C(i,j+joffset) , K - St(j+joffset));
    end % for j
end % for i
% put value equals C with i=0 and j=0:
American_put=max( pu*C(1,0+1+joffset) + pm*C(1,0+joffset) + pd*C(1,0-1+joffset) , K - St(0+joffset))
toc