Archive

Posts Tagged ‘monte carlo simulation’

Multithreading a Monte Carlo simulation

July 11, 2012 Leave a comment

One of the more nebulous areas of computer science is concurrency, with multi-threaded programming being particularly difficult to understand if you haven’t had the benefit of lots of experience. Of course, the precise degree to which your multi-threaded code is “difficult” depends hugely on the problem that the code is intended to solve (multi-threaded graphical user interfaces, for instance, are notoriously difficult to get working). There are, however, a number of problem domains that are eminently suited to multi-threading.

One of these is the problem of Monte Carlo simulation. I had a discussion about this yesterday with someone who is in the process of learning mathematical finance. He loves the subject but is frustrated by the fact that he’s essentially using only one of the four cores on his machine whenever he writes code to, say, price an option using a Monte Carlo algorithm. Surely, he said, there must be an easy way to get a multi-threaded Monte Carlo routing going?

Given that he’s been reading Mark Joshi’s excellent books on mathematical finance and computational finance, I thought that I’d run through the simplest possible example of how multi-threading his routine would result in faster run times. Here’s the code I knocked together in a few minutes based on a trivial Monte Carlo pricer for a call option. (Note: this was for demonstration purposes only, intended to illustrate the ideas behind threading and the Boost Thread library in particular; writing code such as this in a bank will get you fired at your first code review!)

#include <iostream>
using namespace std;

#include <boost/thread.hpp>
#include <boost/timer.hpp>

#if defined(_USE_BOOST_MATH)
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/normal_distribution.hpp>
#include <boost/random/variate_generator.hpp>

typedef boost::mt19937                                      T_base_prng;
typedef boost::normal_distribution<>                        T_norm_dist;
typedef boost::variate_generator<T_base_prng&, T_norm_dist> T_norm_varg;
#endif

#if !defined(_MSC_VER)
#define EXIT_SUCCESS 0
#define EXIT_FAILURE 1
#endif

double getOneGaussianByBoxMuller();

double SimpleMC(double expiry,
                double strike,
                double spot,
                double rate,
                double vol,
                unsigned long numPaths,
                unsigned long base_seed = 42);

void SimpleMC_threaded(double expiry,
                       double strike,
                       double spot,
                       double rate,
                       double vol,
                       unsigned long numPaths,
                       double* target,
                       unsigned long base_seed = 42);

int main(int argc, char* argv[])
{
    try {

        double expiry          = 1;
        double strike          = 100;
        double spot            = 100;
        double rate            = 0.05;
        double vol             = 0.20;
        unsigned long numPaths = 1e8;

        boost::timer timer;
        double result_unthreaded = SimpleMC(expiry, strike, spot, rate,
                                            vol, numPaths);
        double time_unthreaded = timer.elapsed();

        cout << "Single-threaded run\n"
             << "-------------------\n"
             << "Call: " << result_unthreaded << endl
             << "Time: " << time_unthreaded
             << "\n\n";

        double result1;
        double result2;
        double result3;
        double result4;
        boost::timer timer_threaded;

        boost::thread worker1(SimpleMC_threaded, expiry, strike, spot,
                              rate, vol, numPaths / 4, &result1, 41);
        boost::thread worker2(SimpleMC_threaded, expiry, strike, spot,
                              rate, vol, numPaths / 4, &result2, 42);
        boost::thread worker3(SimpleMC_threaded, expiry, strike, spot,
                              rate, vol, numPaths / 4, &result3, 43);
        boost::thread worker4(SimpleMC_threaded, expiry, strike, spot,
                              rate, vol, numPaths / 4, &result4, 44);
        worker1.join();
        worker2.join();
        worker3.join();
        worker4.join();

        double result_threaded = (result1 + result2 + result3 + result4) / 4;
        double time_threaded   = timer_threaded.elapsed();

        cout << "Multi-threaded run\n"
             << "------------------\n"
             << "Call: " << result_threaded << endl
             << "Time: " << time_threaded
             << "\n\n";

#if !defined(_MSC_VER)
        cout << "Press Enter to continue.";
        cin.get();
#else
        system("PAUSE");
#endif

        return (EXIT_SUCCESS);
    }
    catch (exception& e) {
        cerr << e.what() << endl;
        return (EXIT_FAILURE);
    }
    catch (...) {
        cerr << "Non-standard exception caught... " << endl;
        return (EXIT_FAILURE);
    }
}

