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

Short straws

September 11, 2010 Leave a comment

Phew! Just my luck: I finally decide to commit to starting a blog when all hell breaks loose with work and I have to pull eighteen-hour days just to keep on top of things. What’s even more annoying is that I actually have a ton of posts in preparation but I’ve discovered that – just like when I’m trying to write a paper – I’m never happy with the results and keep trying to hone and tweak them to an inevitably unattainable level of perfection.

The moral of all this is of course that chasing perfection is a great idea when you’re trying to get a paper published, but it’s an absolutely terrible idea when you’re trying to get a blog off the ground. I shall have to pay attention to this in future.

Anyway, while we’re waiting for the good stuff to come along, I thought it might be worthwhile mentioning a conversation I had at lunch yesterday. I got into a discussion with a colleague about finance, a discussion which took rather a surprising turn. He’s someone whose knowledge on all things physics related genuinely impresses me but he hit me with somewhat of a bombshell half way through our conversation. We were talking about how large – sometimes stupidly large – amounts of money are made each day on stock exchanges by trading in shares. He understood how someone could make money by buying a stock and holding it in the expectation that its price could rise in value. But – and this is a big but – he had no idea how short selling works.

I’ve talked to some other people since then and almost all of them seem to have the same difficulty. I admit, short-selling stocks is a weird idea when you first encounter it, but it really couldn’t be simpler. Here’s how I explained it.

Suppose your next-door neighbour loves video games. In particular, he’s very fond of his Playstation 3 and tries to keep up with all of the latest game releases. Several months ago, he went to Amazon and pre-ordered a copy of an eagerly awaited game so that he’d be sure of having it delivered on its day of release. He waits and waits until, eventually, the release date arrives, his credit card is charged £40, and the postman pops a copy of it through his letterbox. However, being a busy man, he’s due to go on a business trip later that day and he won’t be able to play the game until he gets back in, say, a month’s time.

You’ve been watching his growing excitement about the release of the game with interest. You think to yourself, “Hey, if he’s this excited about the game surely others are too. And if others are excited about it, perhaps there’s a chance to make a few bob here.” (You say this to yourself in your best Del Boy voice.) Knowing that he’s going away for a while – and knowing that the game is in such high demand that there are actually shortages in stores – you ask him if he would mind you borrowing it while he’s gone, promising to give it back to him the moment he arrives back from his trip. In return for allowing you to borrow the game, you might decide to thank him by promising to keep an eye on his house while he’s gone.

Notice that you promised only to give it back to him when he comes home again; in particular, you didn’t actually tell him what you were going to do with the game! Instead of sitting down to play it, you decide to put it to work for you. You go to ebay, say, and place the copy of the game for sale there. You find to your delight that such is the enormous demand for the game, people are willing to pay a huge premium in order to get their hands on it. The ebay auction actually ends with one poor sucker bidder buying the game from you for £80 – double what your friend paid! You smile as your PayPal account grows by £80 and pop the game in the post to the buyer.

A month passes and your neighbour is due back in a few hours. Remembering that you promised to give him back a copy of the game, you stroll down to your local video games store and find – surprise! – that plenty of copies are now available. Even better, they’re available for only £30 now that the initial excitement has worn off. You pick up a copy and on the way home you buy a couple of six packs of beer, one of which is to thank your friend for letting you have the game.

That evening, you smile to yourself as you enjoy your beer.

Contrived though this example is, it captures essentially everything you need to know about short selling stocks. Broadly speaking, a trader will short sell a stock if he or she believes that the stock is over valued or that the stock price will drop over a certain time frame. Although different regulatory requirements mean that the details can sometimes change, the basic procedure is always the same: the trader seeks out someone who already has the stock and asks to borrow it from them for a fixed period of time in return for a small fee. He then sells the stock on the open market and puts the cash in his account. When the time comes for him to return the stock to the person who lent it to him, the trader buys the stock from the market, hoping that its price has gone down in the meantime. If it has, he makes a profit; if it hasn’t, he has to bear the loss involved with buying the stock and giving it back to its original owner.

There are, of course, plenty of reasons why shorting stocks in practice isn’t as simple as this (regulatory requirements, temporary freezes on shorts, managing the dividend payments on stocks you’ve shorted, stock splits that happen during a short, and so forth). Shorting also has rather a bad image among the public. But it’s really quite a straightforward – and ultimately beneficial – practice.

And it’s all done so that you can enjoy that beer when your short comes home!

Categories: Finance

Loose language

August 10, 2010 3 comments

Most of my work in physics is done from the perspective of a theoretical physicist. For some reason or other, when I was younger I found the clarity and precision of theoretical physics to be much more attractive than the rough and ready approach often found in experimental physics.  Not only was I afflicted with the unfortunate ability to render an experiment senseless simply by being in the same room as it, I also found experimental physics to be too subjective, too open to interpretation for my liking. Thus, at some point I decided to stick to theory and I’ve been there ever since.

