/* --------------------------------------------------------------------------*/
/*
                           BLACKSHL.PRG

THIS IS THE CURRENT VERSION

program to evaluate option pricing bounds. Call option, black-scholes
environment, no trade until expiration.

Produces figure 1 of the paper.

*/
/* --------------------------------------------------------------------------*/

library pgraph,optmum;
optset;
output file=blackshl.out reset;

/*-------------------------*/
/* switches to control run */
/*-------------------------*/

/* Calibrated to gross return but we need continuous time parameters.
/*
   The continuous time drift is defined as mu - 1/2 sig^2 so that the
   mean gross return is e^ mu T not e ^ mu T + 1/2 sig2 T */

   E(R) = E(e^r) = e^( Er + 1/2 sig2 ) = e^( mu T) ;
   sig2(R) = E(e^2r) - E(e^r)^2 =  e^( 2Er + 2sig2 ) - e^( E2r + sig2 )
           = e^(2Er +sig2)(e^sig2 - 1);
           = E(R)^2(e^sig2 -1)
   thus,

   sig2 = ln(1 + sig2(R)/E(R)^2)
   Er = lnE(R) - 1/2sig2;
   mu = lnE(R)
*/

/* assumptions */

justgrf = 0;          /* don't calculate, just play with graph */
r   = ln(1.05);       /* risk free rate */
T   = 1/4;            /* time to expiration (years) */
mu  = ln(1.13);       /* stock drift */
sig = (ln( 1 + (0.16^2)/(1.13^2) ))^0.5;
                      /* stock volatility:
                         this produces gross return sigma 16% */
K   = 100;            /* strike price */
facts = 2;            /* sharpe rat. bound = facts * stock+bond sharpe ratio. */

_mupwts = 0|0|0;
_mlowts = 0|0|0;      /* a global passed betwween procs.
                          weights of m on stock bond and option.
                          the positivity proc searches over these,
                          parameterized as lambda and delta.
                          the exact solution without positivity produces
                          starting values. */

/* derived parameters */

Rf = exp(r*T);                         /* discrete time gross risk free rate */
Ers = exp(mu*T) ;                      /*      "     gross stock exp. return */
sigR = Ers*(exp(T*sig^2)-1)^0.5;       /*      "     gross stock std. dev.   */
sr0 = (Ers-Rf)/sigR ;                  /* stock-bond sharpe ratio      */

srdes = facts*sr0;                     /* sharpe ratio target; h in paper */
biga = ((1+(srdes)^2)^0.5)/Rf;         /* second moment target E(m^2); A paper*/

/* echo assumptions to output */

format /RD 10,4;
"Assumptions"  ;
"  r" r;
" Rf" Rf;
"  T" T;
"sig" sig;
" mu" mu;
" ER" Ers;
" sR" sigR;
"  K" K;

"Stock and bond sharpe ratio" sr0;

/*---------------------------------------------------------*/
/* graph call option prices as function of stock prices */
/*---------------------------------------------------------*/

