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:   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);
 45:   PetscRandomCreate(PETSC_COMM_WORLD, &ran);
 46:   PetscRandomSetFromOptions(ran);

 48:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 49:   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:   PetscOptionsGetInt(NULL, NULL, "-num_of_stocks", &(hinfo.n), NULL);
 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 0;
 93: }

 95: PetscErrorCode stdNormalArray(PetscReal *eps, PetscInt numdim, PetscRandom ran)
 96: {
 97:   PetscInt    i;
 98:   PetscScalar u1, u2;
 99:   PetscReal   t;

101:   for (i = 0; i < numdim; i += 2) {
102:     PetscRandomGetValue(ran, &u1);
103:     PetscRandomGetValue(ran, &u2);

105:     t          = PetscSqrtReal(-2 * PetscLogReal(PetscRealPart(u1)));
106:     eps[i]     = t * PetscCosReal(2 * PETSC_PI * PetscRealPart(u2));
107:     eps[i + 1] = t * PetscSinReal(2 * PETSC_PI * PetscRealPart(u2));
108:   }
109:   return 0;
110: }

112: PetscReal basketPayoff(PetscReal vol[], PetscReal St0[], PetscInt n, PetscReal r, PetscReal dt, PetscReal eps[])
113: {
114:   PetscReal Stk[PETSC_MAXBSIZE], temp;
115:   PetscReal payoff;
116:   PetscInt  maxk, i, j;
117:   PetscInt  pointcount = 0;

119:   for (i = 0; i < n; i++) Stk[i] = St0[i];

121:   for (i = 0; i < n; i++) {
122:     maxk = 0;
123:     for (j = 0; j < (n - i); j++) {
124:       Stk[j] = mcVal(Stk[j], r, vol[j], dt, eps[pointcount++]);
125:       if ((Stk[j] / St0[j]) > (Stk[maxk] / St0[maxk])) maxk = j;
126:     }
127:     exchangeVal(Stk + j - 1, Stk + maxk);
128:     exchangeVal(St0 + j - 1, St0 + maxk);
129:     exchangeVal(vol + j - 1, vol + maxk);
130:   }

132:   payoff = 0;
133:   for (i = 0; i < n; i++) {
134:     temp = (Stk[i] / St0[i]) - 1;
135:     if (temp > 0) payoff += temp;
136:   }
137:   return payoff;
138: }

140: PetscErrorCode readData(MPI_Comm comm, himaInfo *hinfo)
141: {
142:   PetscInt    i;
143:   FILE       *fd;
144:   char        temp[50];
145:   PetscMPIInt rank;
146:   PetscReal  *v = hinfo->vol, *t = hinfo->St0;
147:   PetscInt    num = hinfo->n;

150:   MPI_Comm_rank(comm, &rank);
151:   if (rank == 0) {
152:     PetscFOpen(PETSC_COMM_SELF, DATAFILENAME, "r", &fd);
153:     for (i = 0; i < num; i++) {
154:       double vv, tt;
156:       v[i] = vv;
157:       t[i] = tt;
158:     }
159:     fclose(fd);
160:   }
161:   MPI_Bcast(v, 2 * num, MPIU_REAL, 0, PETSC_COMM_WORLD);
162:   /* 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]); */
163:   return 0;
164: }

166: void exchangeVal(PetscReal *a, PetscReal *b)
167: {
168:   PetscReal t;

170:   t  = *a;
171:   *a = *b;
172:   *b = t;
173: }

175: PetscReal mcVal(PetscReal St, PetscReal r, PetscReal vol, PetscReal dt, PetscReal eps)
176: {
177:   return (St * PetscExpReal((r - 0.5 * vol * vol) * dt + vol * PetscSqrtReal(dt) * eps));
178: }

180: PetscInt divWork(PetscMPIInt id, PetscInt num, PetscMPIInt size)
181: {
182:   PetscInt numit;

184:   numit = (PetscInt)(((PetscReal)num) / size);
185:   numit++;
186:   return numit;
187: }

189: /*TEST

191:    test:
192:       nsize: 2
193:       output_file: output/ex1_1.out
194:       localrunfiles: ex2_stock.txt

196: TEST*/