As wonderful as I find this particular approach to physics, it’s not without its drawbacks. For one thing, in order to be a good theorist, you’ve got to have a good understanding not only of physics but also of the language in which a theorist speaks about physics. This invariably means that you’ve got to be a good mathematician. Early on, I started to think of myself as mediocre — at best — at mathematics. This lasted for several years until I figured out why I was mediocre: I was attempting to understand mathematics largely by reading papers written by mathematicians. Now, I love mathematicians as much as the next man but, by God, they’re awful at writing in a manner that’s clear to a non-expert. At some point during the last century, the Bourbaki crowd decided it would be a great idea to present mathematics in as dense and impenetrable a style as possible; this caught on and we ended up being left with several generations of mathematicians who believed that it’s fine to write about mathematics in a manner that makes it almost impossible for students to tell what’s going on. I’m all for rigour when presenting and proving results, but rigour is generally of little use in a pedagogical sense.

Of course, one can also go too far the other way. Physicists, for instance, tend to speak about mathematics in a manner that horrifies most mathematicians. In addition, many of the “intuitive” ways in which physicists present ideas really aren’t intuitive at all since they often don’t stand up to even the merest scrutiny or immediately run into difficulties once one attempts to extend intuitive pictures beyond simple cases. Perhaps the canonical examples of this sort of problem are to be found in differential geometry, a field which you really must understand if you want to consider yourself a decent theoretical physicist.

One of the things that used to irritate me considerably when I first attempted to learn differential geometry was the fact that I often found the mathematician’s perspective on the subject to be too abstract and needlessly encumbered by Bourbakization. On the other hand, physicists tend to speak about differential geometry using analogies that are often too loose to be of much help. A good example of this lies in the example of a tangent space to a differentiable manifold. If you’re taught the subject by a physicist or happen to learn about it from one of the many books written by physicists, you’ll often first be presented with a pictorial representation of what a physicist thinks a tangent space “looks like”. For instance, here’s a commonly encountered way of representing the tangent space to S^2, a two-sphere:

This is a nice picture to keep in your head if you find such things helpful. However, it’s not at all representative of what a tangent space actually is. For one thing, representing the tangent space to a point X\in S^2 like this really presupposes that you have embedded S^2 in some higher-dimensional flat space; the tangent space at X is then regarded as a hyperplane in this space. However, differential geometry is largely concerned with the intrinsic geometry of manifolds, that is, the study of their geometry without assuming that they are embedded in a higher-dimensional manifold. The picture above makes sense only if we are talking about the extrinsic geometry; we regard S^2 as being embedded in, say, \mathbb{R}^3 and consider the tangent space at X to be a two-dimensional hyperplane since this is the only real way in which the linearity of the tangent space (it is, after all, just some sort of vector space) can be related to the “flatness” of the plane in the picture.

A further problem with the way physicists speak about differential geometry is that the notation one uses when performing practical calculations isn’t necessarily the best notation to use when teaching the subject. To get an idea of what I mean by this, it’s worth bearing in mind that physicists tend as a rule to play fast and loose with notation; mathematicians, on the other hand, tend to be so precise about their notation that it sometimes induces seizures. For instance, if you check out the book I linked to earlier, you’ll find that the author freely admits that he constantly abuses notation when introducing the concept of a tangent vector as a directional derivative of a function along a curve in a manifold. Typically, this means that physicists often begin by saying that f:M\to\mathbb{R} is a real-valued map on a manifold M but later, when it comes to actually performing interesting calculations, they’ll happily use f to denote what is actually the coordinate presentation f\circ\varphi^{-1}_U of the map in some chart (U,\varphi_U). This irritates me no end and causes huge problems for students when they’re initially learning this stuff.

One of the unhappy consequences of looking at differential geometry through the lens of these loose analogies is that it sometimes convinces people to think that the best way to introduce a topic is to rely heavily on what they consider to be the intuitive aspects. This, however, is often a mistake. The concept of a tangent space, for instance, is a great example of how leaning too heavily on intuitive pictures can actually result in you overlooking the much more important algebraic structure that underlies everything. Let me say something slightly controversial:

Think about a tangent space first and foremost as something algebraic in nature. More specifically, tangent spaces are what you find when you play around with the linear spaces of germs of smooth functions in a neighbourhood of a point p\in M.

Once you’ve got the algebraic nature of tangent spaces under your belt, you can then consider tangent spaces as collections of equivalence classes of smooth curves through p. You’ll then understand not only the geometrical picture of what a tangent space is, but also the algebraic reasons that underlie why it is. I’ll deal with this more in a future post.

Physics, but from a hedgerow

August 9, 2010 Leave a comment

Welcome one and all to physics from a hedgerow. Sure, the name is unwieldy but it has a significance I might explain one day.

For the moment I’m still trying to get used to wordpress so don’t expect a great deal of activity over the coming week. Beyond that, however, we’ll be taking some — hopefully interesting — trips into physics, artificial intelligence, and money. Mostly physics of course but rather a lot about the others too.

Indeed.