/* BSVAL.G  Obtain European option, Black scholes values based on known
            (possibly smile dependent) implied volatilities

Format: value=bsval(iv,k,s0,tau,rfr,coc)

 Where: value is a Nx1 vector of Black Scholes European prices
        iv  is a Nx1 vector of (implied) volatilities
        k   is a Nx1 vector of strike prices (negative for puts)
        tau is a 1x1 scaler or Nx1 vector of time to exp. in yrs
        rfr is the continuously compounded annual risk free rate
        coc is the cost to carry

        Note: coc=rfr for a non-dividend paying stock
              coc=rfr-anndivyield for dividend paying stock
              coc=rdom-rfor for currency options
              coc=0 for options on futures


*/
proc bsval(iv,k,s0,tau,rfr,coc);
local n,scvec,rowind,mat,relrows,callp,strike,dep,beta,ind,sighat;
local d1,d2,cs,calldum,ep;
n=maxc(rows(iv)|rows(k)|rows(s0)|rows(tau)|rows(rfr)|rows(coc));
if rows(s0)==1;
   s0=s0.*ones(n,1);
endif;
if rows(iv)==1;
   iv=iv.*ones(n,1);
endif;
if rows(tau)==1;
   tau=tau.*ones(n,1);
endif;
if rows(rfr)==1;
   rfr=rfr.*ones(n,1);
endif;
 if rows(coc)==1;
   coc=coc.*ones(n,1);
 endif;
 calldum=k.>0;
 ep=10^(-10);
 d1=(ln(abs((s0+ep)./(abs(k)+ep)+ep))+(coc+(iv.^2)./2).*tau)./
                   (iv.*sqrt(tau)+ep);
 d2=d1-iv.*sqrt(tau);
 cs=calldum.*(s0.*exp((coc-rfr).*tau).*cdfn(real(d1))-
    (k).*(exp(-rfr.*tau)).*cdfn(real(d2)));
 cs=cs+(1-calldum).*((abs(k)).*exp(-rfr.*tau).*cdfn(-real(d2))-
                   s0.*exp((coc-rfr).*tau).*cdfn(-real(d1)));
retp(cs);
endp;