/* EURIV.G

   procedure to calculate implied volatilities of a set of European option 
   quotes using a Newton Raphson type Method.

   See AMERIV for American option equivalent


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

 Where: iv  is a Nx1 vector or scaler of implied volatilities

        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 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 euriv(cp,k,s0,tau,rfr,coc,niters);
local rcp,rs0,rs,rt,rr,rc,rowvec,n;
local curiv,ep,ep2,iter,cfit,cfit2,vega;
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==1;
      cp=cp.*ones(n,1);
   endif;
endif;  
curiv=0.2.*ones(rows(cp),1); /* initial guess */
ep=10^(-10); ep2=10^(-200);
iter=1; 
do while iter le niters;
 {cfit}=bsval(curiv,k,s0,tau,rfr,coc);
 {cfit2}=bsval(curiv+ep,k,s0,tau,rfr,coc);
 vega=(cfit2-cfit)./ep;
 curiv=curiv+(cp-cfit)./(vega+ep2);
 iter=iter+1;
endo;
retp(curiv);
endp;