Actual source code: ex2.c
petsc-3.4.5 2014-06-29
1: static char help[] = "Tests PetscRandom functions.\n\n";
3: #include <petscsys.h>
5: #define PETSC_MAXBSIZE 40
6: #define DATAFILENAME "ex2_stock.txt"
8: struct himaInfoTag {
9: int n;
10: double r;
11: double dt;
12: int totalNumSim;
13: double *St0;
14: double *vol;
15: };
16: typedef struct himaInfoTag himaInfo;
18: /* function protype */
19: PetscErrorCode readData(MPI_Comm comm,himaInfo *hinfo);
20: double mcVal(double St, double r, double vol, double dt, double eps);
21: void exchange(double *a, double *b);
22: double basketPayoff(double vol[], double St0[], int n, double r,double dt, double eps[]);
23: void stdNormalArray(double *eps, int size,PetscRandom ran);
24: unsigned long divWork(int id, unsigned long num, int np);
26: /*
27: Contributed by Xiaoyan Zeng <zengxia@iit.edu> and Liu, Kwong Ip" <kiliu@math.hkbu.edu.hk>
29: Example of usage:
30: mpiexec -n 4 ./ex2 -num_of_stocks 30 -interest_rate 0.4 -time_interval 0.01 -num_of_simulations 10000
31: */
35: int main(int argc, char *argv[])
36: {
37: /* double payoff; */
38: double r,dt;
39: int n;
40: unsigned long i,myNumSim,totalNumSim,numdim;
41: double *vol, *St0, x, totalx;
42: int np,myid;
43: double *eps;
44: himaInfo hinfo;
45: PetscRandom ran;
47: PetscBool flg;
49: PetscInitialize(&argc,&argv,(char*)0,help);
50: #if defined(PETSC_USE_COMPLEX)
51: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This example does not work for scalar type complex\n");
52: #endif
53: PetscRandomCreate(PETSC_COMM_WORLD,&ran);
54: #if defined(PETSC_HAVE_SPRNG)
55: PetscRandomSetType(ran,PETSCSPRNG);
56: #elif defined(PETSC_HAVE_RAND)
57: PetscRandomSetType(ran,PETSCRAND);
58: #endif
59: PetscRandomSetFromOptions(ran);
61: MPI_Comm_size(PETSC_COMM_WORLD, &np); /* number of nodes */
62: MPI_Comm_rank(PETSC_COMM_WORLD, &myid); /* my ranking */
64: PetscOptionsHasName(NULL, "-check_generators", &flg);
65: if (flg) {
66: PetscRandomGetValue(ran,(PetscScalar*)&r);
67: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] rval: %g\n",myid,r);
68: PetscSynchronizedFlush(PETSC_COMM_WORLD);
69: }
71: hinfo.n = 31;
72: hinfo.r = 0.04;
73: hinfo.dt = 1.0/12; /* a month as a period */
74: hinfo.totalNumSim = 1000;
76: PetscOptionsGetInt(NULL,"-num_of_stocks",&(hinfo.n),NULL);
77: if (hinfo.n <1 || hinfo.n > 31) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only 31 stocks listed in stock.txt. num_of_stocks %D must between 1 and 31",hinfo.n);
78: PetscOptionsGetReal(NULL,"-interest_rate",&(hinfo.r),NULL);
79: PetscOptionsGetReal(NULL,"-time_interval",&(hinfo.dt),NULL);
80: PetscOptionsGetInt(NULL,"-num_of_simulations",&(hinfo.totalNumSim),NULL);
82: n = hinfo.n;
83: r = hinfo.r;
84: dt = hinfo.dt;
85: totalNumSim = hinfo.totalNumSim;
86: vol = hinfo.vol = (double*)malloc(sizeof(double)*(2*n+1));
87: St0 = hinfo.St0 = hinfo.vol + n;
88: readData(PETSC_COMM_WORLD,&hinfo);
90: numdim = n*(n+1)/2;
91: if (numdim%2 == 1) numdim++;
92: eps = (double*)malloc(sizeof(double)*numdim);
94: myNumSim = divWork(myid,totalNumSim,np);
96: x = 0;
97: for (i=0; i<myNumSim; i++) {
98: stdNormalArray(eps,numdim,ran);
99: x += basketPayoff(vol,St0,n,r,dt,eps);
100: }
102: MPI_Reduce(&x, &totalx, 1, MPI_DOUBLE, MPI_SUM,0,PETSC_COMM_WORLD);
103: /* payoff = exp(-r*dt*n)*(totalx/totalNumSim);
104: PetscPrintf(PETSC_COMM_WORLD,"Option price = $%.3f using %ds of %s computation with %d %s for %d stocks, %d trading period per year, %.2f%% interest rate\n",
105: payoff,(int)(stop - start),"parallel",np,"processors",n,(int)(1/dt),r); */
107: free(vol);
108: free(eps);
109: PetscRandomDestroy(&ran);
110: PetscFinalize();
111: return 0;
112: }
114: void stdNormalArray(double *eps, int size, PetscRandom ran)
115: {
116: int i;
117: double u1,u2,t;
120: for (i=0; i<size; i+=2) {
121: PetscRandomGetValue(ran,(PetscScalar*)&u1);CHKERRABORT(PETSC_COMM_WORLD,ierr);
122: PetscRandomGetValue(ran,(PetscScalar*)&u2);CHKERRABORT(PETSC_COMM_WORLD,ierr);
124: t = sqrt(-2*log(u1));
125: eps[i] = t * cos(2*PETSC_PI*u2);
126: eps[i+1] = t * sin(2*PETSC_PI*u2);
127: }
128: }
131: double basketPayoff(double vol[], double St0[], int n, double r,double dt, double eps[])
132: {
133: double Stk[PETSC_MAXBSIZE], temp;
134: double payoff;
135: int maxk,i,j;
136: int pointcount=0;
138: for (i=0;i<n;i++) Stk[i] = St0[i];
140: for (i=0;i<n;i++) {
141: maxk = 0;
142: for (j=0;j<(n-i);j++) {
143: Stk[j] = mcVal(Stk[j],r,vol[j],dt,eps[pointcount++]);
144: if ((Stk[j]/St0[j]) > (Stk[maxk]/St0[maxk])) maxk = j;
145: }
146: exchange(Stk+j-1,Stk+maxk);
147: exchange(St0+j-1,St0+maxk);
148: exchange(vol+j-1,vol+maxk);
149: }
151: payoff = 0;
152: for (i=0; i<n; i++) {
153: temp = (Stk[i]/St0[i]) - 1;
154: if (temp > 0) payoff += temp;
155: }
156: return payoff;
157: }
161: PetscErrorCode readData(MPI_Comm comm,himaInfo *hinfo)
162: {
163: int i;
164: FILE *fd;
165: char temp[50];
167: PetscMPIInt rank;
168: double *v = hinfo->vol, *t = hinfo->St0;
169: int num=hinfo->n;
172: MPI_Comm_rank(comm,&rank);
173: if (!rank) {
174: PetscFOpen(PETSC_COMM_SELF,DATAFILENAME,"r",&fd);
175: for (i=0;i<num;i++) fscanf(fd,"%s%lf%lf",temp,v+i,t+i);
176: fclose(fd);
177: }
178: MPI_Bcast(v,2*num,MPI_DOUBLE,0,PETSC_COMM_WORLD);
179: /* PetscPrintf(PETSC_COMM_SELF,"[%d] vol %g, ... %g; St0 %g, ... %g\n",rank,hinfo->vol[0],hinfo->vol[num-1],hinfo->St0 [0],hinfo->St0[num-1]); */
180: return(0);
181: }
183: void exchange(double *a, double *b)
184: {
185: double t;
187: t = *a;
188: *a = *b;
189: *b = t;
190: }
192: double mcVal(double St, double r, double vol, double dt, double eps)
193: {
194: return (St * exp((r-0.5*vol*vol)*dt + vol*sqrt(dt)*eps));
195: }
197: unsigned long divWork(int id, unsigned long num, int np)
198: {
199: unsigned long numit;
201: numit = (unsigned long)(((double)num)/np);
202: numit++;
203: return numit;
204: }