%      Desmond Higham's Program
%
% Implements binomial method for European call

%%%%%%%%%% Problem and method parameters %%%%%%%%%%%%%
S = 3; E = 2; T = 1; r = 0.05; sigma = 0.3; 
M = 400; dt = T/M; p =0.5;
u = exp(sigma*sqrt(dt) + (r-0.5*sigma^2)*dt); 
d = exp(-sigma*sqrt(dt) + (r-0.5*sigma^2)*dt); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Time T option values 
W = max(S*d.^([M:-1:0]').*u.^([0:M]')-E,0);

% Work back to option value at time zero
for i = M:-1:1
    W = exp(-r*dt)*(p*W(2:i+1) + (1-p)*W(1:i));
end

disp('Option value is'), disp(W)