QuantLib: Random Number Generators

QuantLib uses random numbers in many places, such as in Pricing Engines to calculate the NPV of an instrument, in the Brownian Generator that is used in market-model simulations and in the Monte Carlo Model for path samples. In financial calculations the quality of the random number generator can be important and QuantLib provides a number of different generators that generate pseudo-random numbers and quasi-random numbers with uniform distribution and mechanisms to transform these to arbitrary probability distributions. The library distinguishes between Random Number Generators and Random Sequence Generators. Random number generators provide individual numbers while random sequence generators fill a std::vector with numbers. While random number generators can easily be used to create a random sequence, some algorithms, notably the low discrepancy sequence generators, can only be used to generate sequences and not individual random numbers.

Uniform Random Number Generators

QuantLib defines four pseudo-random number generators with uniform distribution. These are the Mersenne Twister 19937, the  L’Ecuyer generator with added Bays-Durham shuffle, Knuth’s random number generator and Luescher’s  Luxury random number generator. The last generator in this list is not implemented inside QuantLib but is a wrapper for the random number generator implemented in the Boost Random library. There are two versions of the luxury random number generator, Ranlux 3 and Ranlux 4. The second generator is slower than the first but provides more “luxury” (see M. Luescher, Computer Physics Communications, 79 (1994) pp 100-110). The number generators are defined in the following classes.

class MersenneTwisterUniformRng;
class LecuyerUniformRng;
class KnuthUniformRng;
class Ranlux3UniformRng;
class Ranlux4UniformRng;

All of these generators have constructors that take an optional initial seed as argument.

explicit MersenneTwisterUniformRng(unsigned long seed = 0);
explicit LecuyerUniformRng(long seed = 0);
explicit KnuthUniformRng(long seed = 0);
explicit Ranlux3UniformRng(Size seed = 19780503u);
explicit Ranlux4UniformRng(Size seed = 19780503u);

As you can see, different generators use slightly different integral types as seeds. The types depend on the internals of the generators. The default seed for the Mersenne Twister, L’Ecuyer and Knuth generators will cause a seed to be generated based on the system clock. The seed is calculated from the return value of the time() function using the Mersenne Twister. All this is carried out inside the SeedGenerator class.

While the details of the seeding might not be important to you there is one thing to take away from this. The seed is calculated in a deterministic way from the result of the time() function but this function has a granularity of one second. If you run the same code multiple times during one second, the random number generator is likely to produce the same values.

MersenneTwisterUniformRng provides another constructor that takes a vector of seeds which will initialise the internal states of the Mersenne Twister.

explicit MersenneTwisterUniformRng(
    const std::vector& seeds);

Each of the four random number generators has a method to get the next random number. The random number returned is a Real number between, but not including, 0 and 1.

Sample<Real> next() const;

The return type is actually not just a simple Real but a Real value together with a weight wrapped inside the Sample struct.

