Archive

Archive for the ‘Programming’ Category

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;
}

Are computers causing us to produce stupid physicists?

September 27, 2010 Leave a comment

An interesting post over at the Wolfram blog asks whether the use of computer software such as Mathematica contributes to the (perceived) decline in standards in the mathematical abilities of students. It’s an interesting question, and one that has been asked for almost as long as programs such as Mathematica have been used in educating new students.

On the face of it, there’s much to recommend the question. After all, we all know older mathematicians or physicists who can apparently solve hideously complicated problems in their heads. Since this is a skill that younger mathematicians/physicists appear to lack, is it not reasonable to suspect that computers play some part in this?

Well, no. The post I linked to reflects my views on the subject quite closely and is well worth a read if this sort of thing is of interest to you.

Cute stuff

September 27, 2010 Leave a comment

Although most of the programming I do is geared towards the purely numerical end of the spectrum, there are times when I find myself faced with a task that would appear to benefit from having a handy GUI. Often, MATLAB provides me with sufficient tools for this since it comes with a built-in GUI builder that you can use to construct graphical interfaces to your code. However, this suffers from two drawbacks:

1) If MATLAB isn’t great at tackling the underlying problem you want to slap a GUI on, the MATLAB GUI builder probably won’t be a great fit for your needs. For instance, if you’d like to build, say, a graphical application that presents the results of monitoring some simulation that you’re running remotely, MATLAB won’t be a good choice since its networking features are, frankly, terrible.

2) MATLAB’s GUI builder uses Java. Enough said.

Having found myself faced with this kind of problem more times than I care to remember, at one point I decided to do something about it .There are several excellent GUI libraries out there that allow you to build stable and functional graphical interfaces to your code but for me there is one in particular that stands out: Qt. Not only is it a breeze to learn if you’re already familiar with C++, it’s also cross-platform, meaning that running your application on different operating systems is (usually, at least) a case of a simple rebuild. It also allows me to hook into existing C++ code very easily, which is great considering how much of the stuff I’ve got lying around.

If you’re interested in learning a GUI framework for similar reasons to the ones I’ve outlined above, Qt may be a good fit for you. The Qt documentation is great and there are plenty of tutorials online. There’s also a rather good book available on the subject, with a more advanced one due out soon. You can even read up on how to build your Qt applications with Python if that’s more your sort of thing. Indeed, the first of these books has been so useful to me when learning Qt that I’m going to shamelessly steal from it during the next post, where we’ll show just how easy it is to get a Qt application up and running. We’ll also cover one or two of the subtleties you need to take into account if you’re using a new version of Qt since the book deals only with an earlier version of Qt 4.

Categories: C++, MATLAB, Programming, Qt