function [V,av] = xbinom_am_call(rho,S0,u,d,N,K,m)
%    V = binom_am_call(rho,S0,u,d,N,K,m)
% find American call option price using extended binomial tree
% plot option price vs initial stock price
% plot option price vs strike 
%
% 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
%      K:  strike
%      m:  number of extra option prices for lower and higher values of S0
% output:
%      V:  vector of option prices for initial stock prices av*S0
%     av:  vector (u/d)^(-m:m)
%
% example: rho=.01, S0=10, u=1.1, d=0.9, N=10, K=11, m=3
%   S0=10; K=11;
%   [V,av] = xbinom_eu_call(.01,S0,1.1,.9,10,K,3);
%   [S0*av;V]      % table with initial stock prices and option prices for K=11
%   [K./av;V./av]  % table with strike prices and option prices for S0=10
q = (1+rho-d)/(u-d);             % risk-neutral probability
S = S0*d^N*(u/d).^(-m:N+m);      % stock prices at time t_N, m extra values above & below
V = max(S-K,0);                  % option prices at time t_N for call option
for j=N-1:-1:0                   % find option prices at time t_j, j=N-1,...,0
  S = S0*d^j*(u/d).^(-m:j+m);    % stock prices at time t_j
  Vcont = (q*V(2:end)+(1-q)*V(1:end-1))/(1+rho); % expectation for continuing
  V = max(max(S-K,0),Vcont);     % stopping gives H(S), continuing gives Vcont. Pick max.
end
av = (u/d).^(-m:m);              % vector with values (u/d)^k

figure(1); plot(S0*av,V,'.-')    % plot option price vs. initial stock price
xlabel('initial stock price S_0'); 
title(sprintf('American call option price for strike K=%g vs. S_0',K))
hold on
s = linspace(S0*av(1),S0*av(end),500);
plot(s,max(s-K,0),'r',s,max(s-K/(1+rho)^N,0),'r--')
hold off
legend('option price','payoff','payoff with discounted strike')

figure(2); plot(K./av,V./av,'.-') % plot option price vs. strike
xlabel('strike price K'); 
title(sprintf('American call option price for S_0=%g vs. K',S0))
hold on
Kv = linspace(K/av(1),K/av(end),500);
plot(Kv,max(S0-Kv,0),'r',Kv,max(S0-Kv/(1+rho)^N,0),'r--')
hold off
legend('option price','current payoff','current payoff with discounted strike')