template <class T>
struct Sample {
    typedef T value_type;
    Sample(const T& value, Real weight)
      : value(value), weight(weight) {}
    T value;
    Real weight;

This struct is used throughout the QuantLib random number module to represent random numbers. For all samples returned by the random number generators, the weight property is set to 1.0.  To illustrate how to use the random number generators to create and print a sequence of random numbers, let’s look at a simple example.

MersenneTwisterUniformRng generator(12345);
for (int i=0; i<100; ++i) {
  std::cout << generator.next().value << std::endl;

The code above will create a Mersenne Twister and seed it with the number 12345. It will then create and print out 100 numbers.

Generic Random Sequence Generator

You can easily create a random sequence from a random number generator with the generic RandomSequenceGenerator. This class is templated with a random number generator class so it can use any generator to create the sequence.

template<class RNG>
class RandomSequenceGenerator;

The class defines a type called sample_type that defines the return type of the sequence generator. In general all sequence generators define the sample_type in the following way.

typedef Sample<std::vector<Real> > sample_type;

The generic random sequence generator has two constructors. The first argument of both constructors is the dimension, i.e. the length of the sequence that the generator will return. The second parameter is either a random number generator of the type specified by RNG or a seed number that is used to initialise the random number generator. The seed is optional and defaults to zero if not specified.

  (Size dimensionality, const RNG& rng);
  (Size dimensionality, BigNatural seed = 0);

Two methods allow access to the sequence.

const sample_type& nextSequence() const;
const sample_type& lastSequence() const;

The method nextSequence generates a new sequence with a length specified by dimensionality and returns a reference to it. Repeated access to the last sequence, without generating it again, is granted by lastSequence. The following example shows how the sequence generator might be used.

  generator(100, 12345);
std::vector<Real> values = generator.nextSequence().value;
for (int i=0; i<values.size(); ++i) { 
 std::cout << values[i] << std::endl; 

Gaussian and Other Distributions

All the previous random number and random sequence generators produce uniformly distributed random numbers. In many cases random numbers with different probability distributions are needed. The most common distribution is the Gaussian distributions. QuantLib provides two classes that generate Gaussian distributed random numbers from uniformly distributed random numbers. The BoxMullerGaussianRng uses the Box Muller transform that generates pairs of normally distributed values from pairs of uniformly distributed values. Because the Box Muller transform needs to evaluate the square root and the logarithm for every number produced it can be slow for some applications. An alternative is the application of the central limit theorem, implemented in CLGaussianRng. This generator uses the fact that the sum of 12 uniformly random numbers is already a good approximation for a normal distributed random variable. Both classes are templated with the type of the underlying random number generator.

template <class RNG>
class BoxMullerGaussianRng;

template <class RNG>
class CLGaussianRng;

Both classes are simply constructed by passing the random number generator to the constructor.

explicit BoxMullerGaussianRng::BoxMullerGaussianRng
  (const RNG& uniformGenerator);
explicit CLGaussianRng::CLGaussianRng
  (const RNG& uniformGenerator);

In addition, both classes follow the same pattern as the uniform random number generators for accessing the numbers and define the next method.

Sample<Real> next() const;

In this way, both Gaussian generators can be plugged into the RandomSequenceGenerator in the same way as any of the uniform generators. The following example will illustrate this.

typedef BoxMullerGaussianRng<MersenneTwisterUniformRng> BoxMuller;
RandomSequenceGenerator<BoxMuller> generator(BoxMuller(100, 12345)); 
std::vector<Real> values = generator.nextSequence().value; 
for (int i=0; i<values.size(); ++i) { 
  std::cout << values[i] << std::endl; 

The code above creates a vector filled with 100 normally distributed values. Inside the loop the values are then printed out.

For distributions other than Gaussian distributions QuantLib defines a general class called InverseCumulativeRng. The class creates random numbers with an arbitrary distribution from its inverse cumulative distribution. Assume that the desired probability distribution of the random numbers is given by the function \(f(x)\). Then the cumulative distribution function is then given by

$$F(x) = \int_{-\infty}^x f(x’) dx’$$

All that is needed is a function object that evaluates the inverse of the cumulative distribution \(F^{-1}(x)\). The declaration of the InverseCumulativeRng class looks like this.

template <class RNG, class IC>
class InverseCumulativeRng {
    typedef Sample<Real> sample_type;
    typedef RNG urng_type;
    explicit InverseCumulativeRng(const RNG& uniformGenerator);
    sample_type next() const;

The class IC must have a default constructor and must overload the function operator to return the inverse cumulative distribution.

Real IC::operator()(Real x);

The following example may be helpful in illustrating how to use the InverseCumulativeRng class. Let’s assume that you want to create a random variable with a probability distribution given by the logistic distribution with mean zero and scale parameter 1,

$$f(x) = \frac{1}{\left(e^{x/2} + e^{-x/2}\right)}$$

Then the cumulative distribution function is given by

$$F(x) = \int_{-\infty}^x \frac{1}{\left(e^{x’/2} + e^{-x’/2}\right)} dx’ = \frac{1}{1+e^{-x}}$$

The inverse of the function can be found by solving \(y=F(x)\) for \(x\) which results in

$$F^{-1}(y) = \mathrm{ln}\left(\frac{y}{1-y}\right)$$

All we need to do is to create a class that evaluates this function.

class ICLogistics {
    Real operator()(Real y) {
      return log(y/(1-y));

Now we can define the following random number generator.

typedef InverseCumulativeRNG<
  ICLogistics > LogisticsRNG;

In the example above we can now simply replace BoxMuller with LogisticsRNG to produce a sequence of random numbers with the logistics distribution.

Follow the author