Actual source code: ex2.c
petsc-3.14.6 2021-03-30
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: PetscInt n;
10: PetscReal r;
11: PetscReal dt;
12: PetscInt totalNumSim;
13: PetscReal *St0;
14: PetscReal *vol;
15: };
16: typedef struct himaInfoTag himaInfo;
18: PetscErrorCode readData(MPI_Comm,himaInfo *);
19: PetscReal mcVal(PetscReal, PetscReal, PetscReal, PetscReal, PetscReal);
20: void exchangeVal(PetscReal*, PetscReal*);
21: PetscReal basketPayoff(PetscReal[], PetscReal[], PetscInt, PetscReal,PetscReal, PetscReal[]);
22: PetscErrorCode stdNormalArray(PetscReal*, PetscInt,PetscRandom);
23: PetscInt divWork(PetscMPIInt, PetscInt, PetscMPIInt);
25: /*
26: Contributed by Xiaoyan Zeng <zengxia@iit.edu> and Liu, Kwong Ip" <kiliu@math.hkbu.edu.hk>
28: Example of usage:
29: mpiexec -n 4 ./ex2 -num_of_stocks 30 -interest_rate 0.4 -time_interval 0.01 -num_of_simulations 10000
30: */
32: int main(int argc, char *argv[])
33: {
34: PetscReal r,dt;
35: PetscInt n;
36: unsigned long i,myNumSim,totalNumSim,numdim;
37: PetscReal *vol, *St0, x, totalx;
38: PetscMPIInt size,rank;
39: PetscReal *eps;
40: himaInfo hinfo;
41: PetscRandom ran;
44: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
45: PetscRandomCreate(PETSC_COMM_WORLD,&ran);
46: PetscRandomSetFromOptions(ran);
48: MPI_Comm_size(PETSC_COMM_WORLD, &size); /* number of nodes */
49: MPI_Comm_rank(PETSC_COMM_WORLD, &rank); /* my ranking */
51: hinfo.n = 31;
52: hinfo.r = 0.04;
53: hinfo.dt = 1.0/12; /* a month as a period */
54: hinfo.totalNumSim = 1000;
56: PetscOptionsGetInt(NULL,NULL,"-num_of_stocks",&(hinfo.n),NULL);
57: 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);
58: PetscOptionsGetReal(NULL,NULL,"-interest_rate",&(hinfo.r),NULL);
59: PetscOptionsGetReal(NULL,NULL,"-time_interval",&(hinfo.dt),NULL);
60: PetscOptionsGetInt(NULL,NULL,"-num_of_simulations",&(hinfo.totalNumSim),NULL);
62: n = hinfo.n;
63: r = hinfo.r;
64: dt = hinfo.dt;
65: totalNumSim = hinfo.totalNumSim;
66: PetscMalloc1(2*n+1,&hinfo.vol);
67: vol = hinfo.vol;
68: St0 = hinfo.St0 = hinfo.vol + n;
69: readData(PETSC_COMM_WORLD,&hinfo);
71: numdim = n*(n+1)/2;
72: if (numdim%2 == 1) numdim++;
73: PetscMalloc1(numdim,&eps);
75: myNumSim = divWork(rank,totalNumSim,size);
77: x = 0;
78: for (i=0; i<myNumSim; i++) {
79: stdNormalArray(eps,numdim,ran);
80: x += basketPayoff(vol,St0,n,r,dt,eps);
81: }
83: MPI_Reduce(&x, &totalx, 1, MPIU_REAL, MPIU_SUM,0,PETSC_COMM_WORLD);
84: /* payoff = exp(-r*dt*n)*(totalx/totalNumSim);
85: 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",
86: payoff,(int)(stop - start),"parallel",size,"processors",n,(int)(1/dt),r); */
88: PetscFree(vol);
89: PetscFree(eps);
90: PetscRandomDestroy(&ran);
91: PetscFinalize();
92: return ierr;
93: }
95: PetscErrorCode stdNormalArray(PetscReal *eps, PetscInt numdim, PetscRandom ran)
96: {
97: PetscInt i;
98: PetscScalar u1,u2;
99: PetscReal t;
103: for (i=0; i<numdim; i+=2) {
104: PetscRandomGetValue(ran,&u1);
105: PetscRandomGetValue(ran,&u2);
107: t = PetscSqrtReal(-2*PetscLogReal(PetscRealPart(u1)));
108: eps[i] = t * PetscCosReal(2*PETSC_PI*PetscRealPart(u2));
109: eps[i+1] = t * PetscSinReal(2*PETSC_PI*PetscRealPart(u2));
110: }
111: return(0);
112: }
115: PetscReal basketPayoff(PetscReal vol[], PetscReal St0[], PetscInt n, PetscReal r,PetscReal dt, PetscReal eps[])
116: {
117: PetscReal Stk[PETSC_MAXBSIZE], temp;
118: PetscReal payoff;
119: PetscInt maxk,i,j;
120: PetscInt pointcount=0;
122: for (i=0;i<n;i++) Stk[i] = St0[i];
124: for (i=0;i<n;i++) {
125: maxk = 0;
126: for (j=0;j<(n-i);j++) {
127: Stk[j] = mcVal(Stk[j],r,vol[j],dt,eps[pointcount++]);
128: if ((Stk[j]/St0[j]) > (Stk[maxk]/St0[maxk])) maxk = j;
129: }
130: exchangeVal(Stk+j-1,Stk+maxk);
131: exchangeVal(St0+j-1,St0+maxk);
132: exchangeVal(vol+j-1,vol+maxk);
133: }
135: payoff = 0;
136: for (i=0; i<n; i++) {
137: temp = (Stk[i]/St0[i]) - 1;
138: if (temp > 0) payoff += temp;
139: }
140: return payoff;
141: }
143: PetscErrorCode readData(MPI_Comm comm,himaInfo *hinfo)
144: {
145: PetscInt i;
146: FILE *fd;
147: char temp[50];
149: PetscMPIInt rank;
150: PetscReal *v = hinfo->vol, *t = hinfo->St0;
151: PetscInt num=hinfo->n;
154: MPI_Comm_rank(comm,&rank);
155: if (!rank) {
156: PetscFOpen(PETSC_COMM_SELF,DATAFILENAME,"r",&fd);
157: for (i=0;i<num;i++) {
158: double vv,tt;
159: if (fscanf(fd,"%s%lf%lf",temp,&vv,&tt) != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file\n");
160: v[i] = vv;
161: t[i] = tt;
162: }
163: fclose(fd);
164: }
165: MPI_Bcast(v,2*num,MPIU_REAL,0,PETSC_COMM_WORLD);
166: /* 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]); */
167: return(0);
168: }
170: void exchangeVal(PetscReal *a, PetscReal *b)
171: {
172: PetscReal t;
174: t = *a;
175: *a = *b;
176: *b = t;
177: }
179: PetscReal mcVal(PetscReal St, PetscReal r, PetscReal vol, PetscReal dt, PetscReal eps)
180: {
181: return (St * PetscExpReal((r-0.5*vol*vol)*dt + vol*PetscSqrtReal(dt)*eps));
182: }
184: PetscInt divWork(PetscMPIInt id, PetscInt num, PetscMPIInt size)
185: {
186: PetscInt numit;
188: numit = (PetscInt)(((PetscReal)num)/size);
189: numit++;
190: return numit;
191: }
193: /*TEST
195: build:
196: requires: !comple
197: output_file: output/ex1_1.out
199: test:
200: nsize: 2
201: output_file: output/ex1_1.out
202: localrunfiles: ex2_stock.txt
204: TEST*/