if not justgrf;     /*  allows you to just play with graph */


    sv = seqb(85,110,1);           /* stock grid */
    sv = sortc(sv|(K/Rf),1);         /* add K/Rf so lower arb will look ok*/
    cv = zeros(rows(sv),2);          /* place to store call prices*/
    cvpflg = cv;                     /* place to store flag that positive
                                        constraint is not binding */
    cvp = cv;                        /* place to store call prices w. pos.*/
    cbsv = zeros(rows(cv),1);        /* place to store bs price */

    arbnds = sv~maxc( ((sv-K/Rf)')|(0*sv') );
    j = 1;
    do while j <= rows(sv);
        S = sv[j];
        "S" S;
                /* no positivity answer*/
       {c,cpflg} = findcb(K,S,T,sig,r,Rf,mu,bigA);
        cv[j,.] = c';
        cvpflg[j,.] = cpflg';

       /* both binding. The trans proc keeps delta > 0 so
          you do not have to check arbitrage first. */

       _nopos = 0;

       /* translate from weights on xc, R, Rf, to lambda and delta */

       del = -1./_mlowts[3]; /* weight on option is -1/delta */
       lam0 = -_mlowts[1]*del; /* weight on stock */
       lam1 = -_mlowts[2]*del;  /* weight on bond */

       start = (lam0|lam1|del);
       _opmiter = 50;
       __output = 0;
       _minit = 1;                               /* do the lower bound */
       { lamdfin,clo,g,retcode } = optmum(&optnval,trans(start));
       if (retcode /= 0) and abs(clo)> 1E-3; /* wont converge when c=0 */
            "optmum retcode not 0, S,C=" S clo;
       endif;

       del = -1./_mupwts[3]; /* weight on option is -1/delta */
       lam0 = -_mupwts[1]*del; /* weight on stock */
       lam1 = -_mupwts[2]*del;  /* weight on bond */
       start = (lam0|lam1|del);

       _minit = 0;                                /* do the upper bound */

       { lamdfin,cup,g,retcode } = optmum(&optnval,trans(start));
       if (retcode /= 0) and abs(cup)> 1E-3; /* wont converge when c=0 */
           "optmum retcode not 0 S,C=" S cup;
       endif;
       cvp[j,.] = (cup|(-clo))';

       /* black-scholes for comparison */
       d1 = (ln(S/K) + (r + 1/2*sig^2)*T )/(sig*T^0.5);
       Cbsv[j] = S*cdfn(d1) - K/Rf*cdfn(d1-sig*T^0.5);

       j = j+1;
    endo;

endif;

"stock price, range where pos is slack, bounds.";
sv~cvpflg~cvp;



    graphjc;
    begwind;
    makewind(0,0,0,0,1);
    makewind(0,0,0,0,1);

    @title("Call price bounds as function of stock price, K=100");@
    xlabel("Stock price");
    ylabel("Call value");
    xtics(85,110,5,0);
    ytics(-2,12,2,0);
    px = sv;
    py = cvp[.,1]~cbsv~arbnds[.,2];
    _plwidth = 2;
    _pstype = 3;
    _psymsiz = 2;
    _pltype = 6|3|6;
    _plctrl = 2|0|0;

    _plegstr = "Good-deal bound\000Black-Scholes\000Arbitrage bound";
    _plegctl = 1|5|87|7;

    _pcolor = 15;

    xy(px,py);
    nextwind;

    py = cvp[.,2];
    _pstype = 3;
    _pltype = 6;
    _plctrl = 2;
    _plegctl = 0;
    xy(px,py);

    endwind;

    dos copy graphic.tkf blacksh.tkf;


output off;


/* ======================================================================== */
/*                                    PROCS                                 */
/* ======================================================================== */


/*-----------------------------------------------------------------*/
/* proc to calculate exact solutions with positivity assumed slack.

Returns:
   upper | lower call option bounds
   up | lo flag for positivity check
*/
/*------------------------------------------------------------------*/

