% implicitfd_Amer_put.m
% Implicit finite difference pricing of a American put option, see
% page 69 of Clewlow and Strickland "Implementing derivatives models".
% Includes the storage efficiency improvement as in the code of
% explicitfd_Amer_put2.m
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

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) 
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(2,2*Nj+1);
ioffset=1; % use indices 0 and 1 of pseudo-code, adding ioffset to each.
for j=-Nj:Nj
    C(0+ioffset,j+joffset) = max( 0 , K - St(j+joffset) ); % put payoff
end
% compute derivative boundary conditions:
lambda_L = -1*( St( -Nj+1+joffset) - St(-Nj+joffset) );
lambda_U = 0.0;
% step back through lattice:
for i=N-1:-1:0

    solve_implicit_tridiagonal_system; % calls triadiagonal solver
    
    % early exercise condition:
    for j=-Nj:Nj
        C(0+ioffset,j+joffset) = max( C(1+ioffset,j+joffset) , K - St(j+joffset));
    end % for j
end % for i
% put value equals C with i=0 and j=0:
American_put=C(0+ioffset,0+joffset) 
toc