/* AMERIV.G

  procedure to calculate implied volatilities of a set of AMERICAN option
  quotes using a Newton Raphson type Method, and an analytical approx. to
  American option prices as in Barone-Adesi and Whaley

Format: iv=AMERIV(cp,k,s0,tau,rfr,coc,niters)

 Where: iv  is a Nx1 vector or scaler of implied volatilities
            if the code -99 is returned, this indicates the option should be
            trading at its exercisable value.  May also want to error check
            iv afterwards and eliminate any unreasonable values.

        cp  is a Nx1 vector or scaler of call prices (or put prices)

        k   is a Nx1 vector or scaler of strike prices (negative of strike
            for puts)

        s0  is a Nx1 vector or scaler for the current level of the
            underlying asset

        tau is a Nx1 vector or sclaer of time to exp. in yrs

        rfr is a Nx1 vector or scaler for the continuously compounded
            annualized risk free rate
            Note: if you have an effective annual yield (annyield) simply
            pass ln(1+annyield) to the procedure.

        coc is the annualized rate for cost-to-carry
            for non-div. paying stock this is simply the rfr
            for div paying stock this is rfr-divrate
            for currency options this is rdom-rfor

     niters is a scaler indicating the number Newton-Raphson iterations
            to conduct (3 or 4 yields reasonable results)

      note: You may want to order the command disable before calling this
            function, for particularly wild data. Be sure to issue the
            command enable afterwards, to reactivate error checking.

*/
proc AMERIV(cp,k,s0,tau,rfr,coc,niters);
local rcp,rs0,rs,rt,rr,rc,rowvec,n;
local iv0,cfit0,cuts0,ivvec,intrdum,ks;
local valdum,vr,ivvec2;
#include ameriv.dec;
rcp=rows(cp); rs0=rows(s0);
rs=rows(k); rt=rows(tau); rr=rows(rfr); rc=rows(coc);
rowvec=(rcp|rs0|rs|rt|rr|rc);
n=maxc(rowvec);
if n>1 and minc(rowvec) ne n;
   if rs0==1;
      s0=s0.*ones(n,1);
   endif;
   if rs==1;
      k=k.*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 rcp ne n;
      cp=cp.*ones(n,1);
   endif;
endif;  
if rows(_cuts0) ne n;
   _cuts0=_cuts0[1].*ones(n,1);
endif;
/* obtain European iv's first and determine initial cutoffs from
   those values */
valdum=(k.<0).*(cp.>(abs(k)-s0))+(k.>0).*(cp.>(s0-abs(k)));
vr=indexcat(valdum,0.99|1.01);
ivvec=-99.*ones(n,1);
if missrv(vr,-99) ne -99;
 iv0=euriv(cp[vr],k[vr],s0[vr],tau[vr],rfr[vr],coc[vr],niters);
 iv0=minc(iv0'|ones(1,rows(iv0)));
 iv0=maxc(iv0'|(0.01.*ones(1,rows(iv0))));
 {ivvec2}=getivam(cp[vr],k[vr],s0[vr],tau[vr],rfr[vr],coc[vr],iv0,niters);
 ivvec[vr]=ivvec2;
endif;
retp(ivvec);
endp;

proc(1)=getivam(cp,k,s0,tau,rfr,coc,curiv,niters);
local ep,ep2,iter,cfit,cuts,cfit2,vega;
ep=10^(-10); ep2=10^(-200);
iter=1;
do while iter le niters;
 {cfit,cuts}=amerval(s0,k,tau,rfr,coc,curiv);
 {cfit2,cuts}=amerval(s0,k,tau,rfr,coc,curiv+ep);
 _cuts0=cuts;
 vega=(cfit2-cfit)./ep;
 curiv=curiv+(cp-cfit)./(vega+ep2);
 iter=iter+1;
endo;
retp(curiv);
endp;

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;