function [price, lattice] = LatticeAmCallMod(S,X,r,T,sigma,N,D0)
% usage: [price, lattice] = LatticeAmCallMod(S,X,r,T,sigma,N,D0)
% description: Returns the price of an American call option calculated by
% the binomial lattice method where S is the initial security price, X 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 put 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('LatticeAmCall: need six arguments: S,K,r,T,sigma,N');
end
if (length(S) ~= [1,1])
  error('LatticeAmCall: S must be a single number');
end
if (length(K) ~= [1,1])
  error('LatticeAmCall: K must be a single number');
end
if (length(r) ~= [1,1])
  error('LatticeAmCall: r must be a single number');
end
if (length(T) ~= [1,1])
  error('LatticeAmCall: T must be a single number');
end
if (length(sigma) ~= [1,1])
  error('LatticeAmCall: sigma must be a single number');
end
if (length(N) ~= [1,1])
  error('LatticeAmCall: N must be a single number');
end
if (S < 0)
  error('LatticeAmCall: S must be a positive number');
end
if (T < 0)
  error('LatticeAmCall: T must be a positive number');
end
if (sigma < 0)
  error('LatticeAmCall: sigma must be a positive number');
end
if (N < 0 || N > 1000)
  error('LatticeAmCall: N must be a positive number not greater than 1000');
end
if (nargin ~= 7)
  D0 = 0;
end

% actual calculations
deltaT = T/N;
M = 1 + sigma^2*deltaT +exp(2*(r-D0)*deltaT);
u=(M + sqrt(M^2 -4))/2;
d=1/u;
r = r - D0; % adjust for dividend rate
p=(exp(r*deltaT) - d)/(u-d);
lattice = zeros(N+1,N+1);
for j=0:N
   lattice(N+1,j+1)=max(0 , S*(u^j)*(d^(N-j)) - X);
end
for i=N-1:-1:0
   for j=0:i
      % adjust for American call (stay above payoff curve at all times)
      lattice(i+1,j+1) = max(S*u^j*d^(i-j)-X, exp(-r*deltaT) * ...
         (p * lattice(i+2,j+2) + (1-p) * lattice(i+2,j+1)));
   end   
end
price = lattice(1,1);