/* 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,sigma)

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

      sigma   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
*/
proc(2)=amerval(s0,x,tau,rfr,coc,sigma);
local rs0,rs,rt,rr,rc,rst,rowvec,n,valvec,cutvec,i,thisval,thiscut;
local ki,ep;
#include amerval.dec;
ep=10^(-20);
rs0=rows(s0);
rs=rows(x); rt=rows(tau); rr=rows(rfr); rc=rows(coc); rst=rows(sigma);
rowvec=(rs0|rs|rt|rr|rc|rst);
n=maxc(rowvec);
if n>1 and minc(rowvec) ne n;
   if rs0==1;
      s0=s0.*ones(n,1);
   endif;
   if rs==1;
      x=x.*ones(n,1);
   endif;
   if rt==1;
      tau=tau.*ones(n,1);
   endif;
   if rr==1;
      rfr=rfr.*ones(n,1);
   endif;
   if rc==1;
      coc=coc.*ones(n,1);
   endif;
   if rst==1;
      sigma=sigma.*ones(n,1);
   endif;
endif;  
if rows(_cuts0)==1;
      _cuts0=_cuts0[1].*ones(n,1);
endif;
valvec=zeros(n,1);
cutvec=zeros(n,1);
ki=indexcat((x.>0).*(coc.<rfr),1);
if missrv(ki,-99) ne -99;
{thisval,thiscut}=getcall(s0[ki],x[ki],tau[ki],rfr[ki],coc[ki],sigma[ki],
                            _cuts0[ki]);
 valvec[ki]=thisval;
 cutvec[ki]=thiscut;
endif;
ki=indexcat((x.>0).*(coc.ge rfr),1);
if missrv(ki,-99) ne -99;
 {thisval}=bsval(sigma[ki],x[ki],s0[ki]+ep,tau[ki],rfr[ki],coc[ki]);
 valvec[ki]=thisval;
 cutvec[ki]=zeros(rows(ki),1);
endif;
ki=indexcat(x.<0,1);
if missrv(ki,-99) ne -99;
 {thisval,thiscut}=getput(s0[ki],x[ki],tau[ki],rfr[ki],coc[ki],sigma[ki],
                            _cuts0[ki]);
 valvec[ki]=thisval;
 cutvec[ki]=thiscut;
endif;
retp(valvec,cutvec);
endp;

proc(2)=getcall(s,x,t,r,b,sigma,cut0);
local ep,m,n,k,q2,curs,iter,eurval,d1,d2,lhs,rhs,bi,a2,value;
ep=10^(-20);
m=2.*r./(sigma.^2);
n=2.*b./(sigma.^2);
k=1-exp(-r.*t);
q2=(-(n-1)+sqrt((n-1).^2+4.*m./k))./2;
curs=(cut0.==-99).*s+(cut0.ne -99).*cut0;
iter=1;
do while iter le _niters;
   curs=curs+ep;
   eurval=bsval(sigma,x,curs+ep,t,r,b);
   d1=(ln(curs./x)+(b+0.5.*(sigma.^2)).*t)./(sigma.*sqrt(t));
   d2=d1-sigma.*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)./(sigma.*sqrt(t)))./(q2+ep);
   curs=(x+rhs-bi.*curs)./((1-bi)+ep);
   iter=iter+1;
endo;
d1=(ln((curs+ep)./x)+(b+0.5.*(sigma.^2)).*t)./(sigma.*sqrt(t));
a2=(curs./q2).*(1-exp((b-r).*t).*cdfn(d1));
value=(s.<curs).*(bsval(sigma,x,s+ep,t,r,b)+a2.*(s./(curs+ep)).^q2)+
      (s.ge curs).*(s-x);
retp(value,curs);
endp;

proc(2)=getput(s,x,t,r,b,sigma,cut0);
local q1,b2,curs,iter,eurval,d1,d2,lhs,rhs,bi,a1,value;
local ep,m,n,k;
ep=10^(-10);
x=abs(x);
m=2.*r./(sigma.^2);
n=2.*b./(sigma.^2);
k=1-exp(-r.*t);
q1=(-(n-1)-sqrt((n-1).^2+4.*m./k))./2;
b2=minc(b|(0.6.*sigma.*sqrt(t)));
curs=(cut0.==-99).*s+(cut0.ne -99).*cut0;
iter=1;
do while iter le _niters;
   eurval=bsval(sigma,-x,curs+ep,t,r,b);
   d1=(ln(((curs+ep)./(x+ep))+ep)+(b+0.5.*(sigma.^2)).*t)./(sigma.*sqrt(t));
   d2=d1-sigma.*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)./(sigma.*sqrt(t)))./(q1+ep);
   curs=(x-rhs+bi.*curs)./((1+bi)+ep);
   iter=iter+1;
endo;
d1=(ln(((curs+ep)./(x+ep))+ep)+(b+0.5.*(sigma.^2)).*t)./(sigma.*sqrt(t)+ep);
a1=-(curs./(q1+ep)).*(1-exp((b-r).*t).*cdfn(-real(d1)));
value=(s.>curs).*(bsval(sigma,-x,s+ep,t,r,b)+a1.*(abs(s./(curs+ep))).^q1)+
      (s.le curs).*(x-s);
retp(value,curs);
endp;