#include #include #include #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 (StkPr0.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) || m0) {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)= 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); } */