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

% actual calculations
deltaT = T/N;
u=exp(sigma * sqrt(deltaT));
d=1/u;
p=(exp((r-D0)*deltaT) - d)/(u-d); % adjust for dividends
lattice = zeros(N+1,N+1);
discount = exp(-r*deltaT);
for j=0:N
   lattice(N+1,j+1)=max(0 , K-S*(u^j)*(d^(N-j)));
end
% fill in rest of lattice
for i=N-1:-1:0
   for j=0:i
      lattice(i+1,j+1) = max( K-S*u^j*d^(i-j) , ...
         discount *(p * lattice(i+2,j+2) + (1-p) * lattice(i+2,j+1)));
   end
end
price = lattice(1,1);