proc(2) = findcb(K,S,T,sig,r,Rf,mu,bigA);
    local f,Nf,Nfp,Nfm,a0,a1,a2;
    local exx, exx1,exc_x,exc2,clo,cup,p,v,mupwtx,mlowtx,cupflg,cloflg;

    /* terms of lognormal distribution: see appendix */

    f = (ln(K/S) - mu*T)/(sig*T^0.5) - 1/2*sig*T^0.5;
    Nf  = 1 - cdfn(f);
    Nfp = 1 - cdfn(f+sig*T^0.5);
    Nfm = 1 - cdfn(f-sig*T^0.5);
    a0 = Nfp;
    a1 = exp(mu*T)*Nf;
    a2 = exp((2*mu+sig^2)*T)*Nfm;

    /* second moments for formula */

    Exx = ( exp((2*mu+sig^2)*T) ~ exp((r+mu)*T)  )|
          ( exp((r+mu)*T)       ~ exp(2*r*T)     );

    exx1 = inv(exx);
    Exc_x = (S*a2- K*a1) | ( exp(r*T)*(S*a1 - K*a0) );

    Exc2  = (S^2)*a2 - 2*K*S*a1 + K^(2)*a0;

    /* formula for bound */

    p = 1|1;
    clo = p'exx1*exc_x - ((biga^2 - p'Exx1*p)*(Exc2 - Exc_x'exx1*exc_x))^(0.5);
    cup = p'exx1*exc_x + ((biga^2 - p'Exx1*p)*(Exc2 - Exc_x'exx1*exc_x))^(0.5);

    /* check that m > 0, set flag cuppos accordingly. */
    /* also set  m weights on stock bond and option as
      starting values for positivity search. */

    v = ((biga^2 - p'Exx1*p)/(exc2 - exc_x'Exx1*exc_x) )^0.5;
    mupwtx = Exx1*(p-exc_x*v);
    "  mupwtx and v" mupwtx'~v;
    "  b*K/s+c*rf " mupwtx[2]*rf+mupwtx[1]*k/s;
    "  b/S+d      " mupwtx[1]/s+v;

    if (mupwtx[2] < 0) or ((mupwtx[2]*Rf + mupwtx[1]*K/S) < 0 )
            or ((mupwtx[1]/S + v) < 0 );
      cupflg = 0;
    " cup is not positive ";
    else;
    cupflg = 1;
    " cup is positive ";
    endif;

    mlowtx = Exx1*(p + v*exc_x);
    "  mlowtx and v" mlowtx'~v;
    "  b*K/s+c*rf " mlowtx[2]*rf+mlowtx[1]*k/s;
    "  b/S+d      " mlowtx[1]/s-v;
    if (mlowtx[2] < 0) or ((mlowtx[2]*Rf + mlowtx[1]*K/S) < 0 )
            or ((mlowtx[1]/S - v) < 0 );
      cloflg = 0;
    " clo is not positive ";
    else;
    cloflg = 1;
    " clo is positive ";
    endif;

    _mupwts = mupwtx|v;
    _mlowts = mlowtx|(-v);

retp(cup|clo,cupflg|cloflg);
endp;


/* ------------------------------------------------------------------------- */
/*                             SEQB.G                                        */
/*
    An alternative to seqa.  specify instead min, max, and a step size.  If
    max is not min + integer number * step size, the next smaller number is
    used.

    Usage:

    x = seqb(min,max,step);

*/
/* ------------------------------------------------------------------------- */

proc seqb(min,max,step);
    local nsteps, x;
     if max < min ; "seqb: asking for max < min"; endif;
    nsteps = trunc((max-min)/step+1.000000000001);
    x = seqa(min,step,nsteps);
    retp(x);
endp;



/* ====================================================================== */
/* procs for positivity search                                            */
/* ====================================================================== */

/* ---------------------------------------------------------------------- */
/* proc to evaluate

   int[a,b] (alf + bet*R)^2 f(R) dr

for R lognormal with mean mu*T and variance sig^2*T

use a = 9999 for + inf
*/
/*------------------------------------------------------------------------- */

proc(1) = integ(alf,bet,a,b);
    local  T2, botphi1, botphi2, botphi3, topphi1, topphi2, topphi3,
    astr, bstr, ans;

    T2 = T^0.5;

    if b <= a;            /* return 0 for b < a, not - int from b to a */
        retp(0);
    endif;

    if a <= 0;
        botphi1 = 0;
        botphi2 = 0;
        botphi3 = 0;
    elseif a >= 9999;
        botphi1 = 1;
        botphi2 = 1;
        botphi3 = 1;
    else;
        astr = (ln(a) - (mu - 1/2*sig^2)*T) / (sig*T2);
        botphi1 = cdfn(astr);
        botphi2 = cdfn(astr-sig*T2);
        botphi3 = cdfn(astr-2*sig*T2);
    endif;

    if b <= 0;
        topphi1 = 0;
        topphi2 = 0;
        topphi3 = 0;
    elseif b >= 9999;
        topphi1 = 1;
        topphi2 = 1;
        topphi3 = 1;
    else;
        bstr = (ln(b) - (mu - 1/2*sig^2)*T) / (sig*T2);
        topphi1 = cdfn(bstr);
        topphi2 = cdfn(bstr-sig*T2);
        topphi3 = cdfn(bstr-2*sig*T2);
    endif;

    ans = alf^2 * (topphi1 - botphi1) +
          2*alf*bet * exp((mu-1/2*sig^2)*T + sig^2*T/2) * (topphi2 - botphi2) +
          bet^2 * exp(2*((mu-1/2*sig^2)*T + sig^2*T)) * (topphi3 - botphi3) ;

    retp(ans);
endp;

/* ======================================================================= */
/* proc optnval

   returns option price as function of lambda, delta for
   minimizing. uses GLOBALs
    S = stock price
    biga = A^2 = E(m^2)
    Rf= gross risk free rate
    mu = stock drift
    sig = stock diffusion
    T = time to expiration
    K = strike price
    _minit = 1 if want minimum price
    _minit = 0 if want maximum price
    _nopos = 1 if want no positivity (this should agree with exact formulas abov
e)
*/

proc(1) = optnval(tlamdel);
    clearg lamdel, lam0l, lam1l, dell, int1, int2, optval, bot, top, alphl,betl;

    if cols(tlamdel) > 1;
       "rewrite optval to take cols(lamdell ) > 1 ";
    endif;

    lamdel = back(tlamdel);

    lam0l = lamdel[1];
    lam1l = lamdel[2];
    dell  = lamdel[3];

    if _minit and (dell <=0);
        "optnval: to minimize option price, dell must be > 0. dell=" dell;
    endif;
    if (not _minit) and (dell >= 0);
        "optnval: to maximize option price, dell must be < 0. dell=" dell;
    endif;


    if _nopos;            /* no positivity case */
        bot = K/S;
        top = 9999;
        alphl = lam0l*Rf -K;
        betl = (lam1l+S);
        int1 = 1/(dell^2)*integ( alphl, betl, bot, top );

        bot = 0;
        top = (K/S);
        alphl = lam0l*Rf;
        betl = lam1l;
        int2 = 1/(dell^2)*integ( alphl, betl, bot, top );

    endif;

    if not _nopos;

              /* positivity case */
      if (S+lam1l)/dell > 0;
        bot = K/S;
        top = (K-lam0l*Rf)/(lam1l+S);
        alphl = lam0l*Rf -K;
        betl = (lam1l+S);
        int1 = 1/(dell^2)*integ( alphl, betl, bot, top );
      elseif (S+lam1l)/dell < 0;
        bot = maxc((K/S)|((K-lam0l*Rf)/(lam1l+S)));
        top = 9999;
        alphl = lam0l*Rf -K;
        betl = (lam1l+S);
        int1 = 1/(dell^2)*integ( alphl, betl, bot, top );
      elseif (S+lam1l)/dell == 0;
        "figure out S+lam1l)/dell = 0 case";
      endif;

      if lam1l/dell > 0;
        bot = 0;
        top = minc((K/S)|(-lam0l*Rf/lam1l));
        alphl = lam0l*Rf;
        betl = lam1l;
        int2 = 1/(dell^2)*integ( alphl, betl, bot, top );
      elseif lam1l/dell < 0;
        bot = -lam0l*Rf/lam1l;
        top = K/S;
        alphl = lam0l*Rf;
        betl = lam1l;
        int2 = 1/(dell^2)*integ( alphl, betl, bot, top );
      elseif lam1l/dell == 0;
        "figure out lam1l/dell = 0 case";
      endif;
    endif;

    optval = -dell/2 * (int1 + int2) - lam0l - lam1l - dell/2*biga^2;

    if _minit;
        retp(-optval);     /* maximize wrt to lambda, if minimize optionvalue */
    else;
        retp(optval);
    endif;
endp;


/* from lamdel to transformed parameters */
proc(1) = trans(lamdel);
    local lam0l,lam1l,dell,tlam0,tlam1,tdel;
    tlam0 = lamdel[1];
    tlam1 = lamdel[2];
    tdel = abs(lamdel[3])^0.5;
    retp(tlam0|tlam1|tdel);
endp;
proc(1) = back(tlamdel);
    local lam0,lam1,del ;
    lam0 = tlamdel[1];
    lam1 = tlamdel[2];
    if _minit;
        del = tlamdel[3]^2;
    else;
        del = -tlamdel[3]^2;
    endif;
    retp(lam0|lam1|del);
endp;