1 static char help[] = "Tests PetscRandom functions.\n\n";
2
3 #include <petscsys.h>
4
5 #define PETSC_MAXBSIZE 40
6 #define DATAFILENAME "ex2_stock.txt"
7
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;
17
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);
24
25 /*
26 Contributed by Xiaoyan Zeng <zengxia@iit.edu> and Liu, Kwong Ip" <kiliu@math.hkbu.edu.hk>
27
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 */
31
main(int argc,char * argv[])32 int main(int argc, char *argv[])
33 {
34 PetscReal r, dt;
35 PetscInt n;
36 PetscInt i, myNumSim, totalNumSim, numdim;
37 PetscReal *vol, *St0, x, totalx;
38 PetscMPIInt size, rank;
39 PetscReal *eps;
40 himaInfo hinfo;
41 PetscRandom ran;
42
43 PetscFunctionBeginUser;
44 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
45 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &ran));
46 PetscCall(PetscRandomSetFromOptions(ran));
47
48 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
49 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
50
51 hinfo.n = 31;
52 hinfo.r = 0.04;
53 hinfo.dt = 1.0 / 12; /* a month as a period */
54 hinfo.totalNumSim = 1000;
55
56 PetscCall(PetscOptionsGetInt(NULL, NULL, "-num_of_stocks", &hinfo.n, NULL));
57 PetscCheck(hinfo.n >= 1 && hinfo.n <= 31, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only 31 stocks listed in stock.txt. num_of_stocks %" PetscInt_FMT " must between 1 and 31", hinfo.n);
58 PetscCall(PetscOptionsGetReal(NULL, NULL, "-interest_rate", &hinfo.r, NULL));
59 PetscCall(PetscOptionsGetReal(NULL, NULL, "-time_interval", &hinfo.dt, NULL));
60 PetscCall(PetscOptionsGetInt(NULL, NULL, "-num_of_simulations", &hinfo.totalNumSim, NULL));
61
62 n = hinfo.n;
63 r = hinfo.r;
64 dt = hinfo.dt;
65 totalNumSim = hinfo.totalNumSim;
66 PetscCall(PetscMalloc1(2 * n + 1, &hinfo.vol));
67 vol = hinfo.vol;
68 St0 = hinfo.St0 = hinfo.vol + n;
69 PetscCall(readData(PETSC_COMM_WORLD, &hinfo));
70
71 numdim = n * (n + 1) / 2;
72 if (numdim % 2 == 1) numdim++;
73 PetscCall(PetscMalloc1(numdim, &eps));
74
75 myNumSim = divWork(rank, totalNumSim, size);
76
77 x = 0;
78 for (i = 0; i < myNumSim; i++) {
79 PetscCall(stdNormalArray(eps, numdim, ran));
80 x += basketPayoff(vol, St0, n, r, dt, eps);
81 }
82
83 PetscCallMPI(MPI_Reduce(&x, &totalx, 1, MPIU_REAL, MPIU_SUM, 0, PETSC_COMM_WORLD));
84 /* payoff = exp(-r*dt*n)*(totalx/totalNumSim);
85 ierr = 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); */
87
88 PetscCall(PetscFree(vol));
89 PetscCall(PetscFree(eps));
90 PetscCall(PetscRandomDestroy(&ran));
91 PetscCall(PetscFinalize());
92 return 0;
93 }
94
stdNormalArray(PetscReal * eps,PetscInt numdim,PetscRandom ran)95 PetscErrorCode stdNormalArray(PetscReal *eps, PetscInt numdim, PetscRandom ran)
96 {
97 PetscInt i;
98 PetscScalar u1, u2;
99 PetscReal t;
100
101 PetscFunctionBegin;
102 for (i = 0; i < numdim; i += 2) {
103 PetscCall(PetscRandomGetValue(ran, &u1));
104 PetscCall(PetscRandomGetValue(ran, &u2));
105
106 t = PetscSqrtReal(-2 * PetscLogReal(PetscRealPart(u1)));
107 eps[i] = t * PetscCosReal(2 * PETSC_PI * PetscRealPart(u2));
108 eps[i + 1] = t * PetscSinReal(2 * PETSC_PI * PetscRealPart(u2));
109 }
110 PetscFunctionReturn(PETSC_SUCCESS);
111 }
112
basketPayoff(PetscReal vol[],PetscReal St0[],PetscInt n,PetscReal r,PetscReal dt,PetscReal eps[])113 PetscReal basketPayoff(PetscReal vol[], PetscReal St0[], PetscInt n, PetscReal r, PetscReal dt, PetscReal eps[])
114 {
115 PetscReal Stk[PETSC_MAXBSIZE], temp;
116 PetscReal payoff;
117 PetscInt maxk, i, j;
118 PetscInt pointcount = 0;
119
120 for (i = 0; i < n; i++) Stk[i] = St0[i];
121
122 for (i = 0; i < n; i++) {
123 maxk = 0;
124 for (j = 0; j < (n - i); j++) {
125 Stk[j] = mcVal(Stk[j], r, vol[j], dt, eps[pointcount++]);
126 if ((Stk[j] / St0[j]) > (Stk[maxk] / St0[maxk])) maxk = j;
127 }
128 exchangeVal(Stk + j - 1, Stk + maxk);
129 exchangeVal(St0 + j - 1, St0 + maxk);
130 exchangeVal(vol + j - 1, vol + maxk);
131 }
132
133 payoff = 0;
134 for (i = 0; i < n; i++) {
135 temp = (Stk[i] / St0[i]) - 1;
136 if (temp > 0) payoff += temp;
137 }
138 return payoff;
139 }
140
readData(MPI_Comm comm,himaInfo * hinfo)141 PetscErrorCode readData(MPI_Comm comm, himaInfo *hinfo)
142 {
143 PetscInt i;
144 FILE *fd;
145 char temp[50];
146 PetscMPIInt rank;
147 PetscReal *v = hinfo->vol, *t = hinfo->St0;
148 PetscInt num = hinfo->n;
149
150 PetscFunctionBeginUser;
151 PetscCallMPI(MPI_Comm_rank(comm, &rank));
152 if (rank == 0) {
153 PetscCall(PetscFOpen(PETSC_COMM_SELF, DATAFILENAME, "r", &fd));
154 for (i = 0; i < num; i++) {
155 double vv, tt;
156 PetscCheck(fscanf(fd, "%s%lf%lf", temp, &vv, &tt) == 3, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Badly formatted input file");
157 v[i] = vv;
158 t[i] = tt;
159 }
160 fclose(fd);
161 }
162 PetscCallMPI(MPI_Bcast(v, (PetscMPIInt)(2 * num), MPIU_REAL, 0, PETSC_COMM_WORLD));
163 /* ierr = 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]); */
164 PetscFunctionReturn(PETSC_SUCCESS);
165 }
166
exchangeVal(PetscReal * a,PetscReal * b)167 void exchangeVal(PetscReal *a, PetscReal *b)
168 {
169 PetscReal t;
170
171 t = *a;
172 *a = *b;
173 *b = t;
174 }
175
mcVal(PetscReal St,PetscReal r,PetscReal vol,PetscReal dt,PetscReal eps)176 PetscReal mcVal(PetscReal St, PetscReal r, PetscReal vol, PetscReal dt, PetscReal eps)
177 {
178 return St * PetscExpReal((r - 0.5 * vol * vol) * dt + vol * PetscSqrtReal(dt) * eps);
179 }
180
divWork(PetscMPIInt id,PetscInt num,PetscMPIInt size)181 PetscInt divWork(PetscMPIInt id, PetscInt num, PetscMPIInt size)
182 {
183 PetscInt numit;
184
185 numit = (PetscInt)(((PetscReal)num) / size);
186 numit++;
187 return numit;
188 }
189
190 /*TEST
191
192 test:
193 nsize: 2
194 output_file: output/empty.out
195 localrunfiles: ex2_stock.txt
196
197 TEST*/
198