Monday, May 12, 2014

Chapter 6, part 1 of n: random-number generation

Welcome back.

Today's post is the first of a new series that covers chapter 6 of my book. The subject is the Monte Carlo framework of QuantLib, and it's going to take a few posts (probably interspersed with other content).

A plea for help: last month, Richard Stanton reported on the mailing list that compilation of the Python bindings fails on Mac OS X Mavericks. If you had any success doing that, or if you have a Mac and want to help, please post to the mailing list at

In other news, I think I mentioned that you can register for a new edition of my Introduction to QuantLib Development course, and with an early-bird discount to boot. If you haven't been on this blog recently, it's the course that I teach once or twice a year based on the contents of my book. A good part of the learning experience is doing a number of exercises (which you don't get to do on the blog) so you'll have to bring your laptop and get ready to code. It's going to be in London on September 22nd to 24th, and I look forward to see you there.

And finally: in next week's post, I'll need to display some math. I think I'll go with MathJax, which is supposed to work in all browser; but if yours doesn't display correctly bits of math like \( (\mathbf{x},t) \) or \( \Delta t \) here, let me know.

Follow me on Twitter if you want to be notified of new posts, or add me to your circles, or subscribe via RSS: the widgets for that are in the sidebar, at the top right of the page. Also, make sure to check my Training page.

The Monte Carlo framework

The concept behind Monte Carlo methods is deceptively simple—generate random paths, price the desired instrument over such paths, and average the results. However (and even assuming that the tasks I just enumerated were as simple as they seem: they're not) several choices can be made regarding the implementation of the model. Pseudo-random or quasi-random numbers? Which ones? When should the simulation stop? What results should be returned? How about more advanced techniques such as likelihood ratio or pathwise Greeks?

The mission of a library writer, should he accept it, is to make most (possibly all) such choices independently available to the user; to make it possible to model the underlying stochastic process in a polymorphic way; and on top of it all, to provide at least part of the scaffolding needed to put everything together and run the simulation.

The posts in this series will describe the design of our Monte Carlo framework, in what measure it meets (or falls short of) the above requirements, and how it can be integrated with the pricing-engine framework described in chapter 2.

Random-number generation

Among the basic building blocks for our task are random-number generators (RNG), the most basic being uniform RNGs. The library provides quite a few of them, the most notable being the Mersenne Twister among pseudo-random generators and Sobol among low-discrepancy ones. Here, I'll gloss over their implementation, as well as over when and how you should use either kind of generator; such information (as well as everything else related to Monte Carlo methods) is covered in literature much better than I ever could, for instance in Glasserman [1] or Jäckel [2].

The interface of uniform RNGs, shown in listing 6.1, doesn't make for a very exciting reading. Its methods are the ones you would expect: pseudo-random generators sport a constructor that takes an initialization seed and a method, next, that returns the next random value; low-discrepancy generators take an additional constructor argument that specify their dimensionality and return a vector of values from their nextSequence
method. (On second thought, we might have called the nextSequence method simply next. Including the "sequence" token in the method name reminds me of Hungarian notation, which is up there with goto in the list of things considered harmful.) The only feature of note is the Sample struct, used to return the generated values, which also contains a weight. This was done to make it possible for generators to sample values with some kind of bias (e.g., when using importance sampling). However, at this time the library doesn't contain any such generator; all available RNGs return samples with unit weight.

Listing 6.1: Basic random-number generator classes.
    template <class T>
    struct Sample {
        typedef T value_type;
        T value;
        Real weight;

    class MersenneTwisterUniformRng {
        typedef Sample<Real> sample_type;
        explicit MersenneTwisterUniformRng(unsigned long seed = 0);
        sample_type next() const;

    class SobolRsg {
        typedef Sample<std::vector<Real> > sample_type;
        SobolRsg(Size dimensionality,
                 unsigned long seed = 0);
        const sample_type& nextSequence() const;
        Size dimension() const { return dimensionality_; }

So much for uniform \acronym{RNG}s. However, they're only the rawest material; we need to refine them further. This is done by means of classes (wrappers, or combinators, if you will; or again, you might see them as implementations of the Decorator pattern) that take an instance of a uniform RNG and yield a different generator. We need two kinds of decorators. On the one hand, we need tuples of random numbers rather than single ones: namely, M × N numbers in order to evolve M variables for N steps. Low-discrepancy generators return tuples already (of course they need to do so, for their low-discrepancy property to hold at the required dimensionality). Since we want a generic interface, we'll make the pseudo-random generators return tuples as well, by simply drawing from them repeatedly; the class template that does so, called RandomSequenceGenerator, is shown in listing 6.2.

Listing 6.2: Sketch of the RandomSequenceGenerator class template.
    template<class RNG>
    class RandomSequenceGenerator {
        typedef Sample<std::vector<Real> > sample_type;
        RandomSequenceGenerator(Size dimensionality,
                                const RNG& rng)
        : dimensionality_(dimensionality), rng_(rng),
          sequence_(std::vector<Real>(dimensionality), 1.0) {}
        const sample_type& nextSequence() const {
            sequence_.weight = 1.0;
            for (Size i=0; i<dimensionality_; i++) {
                typename RNG::sample_type x(;
                sequence_.value[i] = x.value;
                sequence_.weight  *= x.weight;
            return sequence_;
        Size dimensionality_;
        RNG rng_;
        mutable sample_type sequence_;

On the other hand, we usually want Gaussian random numbers, not uniform ones. There's a number of ways to obtain them, but not all of them might be equally suited for any given situation; for instance, those consuming M uniform random numbers to yield N Gaussian ones don't work with low-discrepancy generators (especially if M varies from draw to draw, like in rejection techniques). If you want to stay generic, the simplest method is probably to feed your uniform numbers to the inverse-cumulative Gaussian function, which yield a Gaussian per uniform; this is also the default choice in the library. The corresponding class template, InverseCumulativeRsg, is shown in listing 6.3. It takes the implementation of the inverse-cumulative function as a template argument, since there's a few approximations and no closed formula for it; most of the times you'll be fine with the one provided by the library, but at times one might choose to change it in order to trade accuracy for speed (or the other way around).

Listing 6.3: Sketch of the InverseCumulativeRsg class template.
    template <class USG, class IC>
    class InverseCumulativeRsg {
        typedef Sample<std::vector<Real> > sample_type;
        InverseCumulativeRsg(const USG& uniformSequenceGenerator,
                             const IC& inverseCumulative);
        const sample_type& nextSequence() const  {
            typename USG::sample_type sample =
            x_.weight = sample.weight;
            for (Size i = 0; i < dimension_; i++) {
                x_.value[i] = ICD_(sample.value[i]);
            return x_;
        USG uniformSequenceGenerator_;
        Size dimension_;
        mutable sample_type x_;
        IC ICD_;

At this point, you might have started to have bad feelings about this framework. I have only covered random numbers, and yet we already have two of three choices to make in order to instantiate a generator. By the time we're done, we might have a dozen of template arguments on our hands. Well, bear with me for a while. In a later post, I'll describe how the library tries to mitigate the problem.

A final note: as you might have noticed, the sequence generators shown in the listings keep a mutable sequence as a data member and fill it with the new draws in order to avoid copies. Of course, this prevents instances of such classes to be shared by different threads—although I'd be a very happy camper if all the problems QuantLib had with multithreading were like this one.

Aside: the road more traveled.

The 2011 C++ standard includes RNGs as part of its standard library, based on the existing Boost implementation. At some point, we'll probably switch to their interface; but I still don't know when nor how we'll do it. We might ditch our implementation and switch to the standard classes; or we could replace the interface of our classes while keeping their current implementation; or, again, we could provide adaptors between the two interfaces. In any case, the changes will break backward compatibility, which means that they'll have to wait for an entirely new revision of QuantLib.


[1] P. Glasserman, Monte Carlo Methods in Financial Engineering. Springer, 2003.
[2] P. Jäckel, Monte Carlo Methods in Finance. John Wiley and Sons, 2002.

Liked this post? Share it: