function V = binom_am(rho,S0,u,d,N,H)
%   V0 = binom_am(rho,S0,u,d,N,H)
% find European option price using binomial tree method 
% input:
%     rho: interest rate per time period
%     S0:  initial stock price
%      u:  factor by which stock price may increase at each step
%      d:  factor by which stock price may decrease at each step 
%      N:  number of time periods to maturity
%      H:  payoff function, e.g. @(s)max(5-s,0) for put option with strike 5
% output:
%     V0:  American option price at time t=0
%
% example:
%   H = @(s) max(7-s,0);         % put option with strike K=7 
%   V0 = binom_am(.01,10,1.1,0.9,10,H)

q = (1+rho-d)/(u-d);             % risk-neutral probability
S = S0*d^N*(u/d).^(0:N);         % stock prices at time t_N
V = H(S);                        % option prices at time t_N given by payoff function H
for j=N-1:-1:0                   % find option prices at time t_j, j=N-1,...,0
  Vcont = (q*V(2:end)+(1-q)*V(1:end-1))/(1+rho); % risk-neutral expectation for continuing
  S = S0*d^j*(u/d).^(0:j);       % stock prices at time t_j
  V = max(H(S),Vcont);           % stopping gives H(S), continuing gives Vcont. Pick max.
end

