#include <stdio.h>
#include <math.h>
#include <string.h>
#define ITMAX 100
#define EPS 3.0e-7
#define M_PI 3.1415926535897932
#define max(a,b) ((a)>(b) ? (a):(b))
#define sgn(x) (fabs(x)/(x))
#define Heaviside(a) ((a)>(0.0) ? (a):(0.0))
#define Dirac(a) ((a)==(0.0) ? (1.0):(0.0))
#define NormDist(x) exp(-(x)*(x)/2)/pow((2*M_PI),.5)
#define NormCDF(x) (1.0 +erf((x)) )/2.0
#define stirlng(n) ((n)+.5)*log(n)-(n)+.5*log(2*M_PI)
#define binom(n,k) exp(stirlng(n)-stirlng(k)-stirlng(((n)-(k))))
#define trinom(n,r,s) exp(stirlng(n)-stirling(r)-stirling(s)-stirling(n-r-s))
#define Nfrx(x,p) pow((1+.33267*(x)),-(p))
#define Nfctr(x) (.43618*Nfrx((x),1)-.12017*Nfrx((x),2)+.9373*Nfrx((x),3))
#define NCDF(x) 1-NormDist(x)*Nfctr(x)
#define sum(n,o,f,s) for(n=o;n<(f+1);n++),+=s
main ()
/* Options trial code. Call 917-287-8087 24hrs for help.
Vasos-Peter John Panagiotopoulos II, Samani Intl Entpses
Program named after Thales who Aristotle says invented options */
{float XrcPr,StkPr,dt,r,sig,NShar,NOpts,OptPr,StCom,OpCom,junk;
int n,m,NPerds,NDvdd; double normal(double); float RollDv,tRoll,StkRll;
float BlkSch(float,float,float,float,float,int),
erf(float),gammp(float,float),gammln(float),
factln(int),bico(int,int),trico(int,int,int);
void gser(float *,float,float,float *),
gcf(float *,float,float,float *); float Divdd=0;
printf("Stock price (try 50) $>");scanf("%f",&StkPr);
printf("Excercise stock price (try 45) $>");scanf("%f",&XrcPr);
printf("Expiration time difference (try.25) >");scanf("%f",&dt);
printf("Interest Rate percent (try 6) >");scanf("%f",&r);r=(r/100);
printf("Volatility (get bs7.62) (try.4472)>");scanf("%f",&sig);
/* Future versions should calc fresh stdev from db sql call
Benninga Num Tec Financ MIT 89 p149 suggests finding underlying
volatility by Newton Raphson iteration of (f(sig)-C) from data */
printf("Discretisation intervals (try 9) >");scanf("%i",&NPerds);
printf("How many dividends during (try 1) >");scanf("%i",&NDvdd);
for (n=1;n<(NDvdd+1);n++) {float thsDiv,thsTim;
printf("\t Enter dividend %i> (try 5)",n);scanf("%f",&thsDiv);
printf("\t\t Enter time to payout> (try .2)");scanf("%f",&thsTim);
if (NDvdd==1) {RollDv=thsDiv;tRoll=thsTim;StkRll=StkPr;};
Divdd+=thsDiv*exp(-r*thsTim);}
{float NetInv,Cap_Gn,TotInc,ROI,AnROI,NetPr,CshRec,OROI,OAnROI;
printf("Shares bought (try 10) >");scanf("%f",&NShar);
printf("Options written (try 7) >");scanf("%f",&NOpts);
printf("Option price (try 5) $>");scanf("%f",&OptPr);
printf("Stock commission (try 1) $>");scanf("%f",&StCom);
printf("Option commission (try .1) $>");scanf("%f",&OpCom);
NetInv=(NShar*StkPr)-(100*NOpts*OptPr)+(StCom+OpCom); /*Divdd w o exp*/
Cap_Gn=NShar*XrcPr-StCom-NetInv;TotInc=Cap_Gn+Divdd;ROI=100*TotInc/NetInv;
AnROI=ROI*dt;NetPr=(100*NOpts*OptPr)-OpCom;CshRec=NetPr+Divdd;
OROI=100*CshRec/NetInv;OAnROI=OROI*dt;
printf("IF CALLED Return %f pct APR %f \n\t NOCALL Return %f pct APR %f",
ROI,AnROI,OROI,OAnROI);
};printf("\n(Theoretical valuation below assumes no quantity or commission)");
/* Basic nondividend Black Scholes and Cox Ross Rubinstein analysis is
from Copeland Weston Financial Thoery AW 1983 pp 252-259. Initial debug
used their test data. Thanks to Pr Nitzan Weiss for 3/83 ref to book.
Also thanks to Martin Liebowitz of Sali for 10/90 ref to Merton.
Dividend and put refinements due to Natenberg Probus 336-343.
Extremal risk anal from NYIF Acct Exec Training manual. */
StkPr-=Divdd; /* Gibson MGH 91 p 194 and Natenberg p336 */
junk=BlkSch(StkPr,XrcPr,dt,r,sig,0); /*wont print*/
/* fprintf("nul:",junk) */
{float UpMult,DnMult,HgProb,rBn;float CBnVal=0.0,PBnVal=0.0;
rBn=1+(pow((1+r),(dt/NPerds))-1);
UpMult=exp(sig*pow(dt/NPerds,.5)); DnMult=1/UpMult;
HgProb=(rBn-DnMult)/(UpMult-DnMult);
for (n=0;n<(NPerds+1);n++) {float bkern1,bkern2;
bkern1=bico((int) NPerds,(int) n)
*pow(HgProb,(float) n)*pow((1-HgProb),(float) (NPerds-n));
bkern2=pow(UpMult,(float) n)*pow(DnMult,(float) (NPerds-n))*StkPr-XrcPr;
CBnVal+=bkern1*max(0.0,bkern2); PBnVal+=bkern1*max(0.0,-bkern2);
};
CBnVal=CBnVal/rBn;PBnVal=PBnVal/rBn;
printf("\n Cox Ross Rubinstein Call value %5.2f Put %5.2f", CBnVal,PBnVal);
};
{/*BEGIN WHALEY MODEL*/
float S1,S2,eigQd2,eigQd1,h,del,RHS,CWhVal,PWhVal;
float b; b=r; /*b is 0 for futures and r for stocks*/ S1=XrcPr;S2=0;
while ((fabs(S1-S2)/S1)>0.04)
{eigQd2=.5*(1-(2*b/(sig*sig)))/2
+.5*pow((pow((2*b/(sig*sig)-1),2)+8*r/(sig*sig*(1-exp(-r*dt)))),.5);
if (S2==0) S1=XrcPr/(1-(1/eigQd2)); else S1=S2;
/* seed not Natenberg 312-341-4432 but Whaley 919-660-7781 JF42#2p310*/
h=log(S1/(XrcPr*exp((b-r)*dt)))/(sig*pow(dt,.5))+sig*pow(dt,.5)/2;
del=exp((b-r)*dt)*((float) normal(h))*(1-(1/eigQd2))
+ (1-(exp((b-r)*dt)*NormDist(h)))/sig/eigQd2/pow(dt,0.5);
/* Black Scholes basis applies only for b=r */
RHS=BlkSch(S1,XrcPr,dt,r,sig,1)
+(1-(exp((b-r)*dt)*((float) normal(h))))*S1/eigQd2;
S2=(XrcPr+RHS-del*S1)/(1-del);
}; /*Natenberg forgot pow eigQd2 below */
if (StkPr<S2) CWhVal=BlkSch(StkPr,XrcPr,dt,r,sig,1)
+pow((StkPr/S2),eigQd2)*(S2/eigQd2)
*(1-(exp((b-r)*dt))*((float) normal(h)));
else CWhVal=StkPr-XrcPr;
S1=XrcPr;S2=0; /* start again for put */
while ((fabs(S1-S2)/S1)>0.04)
{eigQd1=.5*(1-(2*b/(sig*sig)))/2
-.5*pow((pow((2*b/(sig*sig)-1),2)+8*r/(sig*sig*(1-exp(-r*dt)))),.5);
if (S2==0) S1=XrcPr/(1-(1/eigQd1)); else S1=S2;
/* seed not Natenberg 312-341-4432 but Whaley 919-660-7781 JF42#2p310*/
h=log(S1/(XrcPr*exp((b-r)*dt)))/(sig*pow(dt,.5))+sig*pow(dt,.5)/2;
del=-exp((b-r)*dt)*((float) normal(-h))*(1+(1/eigQd1))
- (1-(exp((b-r)*dt)*NormDist(-h)))/sig/eigQd1/pow(dt,0.5);
/* Black Scholes basis applies only for b=r */
RHS=BlkSch(S1,XrcPr,dt,r,sig,-1)
-(1-(exp((b-r)*dt)*((float) normal(-h))))*S1/eigQd1;
S2=(XrcPr-RHS+del*S1)/(1+del);
/* I changed s-x to x-s & Whaley corrxd p372 Smith HbkFE H&R'90 */
};
if (StkPr>S2) PWhVal=BlkSch(StkPr,XrcPr,dt,r,sig,-1)
-pow((StkPr/S2),eigQd1)*(S2/eigQd1)
*(1-(exp((b-r)*dt))*((float) normal(-h)));
else PWhVal=StkPr-XrcPr;
printf("\n Whaley Barone Adesi Call %5.2f Put %5.2f",CWhVal,PWhVal);
}; /*END WHALEY MODEL*/
if ((NDvdd==1) && (0==0)) /*Jarrow Rudd Op Prcg DJI83 p146*/
{ /*BEGIN ROLL MODEL*/
float RollVl,StkMid1,StkMid2,StkMid,RtFct1,RtFct2,a1,a2,b1,b2;
float BvNorm(float,float,float),BlkSch(float,float,float,float,float,int);
double normal(double);
StkMid2=StkRll*1.1;StkMid1=StkRll*.95;StkMid=StkMid2;
do{ RtFct2=StkMid2-XrcPr+RollDv-BlkSch(StkMid2,XrcPr,tRoll,r,sig,+1);
RtFct1=StkMid1-XrcPr+RollDv-BlkSch(StkMid1,XrcPr,tRoll,r,sig,+1);
StkMid+=(StkMid1-StkMid2)*RtFct2/(RtFct2-RtFct1);
StkMid1=StkMid2; StkMid2=StkMid;
} while (fabs((StkMid-StkMid2)/StkMid)>0.04);
printf("\n Roll Underlying Stock Price %5.2f",StkMid);
a1=(log(StkPr/XrcPr)+(r+sig*sig/2)*dt)/sig/pow(dt,.5);
a2=a1-sig*pow(dt,.5);
b1=(log(StkPr/StkMid)+(r+sig*sig/2)*tRoll)/sig/pow(tRoll,.5);
b2=b1-sig*pow(dt,.5);
printf("\n DBG: ROLL PARAMS %5.2f %5.2f %5.2f %5.2f",a1,a2,b1,b2);
RollVl=(StkPr)
*((float) normal(a1)-(float) normal(b1)
+BvNorm(b1,a1,pow((tRoll/dt),.5)))
+XrcPr*exp(-r*dt)*
((float) normal(a2)-BvNorm(a2,b2,pow((tRoll/dt),.5)))
-(XrcPr-RollDv)*exp(r*tRoll)* ((float) normal(b2));
printf("\n Roll Model Call Value %5.2f (coding in progress)",RollVl);
};/*END ROLL MODEL*/
{ /*BEGIN TRUNOM*/
float rTr,prUp,prNo,prDn,DistMean,DistVar,DnMult,UpMult,StepTweak;
float CTrVal=0.0; /* Trunomial Boyle JFQA 3'88 via Pr Thanos Episcopos */
rTr=1+(pow((1+r),(dt/NPerds))-1); DistMean=exp(r*dt/NPerds);
DistVar=pow(DistMean,2)*(exp(pow(sig,2)*dt/NPerds)-1); StepTweak=1.5;
UpMult=exp(StepTweak*sig*pow(dt/NPerds,0.5)); DnMult=1/UpMult;
prUp=((DistVar+pow(DistMean,2)-DistMean)*UpMult-(DistMean-1))
/(UpMult-1)/(pow(UpMult,2)-1);
prDn=(pow(UpMult,2)*(DistVar+pow(DistMean,2)-DistMean)
-pow(UpMult,3)*(DistMean-1))
/(UpMult-1)/(pow(UpMult,2)-1);
prNo=1-prDn-prUp; n=-1;
do{n+=1; for (m=0;m<=NPerds;m++) {float tkern1,tkern2;
tkern1=trico(NPerds,n,m)
*pow(prUp, m)*pow(prDn, n)*pow(prNo,(NPerds-n-m));
tkern2= pow(UpMult,m)*pow(DnMult,n)*pow(1, NPerds-n-m)*StkPr-XrcPr;
CTrVal+=tkern1*max(0.0,tkern2);
};} while (n<=NPerds); CTrVal=CTrVal/(1+r);
printf("\n Trunomial Call value %5.2f ", CTrVal);
}; /* END TRUNOM */
{ /* BEGIN MTCRLO */
float ThisCVal,NxtStPr,LstStPr,rATdt,rPerd,RandVar;/*BoyleMtCrloJFE77*/
float AvgCVal=0;LstStPr=StkPr;m=0;RandVar=sig;rATdt=pow((1+r),dt)-1;
do{ for (n=0;n<=NPerds;n++)
{RandVar=1-fmod((double) (0.5*RandVar+1)*325476, (double) 100)/50;
NxtStPr=LstStPr*exp(rATdt-sig*sig/2+sig*RandVar); /* lognormal */
}; m+=1; ThisCVal=max(0.0,NxtStPr*exp(-rATdt*dt)-XrcPr);
AvgCVal=((m-1)*AvgCVal+ThisCVal)/m;
} while (fabs((ThisCVal-AvgCVal)/(StkPr-XrcPr))>=(r*sig/NPerds)
|| m<NPerds);
printf("\n Monte Carlo Call value %5.2f at %5i iterations", AvgCVal,m);
}; /*END MTCRLO*/
}; /*END MAIN*/
float BlkSch(float StkPr,float XrcPr,float dt,float r,float sig,int act)
{float d1,d2,BSCall,BSPut; double normal(double);
float erf(float),gammp(float,float),gammln(float),factln(int),bico(int,int);
void gser(float *,float,float,float* ),gcf(float *,float,float,float *);
d1=(log(StkPr/XrcPr)+((pow(sig,2)/2)+r)*(dt))/(sig*pow(dt,.5));
d2=d1-sig*pow(dt,.5);
BSCall=StkPr*((float) normal( d1))-XrcPr*((float) normal( d2))*exp(-r*dt);
BSPut=-StkPr*((float) normal(-d1))+XrcPr*((float) normal(-d2))*exp(-r*dt);
if (act==-1) return (BSPut); /* Itsa put */ /*act was a string */
if (act== 1) return (BSCall); /* Itsa call */ /* used strcmp.[1] */
if (act== 0) /* Do it all */
{float gamma,vega,CDelta,PDelta,CRho,PRho,CTheta,PTheta;
printf("\n Black Scholes Call Value is %5.2f Put %5.2f",BSCall,BSPut);
/* sensitivities due to Natenberg also Jarrow Rudd DJI 83 p 117 */
CDelta=(float) normal(d1);PDelta=-(float) normal(-d1);
printf("\n - SENSITIVITIES:");
printf("\n\t Hedge Shares Per Option (delta) Call %f Put %f",
CDelta,PDelta);
gamma=NormDist(d1)/(StkPr*sig*pow(dt,.5));
vega=StkPr*NormDist(d1)*pow(dt,.5);
printf("\n\t Curvature (gamma) %5.2f VolSens (vega) %f",gamma,vega);
/*delta times StkPr/BSCall/StockBeta gives CallBeta JR p121*/
CRho= dt*XrcPr*((float) normal( d2))*exp(-r*dt);
PRho=-dt*XrcPr*((float) normal(-d2))*exp(-r*dt);
printf("\n\t Interest sensitivity (rho) Call %f Put %f",CRho,PRho);
CTheta=StkPr*sig*NormDist(d1)/(2*pow(dt,.5))+r*(CRho/dt);
PTheta=StkPr*sig*NormDist(d1)/(2*pow(dt,.5))+r*(PRho/dt);
printf("\n\t Time Decay (theta) Call %5.2f Put %5.2f",CTheta,PTheta);
printf("\n\t InMoneyRatio %f Maturity %f",((StkPr-XrcPr)/StkPr),dt);
printf("\n\t (High +/- InMoney or low maturity underestimates)");
printf("\n\t (Low positive InMoney or hi maturity overestimates)");
return (0.0);
}; /*Output but return nothing*/
} /* END Black Scholes */
float BvNV(float Q,float H)
/* Integral V used to find biNormal. p120 Greenwood Hartley Guide to Tables
Princeton 1962 refer to D B Owen NBS 1959. IMSL (thanx to NY Fed's Dr Kell
for ref) and Jarrow Rudd use other Owen work. Integrated via MuMath
symbolic manipulation DEFINT and TAYLOR algorithms using methodology of
p266 Bronshtein Semendyayev Math Hbk VNR Nauka 1985 */
{return( Q*H/(4*M_PI) -Q*pow(H,3)/(16*M_PI) +Q*pow(H,5)/(96*M_PI)
-pow(Q,3)*H/(48*M_PI) +pow(Q,3)*pow(H,3)/(144*M_PI)
-pow(Q,3)*pow(H,5)/(768*M_PI) +pow(Q,5)*H/(480*M_PI)
-pow(Q,5)*pow(H,3)/(1280*M_PI) +pow(Q,5)*pow(H,5)/(6400*M_PI) ); }
float BvNorm(float x,float y, float rho)
{double normal(double);
printf("\n DBG: BVNORM CALLED %5.2f %5.2f %5.2f", x, y, rho);
if (rho<0)
{if (x<0) return((float) normal(y) -BvNorm(-x, y,-rho));
if (y<0) return((float) normal(x) -BvNorm( x,-y,-rho));
};
if ((x<0)&&(y<0))
return ((float) normal(x) +(float) normal(y)-1+BvNorm(-x,-y, rho));
if ((x*y*rho)>0)
{float rhox,rhoy,drho;drho=(1-sgn(x)*sgn(y))/4;
rhox=(rho*x-y)*sgn(x)/sqrt(x*x-2*rho*x*y+y*y);
rhoy=(rho*y-x)*sgn(y)/sqrt(x*x-2*rho*x*y+y*y);
return(BvNorm(x,0,rhox)+BvNorm(y,0,rhoy)-drho);
};
{float tau,nu; tau=pow((1-rho*rho),-0.5);nu=rho*tau;
printf("\n DBG: BVNORM KERNEL %5.2f %5.2f",tau,rho);
return( BvNV(x,(tau*y-nu*x))+BvNV(y,(tau*x-nu*y))-0.25
+0.5*((float) normal(x) +(float) normal(y))
+asin(rho)/2/M_PI );
};
}
/*Routines from Press Flannery NumRec Camb 86 and Purdum Toolkit Que 89*/
float erf(float x)
{float gammp(float,float);return x<0.0 ? -gammp(0.5,x*x):gammp(0.5,x*x);}
float gammp(float a,float x)
{float gamser,gammcf,gln;
void gser(float *,float,float,float *),gcf(float *,float,float,float *);
if (x < (a+1.0)) {gser(&gamser,a,x,&gln); return gamser;}
else {gcf(&gammcf,a,x,&gln);return 1.0-gammcf;} }
float gammln(float xx)
{ float x,tmp,ser;int j;
static float cof[6]={76.18009173,-86.50532033,24.01409822,
-1.231739516,0.120858003e-2,-0.536382e-5};
x=xx-1.0;tmp=x+5.5;tmp-=(x+0.5)*log(tmp);ser=1.0;
for (j=0;j<=5;j++) {x+=1.0;ser+=cof[j]/x;}
return -tmp+log(2.50662827465*ser); }
float factln(int n) {static float a[101];float gammln(float);
if (n<=1) return 0.0;
if (n<=100) return a[n]?a[n]:(a[n]=gammln(n+1.0));else return gammln(n+1.0);}
float bico(int n,int k)
{float factln(int);return floor(0.5+exp(factln(n)-factln(k)-factln(n-k)));}
float trico(int n,int k,int m) /* Trinomial Coefficient */
{float factln(int);
return floor(0.5+exp(factln(n)-factln(k)-factln(m)-factln(n-k-m)));}
void gser(float *gamser,float a,float x,float *gln) /* astrsk is by ref */
{int n;float sum,del,ap;float gammln(float); *gln=gammln(a);/*pluck ref w ampsd*/
if (x<=0.0) {if (x<0.0) printf("\a\n NEG x GSER\n");*gamser=0.0;return;}
else {ap=a; del=sum=1.0/a;
for (n=1;n<=ITMAX;n++) { ap += 1.0;del *= x/ap;sum += del;
if (fabs(del)<fabs(sum)*EPS) {*gamser=sum*exp(-x+a*log(x)-(*gln));return;}
} printf("\a\n a too large, ITMAX too small in GSER\n"); return;} }
void gcf(float *gammcf,float a,float x,float *gln)
{int n; float gold=0.0,g,fac=1.0,b1=1.0;
float b0=0.0,anf,ana,an,a1,a0=1.0; float gammln(float); *gln=gammln(a); a1=x;
for (n=1;n<=ITMAX;n++) {
an=(float) n;ana=an-a;a0=(a1+a0*ana)*fac;b0=(b1+b0*ana)*fac;
anf=an*fac;a1=x*a0+anf*a1;b1=x*b0+anf*b1;
if (a1) {fac=1.0/a1;g=b1*fac;
if (fabs((g-gold)/g)<EPS) {*gammcf=exp(-x+a*log(x)-(*gln))*g;return;}
gold=g;}
} printf("\a\n a too large, ITMAX too small in routine GCF\n"); }
double normal(double x) {double w, y, z;
if (x == 0.0) z = 0.0;
else {y=fabs(x)/2.0;
if (y >= 3.0) {z = 1.0;}
else {
if (y < 1.0) {w = y * y;
z = ((((((((0.000124818987 * w
-.001075204047) * w + .005198775019) * w
-.019198292004) * w + .059054035642) * w
-.151968751364) * w + .319152932694) * w
-.531923007300) * w + .797884560593) * y * 2.0;
} else {y = y - 2.0;
z = (((((((((((((-.000045255659 * y
+.000152529290) * y - .000019538132) * y
-.000676904986) * y + .001390604284) * y
-.000794620820) * y - .002034254874) * y
+.006549791214) * y - .010557625006) * y
+.011630447319) * y - .009279453341) * y
+.005353579108) * y - .002141268741) * y
+.000535310849) * y + .999936657524;
}}} if (x>0.0) return (z+1.0)/2.0;else return (1-z)/2.0;}
/*Use date routines to get real dates and find fractions in code*/
/*
int julian(int day, int month, int year)
{static int runsum[] = {0,31,59,90,120,151,181,212,243,273,304,334,365};
int total; total = runsum[month - 1] + day;
if (month > 2) {total += leapyear(year);}return (total); }
int leapyear(int year)
{if (year%4==0 && year%100 != 0 || year%400 == 0) return 1;else return 0;}
char *datename(int day,int month,int year)
{ static char *days[] = {"Sun","Mon","Tue","Wed","Thur","Fri","Sat","Sun"};
int index; if (year<100) year+=1900;
if (month>2) {month-=2;} else {month += 10;year--;}
index=((13*month-1)/5)+day+(year%100)+((year%100)/4)
+((year/100)/4)-2*(year/100)+77; index=index-7*(index/7);
return days[index]; }
unsigned int datediff(int d1, int m1, int y1, int d2, int m2, int y2)
{int i, max, min, t1, t2, years = 0;
t1 = julian(d1, m1, y1); t2 = julian(d2, m2, y2);
if (y1!=y2) {if (y1>y2) {max=y1;min=y2;t2=julian(31,12,y2)-t2;}
else {max=y2;min=y1;t1=julian(31,12,y2)-t1;}
for (i = max; i > min + 1; i--) {years += 365 + leapyear(i);}
} else {t1=julian(d1,m1,y1);t2=julian(d2,m2,y2);return(fabs(t1-t2));}
return (t1 + t2 + years + 1); }
*/