Actual source code: ex2.c
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: PetscInt i, myNumSim, totalNumSim, numdim;
37: PetscReal *vol, *St0, x, totalx;
38: PetscMPIInt size, rank;
39: PetscReal *eps;
40: himaInfo hinfo;
41: PetscRandom ran;
43: PetscFunctionBeginUser;
44: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
45: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &ran));
46: PetscCall(PetscRandomSetFromOptions(ran));
48: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
49: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
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: 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));
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));
71: numdim = n * (n + 1) / 2;
72: if (numdim % 2 == 1) numdim++;
73: PetscCall(PetscMalloc1(numdim, &eps));
75: myNumSim = divWork(rank, totalNumSim, size);
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: }
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); */
88: PetscCall(PetscFree(vol));
89: PetscCall(PetscFree(eps));
90: PetscCall(PetscRandomDestroy(&ran));
91: PetscCall(PetscFinalize());
92: return 0;
93: }
95: PetscErrorCode stdNormalArray(PetscReal *eps, PetscInt numdim, PetscRandom ran)
96: {
97: PetscInt i;
98: PetscScalar u1, u2;
99: PetscReal t;
101: PetscFunctionBegin;
102: for (i = 0; i < numdim; i += 2) {
103: PetscCall(PetscRandomGetValue(ran, &u1));
104: PetscCall(PetscRandomGetValue(ran, &u2));
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: }
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;
120: for (i = 0; i < n; i++) Stk[i] = St0[i];
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: }
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: }
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;
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: }
167: void exchangeVal(PetscReal *a, PetscReal *b)
168: {
169: PetscReal t;
171: t = *a;
172: *a = *b;
173: *b = t;
174: }
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: }
181: PetscInt divWork(PetscMPIInt id, PetscInt num, PetscMPIInt size)
182: {
183: PetscInt numit;
185: numit = (PetscInt)(((PetscReal)num) / size);
186: numit++;
187: return numit;
188: }
190: /*TEST
192: test:
193: nsize: 2
194: output_file: output/ex1_1.out
195: localrunfiles: ex2_stock.txt
197: TEST*/