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 quantlib-users@lists.sourceforge.net.
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 buttons for that are in the footer. 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.
So much for uniform RNGs. 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.
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.
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.
Bibliography
[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.