function [price, lattice] = LatticeEurCall(S,K,r,T,sigma,N,D0)
% usage: [price, lattice] = LatticeEurCall(S,K,r,T,sigma,N,D0)
% description: Returns the call price of an European option calculated by
% the binomial lattice method where S is the initial security price, K is
% the strike price, r is the annual interest rate, T is the time until
% maturity in years, sigma is the annual risk, and N is the number of
% time steps dt = T/N until maturity used to calculate the call price,
% and D0 is the (optional) dividend rate with default zero.
% Optionally returns the entire lattice as a second return argument.

% local variables:
% j: index variable
% i: index variable
% deltaT: size of a time step (in years)
% u: the probability of going up on the lattice
% p: the risk free probability
% lattice: the price of the stock at any point in time

% input checks
if (nargin < 6)
  error('LatticeEurCall: need six arguments: S,K,r,T,sigma,N');
end
if (length(S) ~= [1,1])
  error('LatticeEurCall: S must be a single number');
end
if (length(K) ~= [1,1])
  error('LatticeEurCall: K must be a single number');
end
if (length(r) ~= [1,1])
  error('LatticeEurCall: r must be a single number');
end
if (length(T) ~= [1,1])
  error('LatticeEurCall: T must be a single number');
end
if (length(sigma) ~= [1,1])
  error('LatticeEurCall: sigma must be a single number');
end
if (length(N) ~= [1,1])
  error('LatticeEurCall: N must be a single number');
end
if (S < 0)
  error('LatticeEurCall: S must be a positive number');
end
if (T < 0)
  error('LatticeEurCall: T must be a positive number');
end
if (sigma < 0)
  error('LatticeEurCall: sigma must be a positive number');
end
if (N < 0 || N > 1000)
  error('LatticeEurCall: N must be a positive number not greater than 1000');
end
if (nargin < 7)
  D0 = 0;
end

% actual calculations
deltaT = T/N;
% a technically more accurate formula... 
% B = 0.5*(exp((D0-r)*deltaT) + exp((-D0+r+sigma^2)*deltaT));
% u=B + sqrt(B^2-1);
u = exp(sigma*sqrt(deltaT));
d=1/u;
rtilde = r - D0; % adjust for dividend rate
p = (exp(rtilde*deltaT)-d)/(u-d);
lattice = zeros(N+1,N+1);
discount = exp(-r*deltaT);
% construct the payoff curve at lattice points of expiry
for j=0:N
   lattice(N+1,j+1)=max(0 , S*(u^j)*(d^(N-j)) - K);
end
% fill in rest of lattice
for i=N-1:-1:0
   for j=0:i
      lattice(i+1,j+1) = discount * ...
         (p * lattice(i+2,j+2) + (1-p) * lattice(i+2,j+1));
   end   
end
price = lattice(1,1);
