bad C++ code for fast monte carlo options pricing.

main.cpp

/*
 * ### bad C++ code for fast monte carlo options pricing.
 *
 * ### references
 *
 *  https://news.ycombinator.com/item?id=7450601
 *  http://rawrjustin.github.io/blog/2014/03/18/julia-vs-python-monte-carlo-simulations-of-bitcoin-options/
 *  code originally modelled on https://github.com/arrayfire/arrayfire_examples/blob/master/financial/monte_carlo_options.cpp
 *  the examples of using omp for parallel random number generation from http://numbercrunch.de/trng/
 *
 *
 * ### to compile
 *
 *  g++ -fwhole-program -Wall -Werror -pedantic -O3 -std=c++0x -march=native main.cpp -fopenmp
 *
 * ### to generate random numbers for 100,000 paths and write them to 'randn.dat':
 *
 *  ./a.out g n100000
 *
 * ### to simulate 100,000 paths using prepared random numbers from disk
 *
 *  ./a.out n100000
 *
 */


#include <stdio.h>
#include <string.h>
#include <errno.h>
#include <sys/mman.h>

#include <cmath>

#include <omp.h>

#include <boost/random.hpp>
#include <boost/random/normal_distribution.hpp>


double monte_carlo(int N, double s_0, double K, double t, double sigma, double r, int steps, double *randn) {
    double dt = t / (double)(steps);
    double a = (r - (sigma * sigma * 0.5)) * dt;
    double b = sigma * sqrt(dt);

    int samples = N * steps;

    double acc = 0.0;
#pragma omp parallel reduction(+:acc)
{
    int size = omp_get_num_threads();
    int rank = omp_get_thread_num();
    printf("-- hello from thread %d of %d\n", rank, size);

    int k = (rank*samples/size);

    int i0 = (rank*N/size);
    int i1 = (rank+1)*N/size;

    for (int i=i0; i!=i1; ++i) {
        double log_s = log(s_0);
        for (int j=0; j!=steps; ++j) {
            double noise = randn[k++];
            log_s += (a + b * noise);
        }
        double s = exp(log_s);
        double sk_pos = (s-K > 0.0) ? (s-K) : 0.0; // clamp to be positive
        double result = sk_pos;
        acc += result;
    }
}
    return (acc / N) * exp(-r * t);
}


void generate_random_data(const char * file_name, int samples) {

    boost::mt19937 rng;
    boost::normal_distribution<> nd(0.0, 1.0);
    boost::variate_generator<boost::mt19937&, boost::normal_distribution<> > R(rng, nd);

    double * buffer = new double[samples];
#pragma omp parallel
{
    int size = omp_get_num_threads();
    int rank = omp_get_thread_num();
    printf("-- hello from thread %d of %d\n", rank, size);
    int i0 = (rank*samples/size);
    int i1 = ((rank+1)*samples/size);
    for (int i=i0; i!=i1; ++i) {
        buffer[i] = R();
    }
}

    FILE *f = fopen(file_name, "w");
    fwrite(buffer, sizeof(double), samples, f);
    fclose(f);
    delete[] buffer;
}


int main(int narg, char **argv) {
    // -- bad parsing code
    int N = 100000;
    bool generate = false;
    for (int arg=1; arg!=narg; ++arg) {
        switch (argv[arg][0]) {
            case 'g':
                generate = true;
                break;
            case 'n':
                N = atoi(argv[arg]+1);
                printf("setting N=%d\n", N);
                break;
            default:
                printf("confused by argument: '%s'\n", argv[arg]);
                printf("usage: %s [g] [nN]\n", argv[0]);
                printf("\twhere optional g argument used to generate random numbers\n");
                printf("\tand optional nN argument sets number of paths to sample, e.g. n100 for 100 paths\n");
                return 1;
        }
    }

    const char* data_file_name = "randn.dat";
    int steps = 100;
    int samples = N * steps;

    if (generate) {
        printf("generating %d random samples, saving to '%s'\n", samples, data_file_name);
        generate_random_data(data_file_name, samples);
        printf("ok");
    } else {
        printf("reading samples from '%s'\n", data_file_name);
        FILE *f = fopen(data_file_name, "r");
        if (f == NULL)
            goto fail;
        int fd = fileno(f);
        if (fd == -1)
            goto fail;
        void *p = mmap(NULL, samples * sizeof(double), PROT_READ, MAP_PRIVATE, fd, 0);
        if (p == MAP_FAILED)
            goto fail;

        double *randn = (double*)p;

        printf("simulating...\n");

        double s_0 = 600.0; // current thing price
        double stock_price = 1000.0;
        double t = 1.0;
        double sigma = 2.0;
        double r = 0.02;

        int n_reps = 10;
        for (int i = 0; i < n_reps; ++i) {
            printf("simulating, repetition %d of %d\n", i, n_reps);
            printf("result %g\n", monte_carlo(N, s_0, stock_price, t, sigma, r, steps, randn));
        }

        printf("unloading samples\n");

        int rv;
        rv = munmap(p, samples * sizeof(double));
        if (rv == -1)
            goto fail;
        rv = fclose(f);
        if (rv == EOF)
            goto fail;

        printf("fin\n");
    }
    return 0;
fail:
    int err = errno;
    printf("something broke; errno %d, '%s'\n", err, strerror(err));
}