/* proc impv,blacksh. Finds implied volatilities and bs formula */

        /* proc takes vectors. If s,k,r,t are scalars, same will be used for
                  each element of c */

proc(1) = impv(c,S,K,r,T);
    local rf,siglo,sigmid,sighi,clo,cmid,chi,i,j,iv,bound;

    iv = zeros(rows(c),1);
    if rows(c) > 1;
        if rows(s) == 1;
            s = s*ones(rows(c),1);
        endif;
        if rows(k) == 1;
            k = k*ones(rows(c),1);
        endif;
        if rows(r) == 1;
            r = r*ones(rows(c),1);
        endif;
        if rows(t) == 1;
            t = t*ones(rows(c),1);
        endif;
    endif;


    j = 1;
    do while j <= rows(c);
        bound =  maxc((s[j]-K[j]*exp(-r[j]*t[j]))|0);
        if c[j] <  bound;
            format /RE 10,5;
            "violating arbitrage bound c = " c[j] " bound " bound
                        "diff" c[j]-bound " returning zero ";
            format /RD 10,5;
            iv[j] = 0;
        elseif c[j] <= bound+1E-10;
            iv[j] = 0;
        elseif c[j] < 1e-10;
            iv[j] = 0;
        else;
            siglo = 1e-10;
            sighi = 5;
            i = 1;
            do while ((sighi-siglo) > 1e-7) and (i < 1000);
                clo = blacksh(s[j],k[j],t[j],r[j],siglo);
                chi = blacksh(s[j],k[j],t[j],r[j],sighi);
                sigmid = (siglo+sighi)/2;
                cmid = blacksh(s[j],k[j],t[j],r[j],sigmid);
                if (c[j] > chi) ;
                    "impv:error c above chi";
                elseif c[j] < clo;
                    "impv error. C below clo. c-clo" c[j]-clo "sigmas"
                       siglo~sighi;
                endif;
                if c[j] < cmid;
                    sighi = sigmid;
                else;
                    siglo = sigmid;
                endif;
                i = i+1;
             endo;
             if i >= 1000;
                "error: cant find sigma. diff";
                 format /RE 10,4; (sighi-siglo); format /RD 10,4;
             endif;
             iv[j] = sigmid;
        endif;
        j = j+1;
    endo;
    retp(iv);
endp;