% tree_Euro_call.m
% Trinomial tree pricing of a European call option, see
% page 54 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=3 % number of time steps
dx=0.2 % space step

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*((sig^2*dt+nu^2*dt^2)/dx^2+nu*dt/dx)
pm=1.0-(sig^2*dt+nu^2*dt^2)/dx^2
pd=0.5*((sig^2*dt+nu^2*dt^2)/dx^2-nu*dt/dx)
disc=exp(-r*dt)

tic

% initialize asset prices at maturity:
joffset=N+1; % used to translate from tree j coordinate to MATLAB vector index via index=j+joffset
St(-N+joffset) = S*exp(-N*dx);
for j=-N+1:N
    St(j+joffset)=St(j-1+joffset)*edx;
end 

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

% step back through lattice:
for i=N-1:-1:1
    for j=-i:i
        C(i,j+joffset) = disc*(pu*C(i+1,j+1+joffset) + pm*C(i+1,j+joffset) + pd*C(i+1,j-1+joffset));
    end    
end

% call value equals C with i=0 and j=0:
European_call=disc*(pu*C(1,0+1+joffset) + pm*C(1,0+joffset) + pd*C(1,0-1+joffset))

toc