double SimpleMC(double expiry,
                double strike,
                double spot,
                double rate,
                double vol,
                unsigned long numPaths,
                unsigned long base_seed)
{
    double variance      = vol * vol * expiry;
    double rootVariance  = sqrt(variance);
    double itoCorrection = -0.5 * variance;
    double movedSpot     = spot * exp(rate * expiry + itoCorrection);
    double runningSum    = 0;
    double thisSpot;

#if !defined(_USE_BOOST_MATH)
    for (unsigned long i = 0; i < numPaths; ++i)
    {
        double thisGaussian = getOneGaussianByBoxMuller();
        thisSpot = movedSpot * exp(rootVariance * thisGaussian);
        double thisPayoff = thisSpot - strike;
        thisPayoff = thisPayoff > 0 ? thisPayoff : 0;
        runningSum += thisPayoff;
    }
#else
    T_base_prng base_prng(base_seed);
    T_norm_dist norm_dist(0, 1);
    T_norm_varg norm_varg(base_prng, norm_dist);

    for (unsigned long i = 0; i < numPaths; ++i) 
    {
        thisSpot = movedSpot * exp(rootVariance * norm_varg());
        double thisPayoff = thisSpot - strike;
        thisPayoff = thisPayoff > 0 ? thisPayoff : 0;
        runningSum += thisPayoff;
    }
#endif

    double mean = runningSum / numPaths;
    mean *= exp(-rate * expiry);

    return mean;
}

void SimpleMC_threaded(double expiry,
                       double strike,
                       double spot,
                       double rate,
                       double vol,
                       unsigned long numPaths,
                       double* target,
                       unsigned long base_seed)
{
    double variance      = vol * vol * expiry;
    double rootVariance  = sqrt(variance);
    double itoCorrection = -0.5 * variance;
    double movedSpot     = spot * exp(rate * expiry + itoCorrection);
    double runningSum    = 0;
    double thisSpot;

#if !defined(_USE_BOOST_MATH)
    for (unsigned long i = 0; i < numPaths; ++i)
    {
        double thisGaussian = getOneGaussianByBoxMuller();
        thisSpot = movedSpot * exp(rootVariance * thisGaussian);
        double thisPayoff = thisSpot - strike;
        thisPayoff = thisPayoff > 0 ? thisPayoff : 0;
        runningSum += thisPayoff;
    }
#else
    T_base_prng base_prng(base_seed);
    T_norm_dist norm_dist(0, 1);
    T_norm_varg norm_varg(base_prng, norm_dist);

    for (unsigned long i = 0; i < numPaths; ++i) 
    {
        thisSpot = movedSpot * exp(rootVariance * norm_varg());
        double thisPayoff = thisSpot - strike;
        thisPayoff = thisPayoff > 0 ? thisPayoff : 0;
        runningSum += thisPayoff;
    }
#endif

    double mean = runningSum / numPaths;
    mean *= exp(-rate * expiry);

    *target = mean;
}

double getOneGaussianByBoxMuller()
{
    double result;
    double x;
    double y;
    double sizeSquared;

    do {
        x = 2.0 * rand() / static_cast<double>(RAND_MAX) - 1;
        y = 2.0 * rand() / static_cast<double>(RAND_MAX) - 1;
        sizeSquared = x * x + y * y;
    }
    while (sizeSquared >= 1.0);

    result = x * sqrt(-2 * log(sizeSquared) / sizeSquared);

    return result;
}