/* calculates real option example */

library pgraph;
output file=reaopt.out reset;
convt = 0;

Rfl = 1.05;
ERl = 1.13;
sigRl = 0.16;

rhov = 0.9; /* correlation */
atoh = 2;      /* A/h */

V = seqb(80,110,1);
K = 100;
T  = 1/4;

/* find implied continuous time parameters */

/* ER = exp Er + 1/2 sigr^2
   E(R^2)= exp  2Er +2sigr^2

   lnER  = Er + 1/2 sigr^2
   lnER2 = 2Er + 2sigr^2

   sigr^2 = lnER2 - 2lnER
   er = lnER - 1/2 sigr^2 */

sigr = (ln(ERl^2+sigRl^2) - 2*ln(ERl))^0.5;
Er = ln(ERl) - 1/2*sigr^2;
rf = ln(Rfl);

"Continuous time parameters of stock process" ;
"   Er " er;
"sig r " sigr;
"   rf " rf;

ir = 1;
do while ir <= rows(rhov);
rho = rhov[ir];
"working on rho" rho;

/* Sigvz/sigv = rho. Let sigV = sigS -- same as market volatility. Then
   sigvz = sigv*rho;
   sigvw = (sigv^2-sigvz^2)^0.5
*/
sigv = sigr;
sigvz = sigr*rho;
sigvwh = (sigv^2-sigvz^2)^0.5;

"sigv,sigvz,sigvw";
sigv sigvz sigvw;

/* now GO */


a = 1;
d = ((Er-rf)/sigr - (er-rf)/sigr*(rho-a*(atoh^2-1)^0.5*(1-rho^2)) )*sigv;

"dup" d;

Cup = V.*exp(d*T).*cdfnc((ln(K./V)-(d+rf)*T)./(sigv*T^0.5)-1/2*sigv*T^0.5)-
    K.*exp(-rf*T).*cdfnc((ln(K/V)-(d+rf)*T)./(sigv*T^0.5)+1/2*sigv*T^0.5);

a = -1;
d = ((Er-rf)/sigr - (er-rf)/sigr*(rho-a*(atoh^2-1)^0.5*(1-rho^2)) )*sigv;

"dlo" d;

Clo = V.*exp(d*T).*cdfnc((ln(K./V)-(d+rf)*T)./(sigv*T^0.5)-1/2*sigv*T^0.5)-
    K.*exp(-rf*T).*cdfnc((ln(K/V)-(d+rf)*T)./(sigv*T^0.5)+1/2*sigv*T^0.5);

"upper" cup[1:5]';
"lower" clo[1:5]';
if ir == 1;
    cupv = cup;
    clov = clo;
else;
    cupv = cupv~cup;
    clov = clov~clo;
endif;

ir = ir+1;
endo;

"Upper v";
cupv[1:5,.]';
"lower v";
clov[1:5,.]';

Cbs = V.*cdfnc((ln(K./V)-rf*T)./(sigv*T^0.5)-1/2*sigv*T^0.5)-
    K.*exp(-rf*T).*cdfnc((ln(K/V)-rf*T)./(sigv*T^0.5)+1/2*sigv*T^0.5);


graphjc;
@title("Call option price bounds with imperfect replication");@
xtics(85,110,5,0);
@ytics(0,20,5,0);@
_pcolor = 15;
_pltype = 6|6|6|6|3;
_plctrl = 1;
_plegstr = "Upper\000Lower\000Black-Scholes";
_plegctl = 1;
xlabel("V");
ylabel("Option value");
xy(V,cupv~clov~cbs);
if convt;
    dos pqgrun graphic.tkf -cf=realop1.pic -c=2;
endif;

graphjc;
@title("Call option price bounds / Black Scholes price");@
_pcolor = 15;
_pltype = 6|6|3;
_plwidth = 2;
_plctrl = 0;
_plegstr = "Upper\000Upper\000Lower\000Lower\000Black-Scholes";
_plegctl = 0;
xtics(85,110,5,0);
ytics(-20,30,10,2);
xlabel("V");
ylabel("Bound / Black-Scholes (%)");
xy(V,100*(( (cupv./cbs)~(clov./cbs)~(ones(rows(v),1)) ) -1));
dos copy graphic.tkf realop.tkf;

/* graph in terms of implied volatility */
/* watch out -- there is no reason this should not go below arb bound */

if 0;
lowimpv = (clo);
hiimpv = (cup);
bsimpv = (cbs);

i = 1;
do while i <= rows(clo);
    lowimpv[i] = impv(clo[i],v[i],K,rf,1);
    hiimpv[i]  = impv(cup[i],v[i],K,rf,1);
    bsimpv[i] = impv(cbs[i],v[i],K,rf,1);
    i = i+1;
endo;

graphjc;
title("Implied volatility bounds with imperfect replication");
xtics(85,110,5,0);
_pcolor = 15;
_pltype = 6|6|3;
_plegstr = "Upper\000Lower\000Black-Scholes";
_plegctl = 1|4|90|0.05;
xlabel("V");
ylabel("B-S implied volatility");
xy(V,hiimpv~lowimpv~bsimpv);

endif;

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

    Rf = exp(-r*T);
    if c <= maxc((s-K/Rf)|0);
        "Violating lower arbitrage bound, report 0. c=" C "s-k/rf" s-k/Rf;
        retp(0);
    else;
        siglo = 0.000001;
        sighi = 5;
        i = 1;
        do while ((sighi-siglo) > 1e-7) and (i < 1000);
            clo = blacksh(s,k,t,r,siglo);
            chi = blacksh(s,k,t,r,sighi);
            sigmid = (siglo+sighi)/2;
            cmid = blacksh(s,k,t,r,sigmid);
            if c <= 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;
         retp(sigmid);
    endif;
endp;


proc(1) = blacksh(S,K,T,r,sig);
    local d1,d2,c;
    d1 = (ln(S/K)+(r+sig^2/2)*T)/(sig*T^0.5);
    d2 = d1-sig*T^0.5;
    c = S*cdfn(d1)-K*exp(-r*T)*cdfn(d2);
    retp(c);
endp;



output off;