function [valvec,cutvec]=amerval(s0,x,tau,rfr,coc,sig)
%  [valvec,cutvec]=amerval(s0,x,tau,rfr,coc,sig)
%
%  Add brief description for function in line above for Matlab help searches
%
% NOTE: YOU MUST ISSUE THE COMMAND: "run  amervald.m"
%       TO INITIALIZE GLOBAL VARIABLES BEFORE CALLING THIS PROCEDURE!!! 
%
% program to price American Put or Call options 
%   using analytical approximations as in Barone-Adesi and Whaley
%
%   see also: ameriv.g
%
%{value,cutoffs}=amerval(s0,strike,tau,rfr,coc,sig)
%
%where: value  is the value of the american call
%
%       s0     is the value of the underlying asset
%
%       strike is the option strike price (negative of strike for puts)
%
%       tau    is the time to expiration in years
%
%       rfr    is the risk free rate
%
%       coc    is the cost-to-carry
%              for non-div. paying stocks this is simply the rfr
%              for div paying stock this is rfr-divrate
%              for currency options this is rdom-rfor
%
%      sig   is the assumed volatility
%
%   Globals:
%
%      _cuts0  vector or scaler of cutoff points to override
%              analytically determined starting points
%
%Note: all inputs may be equal length vectors or scalers
%

global CUTS0 NITERS ;
ep=10^(-20);
rs0=rows(s0);
rs=rows(x); rt=rows(tau); rr=rows(rfr); rc=rows(coc); rst=rows(sig);
rowvec=([rs0;rs;rt;rr;rc;rst]);
n=max(rowvec);
if n>1 & min(rowvec)~= n;
    if rs0==1;
      s0=s0.*ones(n,1);
   end;
    if rs==1;
      x=x.*ones(n,1);
   end;
    if rt==1;
      tau=tau.*ones(n,1);
   end;
    if rr==1;
      rfr=rfr.*ones(n,1);
   end;
    if rc==1;
      coc=coc.*ones(n,1);
   end;
    if rst==1;
      sig=sig.*ones(n,1);
   end;
end;   if rows(CUTS0)==1;

      CUTS0=CUTS0(1).*ones(n,1);
end;
valvec=zeros(n,1);
cutvec=zeros(n,1);
ki=indexcat((x>0).*(coc<rfr),1);
if isempty(ki)==0;
[thisval,thiscut]=getcall(s0(ki),x(ki),tau(ki),rfr(ki),coc(ki),sig(ki),CUTS0(ki));

 valvec(ki)=thisval;
 cutvec(ki)=thiscut;
end;
ki=indexcat((x>0).*(coc>= rfr),1);
if isempty(ki)==0;
 [thisval]=bsval(sig(ki),x(ki),s0(ki)+ep,tau(ki),rfr(ki),coc(ki));
 valvec(ki)=thisval;
 cutvec(ki)=zeros(rows(ki),1);
end;
ki=indexcat(x<0,1);
if isempty(ki)==0;
                             [thisval,thiscut]=getput(s0(ki),x(ki),tau...
   (ki),rfr(ki),coc(ki),sig(ki),CUTS0(ki));

 valvec(ki)=thisval;
 cutvec(ki)=thiscut;
end;

% End of function 
function [value,curs]=getcall(s,x,t,r,b,sig,cut0)
global CUTS0 NITERS ;
 
ep=10^(-20);
m=2.*r./(sig.^2);
n=2.*b./(sig.^2);
k=1-exp(-r.*t);
q2=(-(n-1)+sqrt((n-1).^2+4.*m./k))./2;
curs=(cut0==-99).*s+(cut0~=-99).*cut0;
iter=1;
while iter <= NITERS;
   curs=curs+ep;
   eurval=bsval(sig,x,curs+ep,t,r,b);
   d1=(log(curs./x)+(b+0.5.*(sig.^2)).*t)./(sig.*sqrt(t));
   d2=d1-sig.*sqrt(t);
   d1=real(d1); d2=real(d2);
   lhs=curs-x;
   rhs=eurval+(1-exp((b-r).*t).*cdfn(d1)).*curs./(q2+ep);
   bi=exp((b-r).*t).*cdfn(d1).*(1-1./(q2+ep))+(1-exp((b-r).*t)...
    .*pdfn(d1)./(sig.*sqrt(t)))./(q2+ep);

   curs=(x+rhs-bi.*curs)./((1-bi)+ep);
   iter=iter+1;
end;
d1=(log((curs+ep)./x)+(b+0.5.*(sig.^2)).*t)./(sig.*sqrt(t));
a2=(curs./q2).*(1-exp((b-r).*t).*cdfn(d1));
value=(s<curs).*(bsval(sig,x,s+ep,t,r,b)+a2.*(s./(curs+ep)).^q2)+ ...
(s>=curs).*(s-x);


% End of function 
function [value,curs]=getput(s,x,t,r,b,sig,cut0)
global CUTS0 NITERS ;
 

ep=10^(-10);
x=abs(x);
m=2.*r./(sig.^2);
n=2.*b./(sig.^2);
k=1-exp(-r.*t);
q1=(-(n-1)-sqrt((n-1).^2+4.*m./k))./2;
b2=min([b;(0.6.*sig.*sqrt(t))]);
curs=(cut0==-99).*s+(cut0~=-99).*cut0;
iter=1;
while iter <= NITERS;
   eurval=bsval(sig,-x,curs+ep,t,r,b);
   d1=(log(((curs+ep)./(x+ep))+ep)+(b+0.5.*(sig.^2)).*t)./(sig....
   *sqrt(t));
   d2=d1-sig.*sqrt(t);
   d1=real(d1); d2=real(d2);
   lhs=x-curs;
   rhs=eurval-(1-exp((b-r).*t).*cdfn(-d1)).*curs./(q1+ep);
   bi=exp((b-r).*t).*cdfn(-d1).*(1./(q1+ep)-1)-(1+exp((b-r)....
      *t).*pdfn(-d1)./(sig.*sqrt(t)))./(q1+ep);

   curs=(x-rhs+bi.*curs)./((1+bi)+ep);
   iter=iter+1;
end;
d1=(log(((curs+ep)./(x+ep))+ep)+(b+0.5.*(sig.^2)).*t)./(sig.*sqrt...
   (t)+ep);
a1=-(curs./(q1+ep)).*(1-exp((b-r).*t).*cdfn(-real(d1)));
value=(s>curs).*(bsval(sig,-x,s+ep,t,r,b)+a1.*(abs(s./(curs+ep)))....
   ^q1)+(s<= curs).*(x-s);


% End of function