function ivvec= AMERIV(cp,k,s0,tau,rfr,coc,niters)
%  ivvec= AMERIV(cp,k,s0,tau,rfr,coc,niters)
%
%
% NOTE: YOU MUST ISSUE THE COMMAND: "amervald"
%       TO INITIALIZE GLOBAL VARIABLES BEFORE CALLING THIS PROCEDURE!!! 
%
%
%  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.
%
%

global CUTS0 NITERS ;
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=max(rowvec);
if n>1 & min(rowvec)~= n;
    if rs0==1;
      s0=s0.*ones(n,1);
   end;
    if rs==1;
      k=k.*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 rcp ~= n;
      cp=cp.*ones(n,1);
   end;
end;

if rows(CUTS0)~= n;

   CUTS0=CUTS0(1).*ones(n,1);
end;
% 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)~=-99;
 iv0=euriv(cp(vr),k(vr),s0(vr),tau(vr),rfr(vr),coc(vr),niters);
 iv0=min([iv0';ones(1,rows(iv0))]);
 iv0=max([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;
end;

% End of function 
function curiv=getivam(cp,k,s0,tau,rfr,coc,curiv,niters)
global CUTS0 NITERS ;
 
ep=10^(-10); ep2=10^(-200);
iter=1;
while iter <= 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;
end;

% End of function 
function [valvec,cutvec]=amerval(s0,x,tau,rfr,coc,sigma)
global CUTS0 NITERS ;
 

global CUTS0 NITERS ;
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=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;
      sigma=sigma.*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 missrv(ki,-99)~=-99;
[thisval,thiscut]=getcall(s0(ki),x(ki),tau(ki),rfr(ki),coc(ki),sigma...
   (ki),CUTS0(ki));

 valvec(ki)=thisval;
 cutvec(ki)=thiscut;
end;
ki=indexcat((x>0).*(coc>= rfr),1);
if missrv(ki,-99)~=-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);
end;
ki=indexcat(x<0,1);
if missrv(ki,-99)~=-99;
   [thisval,thiscut]=getput(s0(ki),x(ki),tau...
   (ki),rfr(ki),coc(ki),sigma(ki),CUTS0(ki));

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

% End of function 

function [value,curs]=getcall(s,x,t,r,b,sigma,cut0)
global CUTS0 NITERS ;
 
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~=-99).*cut0;
iter=1;
while iter <= NITERS;
   curs=curs+ep;
   eurval=bsval(sigma,x,curs+ep,t,r,b);
   d1=(log(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;
end;
d1=(log((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>= ...
   curs).*(s-x);


% End of function 

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

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=min([b;(0.6.*sigma.*sqrt(t))]);
curs=(cut0==-99).*s+(cut0~=-99).*cut0;
iter=1;
while iter <= NITERS;
   eurval=bsval(sigma,-x,curs+ep,t,r,b);
   d1=(log(((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;
end;
d1=(log(((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<= curs).*(x-s);


% End of function