/* MULTIOPT.PRG program to evaluate option pricing bounds and deltas. Call option, black-scholes environment, no trade until expiration, multiple options. THIS IS THE FINAL VERSION It is copied from multiopt.prg and deltas.prg and updated. Note: you generally can ignore optmum not converging due to maxiterations. This occurs typcially when very near or at an arbitrage bound, since the program does not separately check E(m2) for arbitrage bounds. */ /* --------------------------------------------------------------------------*/ library pgraph,optmum; optset; output file=multiopt.out reset; clearg _a, _d, _g, _MG, _alpha; /* globals created by procs */ /*-------------------------*/ /* switches to control run */ /*-------------------------*/ /* assumptions */ justgrf = 1; /* don't calculate, just do graph */ testnopos = 0; /* don't use pos. Used to check that tigher deltas than arbitrage really is true */ S = 100; XM = {95,100,105}; /* strike prices of traded options for which prices are known. You have to introduce them as an increasing sequence. */ Xx = seqb(XM[1]-5,XM[rows(XM)]+5,.5); /* strikes at which to compute */ dx = 0.0001; /* delta x for delta derivatives wrt strike */ 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. */ /* 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 ancd bond sharpe ratio" sr0; /* ----------------------------------------------------- */ /* graph call option prices as function of strike prices */ /* Multiple option Case */ /* ----------------------------------------------------- */ if (not justgrf); Nopt = rows(XM) + 1; /* no. of options, including one we price */ /* assign black-scholes prices to the traded options */ /* program does not require this asumption */ CM = zeros(rows(XM),1); /* place to put prices of known options */ delBS = zeros(rows(XM),1); /* deltas of known options */ m = 1; do while m < Nopt; d1a = (ln(S/XM[m]) + (r + 1/2*sig^2)*T )/(sig*T^0.5); CM[m] = S*cdfn(d1a) - XM[m]/Rf*cdfn(d1a-sig*T^0.5); delBS[m] = cdfn(d1a); m = m+1; endo; /* add strikes at delta x for later derivatives */ xx = sortc(xx|(xm+dx)|(xm-dx),1); /* Eliminate strikes with known*/ /* prices from calculation (They get put back later)*/ i = 1; do while i <= rows(Xx); if minc(abs(Xx[i]-Xm)) < 1E-4 ; if i == 1; Xx = Xx[2:rows(Xx)]; elseif i == rows(Xx); Xx = Xx[1:rows(Xx)-1]; else; Xx = Xx[1:i-1]|Xx[i+1:rows(Xx)]; endif; endif; i = i+1; endo; n = 1; do while n <= rows(Xx); "K = " xx[n]; /* produce globals for proc. NewX = strikes of known & unknown, in order. NewC = price of knonwn & 0 Inx = index of known */ NewXC = sortc((XM|xx[n])~(CM|0),1); NewX = newXC[.,1]; NewC = newXC[.,2]; Inx = indexcat(Newc,0); if ismiss(inx); "inx contains missings, fix indexcat call"; endif; /* non positivity */ Cp= findmcb(S,T,sig,Rf,mu,(Xx[n]|XM),Nopt,(-1|CM),srdes); "Upper lower non pos" cp'; /* store results */ if n == 1; CCup=Cp[1]; CCdow=Cp[2]; elseif n > 1; CCup=(CCup|Cp[1]); CCdow=(CCdow|Cp[2]); endif; /* positivity */ _nopos = 0; del = -S./_alpha[1,.]'; lam0 = S*_alpha[Nopt+2,.]'./(_alpha[1,.]'.*Rf); lam1 = S*_alpha[Nopt+1,.]'./_alpha[1,.]'; starthi = (lam0[1]|lam1[1]|del[1]| _alpha[2:Nopt,1]./_alpha[1,1]); startlo = (lam0[2]|lam1[2]|del[2]| _alpha[2:Nopt,2]./_alpha[1,2]); _opmiter = 100; __output = 0; /* Computing the upper bound */ _minit = 0; { lamdfin,ch,g,retcode } = optmum(&optmopt,trans(starthi)); if (retcode /= 0) and abs(ch)>1E-3; /* wont converge when c=0 */ " Up pos optmum retcode not 0. Retcode, C=" retcode ch; " final lam del" back(lamdfin)'; endif; "Upper pos" ch; /* store results */ if n == 1; CPup = ch; elseif n > 1; CPup=(CPup|Ch); endif; /* Computing the lower bound */ _minit = 1; _nopos = 0; { lamdfin,cl,g,retcode } = optmum(&optmopt,trans(startlo)); if (retcode /= 0) and abs(cl)>1E-3; /* wont converge when c=0 */ " Low pos optmum retcode not 0. Code C=" retcode cl; " final lam del" back(lamdfin)'; endif; "Lower pos" -cl; /* store result */ if n == 1; CPdow = -Cl; Xv = Xx[n]; elseif n > 1; CPdow = (CPdow|-Cl); Xv = (Xv|Xx[n]); endif; n = n+1; endo; /* add known strikes and prices */ result = sortc((Xv|XM)~(CCup|CM)~(CCdow|CM)~(CPup|CM)~(CPdow|CM),1); xv = result[.,1]; CCup = result[.,2]; CCdow = result[.,3]; CPup = result[.,4]; CPdow = result[.,5]; if testnopos; clo = ccdow; cup = ccup; endif; clo = maxc((ccdow~cpdow)'); /* due to rounding error, a few violated this*/ cup = minc((ccup~cpup)'); /*no violations, but in case you change assn*/ endif; /* temporary ends if grafonly */ /* calculate arbitrage bounds by concavity over strikes */ /* this version is not very elegant and is hard-wired to three strikes */ if rows(Xm) /= 3; "arbitrage bound calculation is hard wired to three strikes. Fix it"; endif; /* to fit a line through x1,y1 and x2,y2 (y-y1)/(x-x1) = (y2-y1)/(x2-x1) y = y1 + (y2-y1)/(x2-x1)*(x-x1) */ xarblo = xv; /* use same strikes for arb bounds to start */ xarblo = sortc(xarblo|seqb(102.9,103.5,0.02),1); /* art near bend */ /* line through first and second strike */ abnd1 = CM[1] + (CM[2]-CM[1])/(XM[2]-XM[1]) * (xarblo-Xm[1]); abnd2 = ones(rows(xarblo),1)*cm[3]; /* horizontal line at third x */ abnd1 = maxc((abnd1~abnd2)'); /* bend at lower right */ abnd1 = maxc((abnd1~(S-xarblo/Rf))'); /* make sure result is > */ abnd1 = maxc((abnd1~(0*xarblo))'); /* regular (S,B) bounds */ /* line through second and third strike */ abnd2 = CM[2] + (CM[3]-CM[2])/(XM[3]-XM[2]) * (xarblo-Xm[2]); abnd2 = maxc((abnd2~(S-xarblo/Rf))'); /* regular (S,B) bounds */ abnd2 = maxc((abnd2~(0*xarblo))'); arblo = selif(abnd1,xarblo.<=xm[1])| /* 2-3 line here */ selif(minc((abnd1~abnd2)'),xarblo .> xm[1]); abnd1 = selif(abnd1, xarblo .> xm[1]); /* goes to inf left of x[1]*/ abnd2 = selif(abnd2, xarblo .> xm[1]); arbup = cm[1]|maxc((abnd1~abnd2)'); xarbup = selif(xarblo, xarblo .>= xm[1]); xarbup = xm[1]|xarbup; /* add point so line goes up */ arbup = 50|arbup; /* calculate delta bounds from derivatives of bounds near known */ /* keep forward and backward derivatives */ /* xv, clo and cup are conformable */ indx = indexcat(xv,xm[1]+(-dx|dx)/3); /* find index of lowest known */ dcdxu1 = (cup[indx]-cup[indx-1])/dx; dcdxu2 = (cup[indx+1]-cup[indx])/dx; dcdxl1 = (clo[indx]-clo[indx-1])/dx; dcdxl2 = (clo[indx+1]-clo[indx])/dx; indx = indexcat(xarblo,xm[1]+(-dx|dx)/10); /* find index of xm[1],lo arb*/ dcdxla1 = (arblo[indx]-arblo[indx-1])/dx; /* lower arb, backwd */ dcdxla2 = (arblo[indx+1]-arblo[indx])/dx; /* lower arb forward */ indx = 2; /* 50, xm[1]-0.01,this one*/ dcdxua1 = -9999;/* infinity */ /* upper arb back */ dcdxua2 = (arbup[indx+1]-arbup[indx])/dx; /* upper arb fwd */ i = 2; do while i <= rows(xm); indx = indexcat(xv,xm[i]+(-dx|dx)/3); /* find index of XM[i]*/ dcdxu1 = dcdxu1|((cup[indx]-cup[indx-1])/dx); dcdxu2 = dcdxu2|((cup[indx+1]-cup[indx])/dx); dcdxl1 = dcdxl1|((clo[indx]-clo[indx-1])/dx); dcdxl2 = dcdxl2|((clo[indx+1]-clo[indx])/dx); indx = indexcat(xarblo,xm[i]+(-dx|dx)/10); /* find index of xm[1] in low arb x */ dcdxla1 = dcdxla1|((arblo[indx]-arblo[indx-1])/dx); /* lo arb, back */ dcdxla2 = dcdxla2|((arblo[indx+1]-arblo[indx])/dx); /* lo arb fowd */ indx = indexcat(xarbup,xm[i]+(-dx|dx)/10); dcdxua1 = dcdxua1|((arbup[indx]-arbup[indx-1])/dx); /* up arb back */ dcdxua2 = dcdxua2|((arbup[indx+1]-arbup[indx])/dx); /* up arb fwd */ i = i+1; endo; /* delta = dc/ds = c/s - dc/dk*k/s */ delfup = cm./S - dcdxu2.*xm./s; /* forward, upper */ delbup = cm./S - dcdxu1.*xm./s; /* backward upper */ delflo = cm./S - dcdxl2.*xm./s; /* forward lower */ delblo = cm./S - dcdxl1.*xm./s; /* backward lower */ delafup = cm./S - dcdxua2.*xm./s; /* forward, upper */ delabup = cm./S - dcdxua1.*xm./s; /* backward upper */ delaflo = cm./S - dcdxla2.*xm./s; /* forward lower */ delablo = cm./S - dcdxla1.*xm./s; /* backward lower */ "deltas: note non-infinity ones are same backward or forward to check"; delafup~delabup~delaflo~delablo; @endif; /* ends just graph if */@ if 1; /* use 0 to skip graph */ graphjc; begwind; makewind(0,0,0,0,1); makewind(0,0,0,0,1); _pcolor = 15; @title("Option Prices as function of strikes. Traded Options, S=100");@ xlabel("Strike Price"); ylabel("Call Value"); px = fill(xv,xM); py = fill(cup,cM); px = fill(px,xarblo); py = fill(py,arblo); _plegstr = "Good-deal bound\000Known options\000Arbitrage bound"; _plegctl = 1|5|100|6; ytics(0,10,2,0); xtics(90,110,5,0); _plwidth = 2; _pstype = 3|9|1; _pltype = 6|3|3; _plctrl = 0|-1|0; xy(px,py); nextwind; _pstype = 3|1; _pltype = 6|3; _plctrl = 0|0; px = fill(Xv,xarbup); py = fill(clo,arbup); _plegctl = 0; xy(px,py); endwind; dos copy graphic.tkf multiopt.tkf; endif; if 1; graphjc; @title("Option Delta Bounds:Hedging with other Options S=100");@ xlabel("Strike Price"); ylabel("Call Delta Bounds"); px = (XM)~(xm)~(xm)~(xm)~xm; /* left bounds and bs */ py = delbup~delblo~delabup~delablo~delBS; _pstype = 13|10|1|1|9; _psymsiz =4.5; _pltype = 6; _plctrl = -1|-1|-1|-1|-1; m = 1; _pline = (1~6~(XM[m])~delablo[m]~(XM[m])~delabup[m]~1~15~2) ; /* make sure top of arrow is at ytics upper setting */ _Parrow = (xm[m]~delablo[m]~xm[m]~1.2~1~0.13~21~15~1~6~2); /* xs ~ ys ~ xe ~ ye ~ hdlw ~ hdsz ~ hdtyp ~ co(10) ~ units(1) ~ lintype ~ thick; note: type e.g. 10 first number: 0:solid 1:empty 2:open 3: closed controls head second number: 0:none 1:@end 2:@start */ m = 2; do while m < Nopt; _pline = _pline| (1~6~(XM[m])~delablo[m]~(XM[m])~delabup[m]~1~15~2) ; m = m+1; endo; _plegstr = "Upper\000Lower\000\Upper Arbitrage"$+ "\000Lower Arbitrage\000Black-Scholes"; _plegctl = 1|5|105|.6; xtics(90,110,5,0); ytics(0,1.2,0.2,0); _pcolor = 15; xy(px,py); dos copy graphic.tkf deltas.tkf; endif; /* _pnum = 1|2; _pdate = 0; _pframe = 0|1; xtics(70,130,10,0); _plegstr = "Upper 2\000Upper 10\000Lower 2\000Lower 10\000Black-Scholes"; _plegctl = 1|5|75|-0.4; xy(px,py); _protate = 1; begwind; makewind(9,4,0,2.,1); _ptitlht=0.25; _paxht=0.25; _pnumht=0.2; _psymsiz=4; xy(px,py); endwind; */ output off; /* ======================================================================== */ /* PROCS */ /* ======================================================================== */ /* ************* FINDCB ***************/ /*-----------------------------------------------------------------*/ /* ONE OPTION CASE */ /* proc to return either sharpe ratio or call price, no positivity . feed -1 for the one you want returned as a function of the other. One of C or srdes MUST be = -1. Returns: 1) if srdes < 0 and c > 0 returns sharpe ratio of stock, bond and option all together | sharpe ratio of option and bond alone. 2) if srdes > 0 and c < 0 returns upper | lower call option bounds */ /*------------------------------------------------------------------*/ proc(1) = findcb(X,S,T,sig,Rf,mu,C,srdes); local f,Nf,Nfp,Nfm,A,B1,B2,B3,B4,B5,B6,B,v,SR,W,T1,T2,T3,disc,cs,erc, sigrc,shrc,a0,a1,a2; f = (ln(X/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; /* matrix A is coefficients on a, d, g, in bond, stock and option formula */ A = zeros(3,3); /* bond formula */ A[1,1] = a1 - X/S*a0; A[1,2] = exp(mu*T); A[1,3] = 1; /* stock formula */ A[2,1] = a2 - X/S*a1; A[2,2] = exp((2*mu+sig^2)*T); A[2,3] = exp(mu*T); /* option formula */ A[3,1] = a2 - 2*X/S*a1 + (X/S)^2*a0; A[3,2] = a2 - (X/S)*a1; A[3,3] = a1 - (X/S)*a0; /* Terms on quadratic form for Sharpe ratio */ B1 = a2 - 2*(X/S)*a1 + (X/S)^2*a0; /* a2 */ B2 = a2 - (X/S)*a1; /* ad */ B3 = a1 - X/S*a0; /* ag */ B4 = exp((2*mu+sig^2)*T); /* d2 */ B5 = exp(mu*T); /* dg */ B6 = 1; /* g2 */ B = (B1~B2~B3)| (B2~B4~B5)| (B3~B5~B6); /* now ready to find sharpe ratio given call price */ if (c >= 0) and (srdes < 0); /* if want sr as function of c */ v = ((C/S)|1|1/Rf); SR=v'*inv(B)*v; SR=RF^2*SR -1; SR=SR^(0.5); Erc = S/C*(a1 - X/S*a0); Sigrc = (S/C)^2 * (a2 - a1^2 + (X/S*a0 - 2*a1)*X/S*(1-Nfp)); if sigrc >= 0; sigrc = sigrc^0.5; else; "formula for variance of call option return not positive"; sigrc = (-sigrc)^0.5; endif; shrc = abs(Erc-Rf)/sigrc; retp(sr|shrc); /* call price for given sharpe ratio */ elseif (c < 0) and (srdes >= 0); W=inv(B); T1 = w[1,1]*Rf^2; T2 = 2*(W[1,2]*Rf^2 + W[1,3]*Rf); T3 = W[2,2]*Rf^2 + 2*W[2,3]*Rf + W[3,3] - 1 - srdes^2; disc = T2^2 - 4*T1*T3; if disc >= 0 ; cs = (-T2 + disc^0.5)/(2*T1); cs = cs| ((-T2 - disc^0.5)/(2*T1)); c = s*cs; clearg adg,v1,v2; v1 = (1/Rf)|1|(C[1]/S); adg = inv(A)*v1; _a = adg[1]; _d = adg[2]; _g = adg[3]; v2 = (1/Rf)|1|(C[2]/S); adg = inv(A)*v2; _a = _a|adg[1]; _d = _d|adg[2]; _g = _g|adg[3]; retp(c); elseif disc > -0.001; cs = (-T2 )/(2*T1); cs = cs| ((-T2)/(2*T1)); c = s*cs; clearg adg,v1,v2; v1 = (1/Rf)|1|(C[1]/S); adg = inv(A)*v1; _a = adg[1]; _d = adg[2]; _g = adg[3]; v2 = (1/Rf)|1|(C[2]/S); adg = inv(A)*v2; _a = _a|adg[1]; _d = _d|adg[2]; _g = _g|adg[3]; retp(c); else; "error: can't solve for call price 1"; retp(0|0); endif; else; "Findcb: error. one of c or srdes must < 0 to indicate return"; endif; endp; /* ************* FINDCB ***************/ /*-----------------------------------------------------------------*/ /* MULTIPLE OPTION CASE */ /* proc to return either sharpe ratio or a call price (MC[1]) when you are allowed to trade on the bond, the stock and MOPT options. (we know the prices of the last MOPT-1 options) (the last MOPT-1 options have to be ordered by increasing strikes) no positivity . feed -1 for the one you want returned as a function of the other. One of MC[1] or srdes MUST be = -1. Returns: 1) if srdes < 0 and MC[1] > 0 returns sharpe ratio of stock, bond and options all together | sharpe ratio of option and bond alone. 2) if srdes > 0 and MC[1] < 0 returns upper | lower bounds for the first call option Inputs: MX: vector of strikes of call options MC: vector of prices of options */ /*------------------------------------------------------------------*/ proc(1) = findmcb(S,T,sig,Rf,mu,MX,MOPT,MC,srdes); local MAXX,c; clearg f,Nf,Nfp,Nfm,A,B1,B2,B3,B4,B5,B6,B,v,SR,W,T1,T2,T3,disc,cs,erc, sigrc,shrc,a0,a1,a2,i,j,k; /* matrix A is E[x'x] */ A = zeros(MOPT+2,MOPT+2); v = zeros(MOPT+2,1); /* options */ j=1; do while j<=MOPT; k=j+1; do while k<=MOPT; f = (ln(MAXC(MX[k]|MX[j])/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; A[j,k] = a2 - ((MX[j]+MX[k])/S)*a1 + ((MX[j]*MX[k])/(S^2))*a0; k=k+1; endo; f = (ln(MX[j]/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; A[j,j] = a2 - 2*(MX[j]/S)*a1 + (MX[j]/S)^2*a0; A[j,MOPT+1] = a2 - (MX[j]/S)*a1; A[j,MOPT+2] = a1 - (MX[j]/S)*a0; j=j+1; endo; /* stock */ A[MOPT+1,MOPT+1] = exp((2*mu+sig^2)*T); A[MOPT+1,MOPT+2] = exp(mu*T); /* bond */ A[MOPT+2,MOPT+2] = 1; j=1; do while j<=(MOPT+2); k=1; do while k<=j; A[j,k]=A[k,j]; k=k+1; endo; j=j+1; endo; j=1; do while j<=MOPT; v[j]=MC[j]/S; j=j+1; endo; v[MOPT+1] = 1; v[MOPT+2] = (1/Rf); /* now ready to find sharpe ratio given call prices */ if (MC[1] >= 0) and (srdes < 0); /* if want sr as function of c */ SR = v'inv(A)*v; SR = Rf^2*SR - 1; SR = SR^(1/2); f = (ln(MX[1]/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; Erc = S/MC[1]*(a1 - MX[1]/S*a0); Sigrc = (S/MC[1])^2 * (a2 - a1^2 + (MX[1]/S*a0 - 2*a1)*MX[1]/S*(1-Nfp)); if sigrc >= 0; sigrc = sigrc^0.5; else; "formula for variance of call option return not positive"; sigrc = (-sigrc)^0.5; endif; shrc = abs(Erc-Rf)/sigrc; retp(sr|shrc); /* call price for given sharpe ratio */ elseif (MC[1] < 0) and (srdes >= 0); BigA = ((1+(srdes)^2)^0.5)/Rf; W = inv(A); T1 = W[1,1]; T2 = 2*W[1,2:MOPT+2]*v[2:MOPT+2]; T3 = (v[2:MOPT+2])'*W[2:MOPT+2,2:MOPT+2]*v[2:MOPT+2]; T3 = T3-BigA^2; disc = T2^2 - 4*T1*T3; if disc >= 0 ; cs = (-T2 + disc^0.5)/(2*T1); cs = cs| ((-T2 - disc^0.5)/(2*T1)); c = s*cs; _alpha=( ( inv(A)*(cs[1]|v[2:Mopt+2]) ) ~ ( inv(A)*(cs[2]|v[2:Mopt+2]) )); retp(c); elseif disc > -0.001; cs = (-T2 )/(2*T1); cs = cs| ((-T2)/(2*T1)); c = s*cs; _alpha=( ( inv(A)*(cs[1]|v[2:Mopt+2]) ) ~ ( inv(A)*(cs[2]|v[2:Mopt+2]) )); retp(c); else; "error: can't solve for call price 2"; retp(0|0); endif; else; "Findcb: error. one of c or srdes must < 0 to indicate return"; endif; endp; /* ------------------------------------ */ /* versions to find strictly positive m */ /* ------------------------------------ */ /* IS THIS SUPERFLUOUS /*------------------------------------------------------------------*/ /* proc findcbp mimimizes r+ to find sharpe ratio does not yet find call prices given sharpe ratio set global _ad starting value before you call this proc the first time. then it will use the last minimum as a starting value. */ /*------------------------------------------------------------------*/ proc(1) = findcbp(X,S,T,sig,Rf,mu,C,srdes); local Erp2,g,retcode,sr; if (srdes < 0) and (c > 0); {_ad,Erp2,g,retcode} = optmum(&rplus,_ad); /* find min E(R+^2) */ if retcode /= 0; "findcbp: optmum terminated abnormally retcode = " retcode; endif; sr = (Rf^2 / Erp2 - 1 )^0.5; retp(sr); elseif (srdes > 0) and (c <= 0); "findbcp: error. not implemented to find c given sr yet"; retp(0|0); endif; endp; /* -------------------------------------------- */ /* subproc rplus for opmtmum to minimize Warning: Optmum does not allow extra parameters in proc so uses globals for S,C,S,Rf everything but a,d ! */ /* -------------------------------------------- */ proc(1) = rplus(ad); local a,d,Rs,ans; /* "rplus: trying a d" ad';*/ a = ad[1]; d = ad[2]; if (a*S/C + d) < 0; /* upper limits of integration */ Rs = (a*X - (1-a-d)*C*Rf)/(a*S+d*C); else; Rs = 10000; endif; ans = d^2*alph(2,0,Rs) + 2*d*(1-a-d)*Rf*alph(1,0,Rs) + (1-a-d)^2*Rf^2*alph(0,0,Rs) + a^2*(S/C)^2* ( alph(2,X/S,Rs) - 2*X/S*alph(1,X/S,Rs) + (X/S)^2*alph(0,X/S,Rs) ) + 2*a*d*S/C*(alph(2,X/S,Rs) - X/S*alph(1,X/S,Rs)) + 2*a*(1-a-d)*S/C*Rf*(alph(1,X/S,Rs) - X/S*alph(0,X/S,Rs)); /*" answer" ans;*/ retp(ans); endp; /* --------------------------------------- */ /* subproc alph(i,a,b) returns alpha_i(a,b) in paper notation. Use b>= 1000000 for an upper limit of infinity WARNING: Uses GLOBAL values of mu, sig T. Make sure these are set right! */ /* ---------------------------------------- */ proc(1) = alph(i,a,b); local top,bot,alph; if b <= a; /* return zero for upper bound less than lower */ retp(0); endif; if b < 1000000; /* code for upper bound at infty */ top = (ln(b) - (mu - 1/2*sig^2)*T)/(sig*T^0.5) - i*sig*T^0.5; top = cdfn(top); else; top = 1; endif; if a > 0; /* don't try to take ln(0)! */ bot = (ln(a) - (mu - 1/2*sig^2)*T)/(sig*T^0.5) - i*sig*T^0.5; bot = cdfn(bot); else; bot = 0; endif; alph = top-bot; if i == 1; alph = exp(mu*T)*alph; elseif i == 2; alph = exp((2*mu+sig^2)*T)*alph; endif; retp(alph); endp; END IS SUPERFLUOUS */ /* ------------------------------------------------------------------------- */ /* 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 = 999999 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 >= 999999; 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 >= 999999; 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 to evaluate int[a,b] (alf + bet*R)*Rf f(R) dr for R lognormal with mean mu*T and variance sig^2*T use a = 999999 for + inf ------------------------------------------------------------------------- */ proc(1) = integ0(alf,bet,a,b); clearg 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 >= 999999; 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 >= 999999; 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*Rf * (topphi1 - botphi1) + Rf*bet * exp((mu-1/2*sig^2)*T + sig^2*T/2) * (topphi2 - botphi2) ; retp(ans); endp; /* ---------------------------------------------------------------------- */ /* proc to evaluate int[a,b] (alf + bet*R)*R f(R) dr for R lognormal with mean mu*T and variance sig^2*T use a = 999999 for + inf ------------------------------------------------------------------------- */ proc(1) = integ1(alf,bet,a,b); clearg 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 >= 999999; 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 >= 999999; 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* exp((mu-1/2*sig^2)*T + sig^2*T/2) * (topphi2 - botphi2) + bet * exp(2*((mu-1/2*sig^2)*T + sig^2*T)) * (topphi3 - botphi3) ; retp(ans); endp; /* from lamdel to transformed parameters */ proc(1) = trans(lamdel); local tlam0,tlam1,tdel; tlam0 = lamdel[1]; tlam1 = lamdel[2]; tdel = abs(lamdel[3])^0.5; if Nopt < 2; retp(tlam0|tlam1|tdel); else; retp(tlam0|tlam1|tdel|lamdel[4:Nopt+2]); endif; 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; if Nopt < 2; retp(lam0|lam1|del); else; retp(lam0|lam1|del|tlamdel[4:Nopt+2]); endif; endp; /**********************************************************************/ /* OPTMOPT */ /**********************************************************************/ /*-----------------------------------------------------------------*/ /* MULTIPLE OPTION CASE */ /* proc to compute option price bounds with Positive M you need to max (min) an objective function when you want to compute a call lower (upper) bound NewC[Inx] and you are allowed to trade on the bond, the stock and MOPT options. (we know the prices of NOPT-1 options NewC) The NOPT options have to be ordered by increasing strikes i.e. NewX increasing Remarks: Indx: index of the element in the vectors MX: vector of strikes of call options MC: vector of prices of options which corresponds to the option whose price bounds we want to determine Furthermore, it has to be the case that : NewC[Inx]=0 ----------------------------------------------------------------------*/ proc(1) = optmopt(mlamdel); /* Globals: NewX NewC Inx*/ clearg lamdel, lam0l, lam1l, dell, int1, int2, value, bot, top, alphl, betl, Newlam,ii, XSold; if cols(mlamdel) > 1; "rewrite optval to take cols(mlamdell ) > 1 "; endif; lamdel = back(mlamdel); if Nopt < 2; lam0l = lamdel[1]; lam1l = lamdel[2]; dell = lamdel[3]; Newlam = 1; else; lam0l = lamdel[1]; lam1l = lamdel[2]; dell = lamdel[3]; Newlam = lamdel[4:Nopt+2]; if Inx == 1; Newlam = (1|Newlam); elseif Inx < Nopt; Newlam = (Newlam[1:inx-1]|1|Newlam[inx:Nopt-1]); elseif Inx == Nopt; Newlam = (Newlam|1); endif; endif; 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 */ ii = Nopt; XSold = 99999999; int1 = 0; do while ii>= 1; bot = NewX[ii]/S; top = XSold; alphl = -Newlam[1:ii]'*NewX[1:ii]+lam0l*Rf; betl = (lam1l+S*(ones(1,ii)*Newlam[1:ii])); int1 = int1 + 1/(dell^2)*integ( alphl, betl, bot, top ); XSold = NewX[ii]/S; ii=ii-1; endo; bot = 0; top = NewX[1]/S; alphl = lam0l*Rf; betl = lam1l; int2 = 1/(dell^2)*integ( alphl, betl, bot, top ); else; /* positivity case */ ii = Nopt; XSold = 999999; int1 = 0; do while ii>= 1; if (S*(ones(1,ii)*Newlam[1:ii])+lam1l)/dell > 0; bot = NewX[ii]/S; top = minc(((Newlam[1:ii]'*NewX[1:ii]-lam0l*Rf)/ (lam1l+S*(ones(1,ii)*Newlam[1:ii]))) |XSold); alphl = -Newlam[1:ii]'*NewX[1:ii]+lam0l*Rf; betl = (lam1l+S*(ones(1,ii)*Newlam[1:ii])); int1 = int1 + 1/(dell^2)*integ( alphl, betl, bot, top ); elseif (S*(ones(1,ii)*Newlam[1:ii])+lam1l)/dell < 0; bot = maxc((NewX[ii]/S)| ((Newlam[1:ii]'*NewX[1:ii]-lam0l*Rf)/ (lam1l+S*ones(1,ii)*Newlam[1:ii]))); top = XSold; alphl = -Newlam[1:ii]'*NewX[1:ii]+lam0l*Rf; betl = (lam1l+S*(ones(1,ii)*Newlam[1:ii])); int1 = int1 + 1/(dell^2)*integ( alphl, betl, bot, top ); elseif (S*(ones(1,ii)*Newlam[1:ii])+lam1l)/dell == 0; "figure out S+lam1l)/dell = 0 case"; endif; XSold = NewX[ii]/S; ii=ii-1; endo; if lam1l/dell > 0; bot = 0; top = minc((NewX[1]/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 = NewX[1]/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; value = -dell/2 * (int1 + int2) - lam0l - lam1l - (NewC'*Newlam)- dell/2*BigA^2; if _minit; retp(-value); /* maximize wrt to lambda, if minimize optionvalue */ else; retp(value); endif; endp; /**********************************************************************/ /* Procedure Graphps: rotates the graphs */ /**********************************************************************/ proc(0) = graphps(px,py); _protate = 1; begwind; makewind(9,4,0,2.,1); _ptitlht=0.25; _paxht=0.25; _pnumht=0.2; xy(px,py); endwind; endp;