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*/