class Parameters

{

public:

Parameters(const ParametersInner& innerObject);

Parameters(const Parameters& original);

Parameters& operator=(const Parameters& original);

virtual ˜Parameters();

60 Bridging with a virtual constructor

inline double Integral(double time1, double time2) const;

inline double IntegralSquare(double time1,

double time2) const;

double Mean(double time1, double time2) const;

double RootMeanSquare(double time1, double time2) const;

private:

ParametersInner* InnerObjectPtr;

};

inline double Parameters::Integral(double time1,

double time2) const

{

return InnerObjectPtr->Integral(time1,time2);

}

inline double Parameters::IntegralSquare(double time1,

double time2) const

{

return InnerObjectPtr->IntegralSquare(time1,time2);

}

class ParametersConstant : public ParametersInner

{

public:

ParametersConstant(double constant);

virtual ParametersInner* clone() const;

virtual double Integral(double time1, double time2) const;

virtual double IntegralSquare(double time1,

double time2) const;

private:

double Constant;

double ConstantSquare;

};

#endif

and the source ¬le

4.7 A parameters class 61

Listing 4.18 (Parameters.cpp)

#include <Parameters.h>

Parameters::Parameters(const ParametersInner& innerObject)

{

InnerObjectPtr = innerObject.clone();

}

Parameters::Parameters(const Parameters& original)

{

InnerObjectPtr = original.InnerObjectPtr->clone();

}

Parameters& Parameters::operator=(const Parameters& original)

{

if (this != &original)

{

delete InnerObjectPtr;

InnerObjectPtr = original.InnerObjectPtr->clone();

}

return *this;

}

Parameters::˜Parameters()

{

delete InnerObjectPtr;

}

double Parameters::Mean(double time1, double time2) const

{

double total = Integral(time1,time2);

return total/(time2-time1);

}

double Parameters::RootMeanSquare(double time1,

double time2) const

{

double total = IntegralSquare(time1,time2);

return total/(time2-time1);

62 Bridging with a virtual constructor

}

ParametersConstant::ParametersConstant(double constant)

{

Constant = constant;

ConstantSquare = Constant*Constant;

}

ParametersInner* ParametersConstant::clone() const

{

return new ParametersConstant(*this);

}

double ParametersConstant::Integral(double time1,

double time2) const

{

return (time2-time1)*Constant;

}

double ParametersConstant::IntegralSquare(double time1,

double time2) const

{

return (time2-time1)*ConstantSquare;

}

We have added a couple of useful, if unnecessary, methods to the Parameters

class. The Mean method returns the average value of the parameter between two

times, and is implemented by calling the Integral method which does all the

work. The RMS method returns the root-mean-square of the parameter between

two times. As it is the root-mean-square of the volatility that is appropriate to

use in the Black“Scholes formula when volatility is variable, this comes in

useful.

The rest of the code works in the same way as for the PayOff class. We give one

inherited class, the simplest possible one, ParametersConstant, which enacts a

constant parameter. As well as storing the value of the constant, we also store its

square in order to minimize time spent computing the square integral. In this case,

the saving is rather trivial, but the principle that time can be saved by computing

values, which may be needed repeatedly, once and for all in the constructor, is

worth remembering.

4.7 A parameters class 63

In SimpleMC6, we give a modi¬ed implementation of the simple Monte Carlo

which uses the new classes. First though we need

Listing 4.19 (SimpleMC6.h)

#ifndef SIMPLEMC6_H

#define SIMPLEMC6_H

#include <Vanilla3.h>

#include <Parameters.h>

double SimpleMonteCarlo4(const VanillaOption& TheOption,

double Spot,

const Parameters& Vol,

const Parameters& r,

unsigned long NumberOfPaths);

#endif

Listing 4.20 (SimpleMC6.cpp)

#include<SimpleMC6.h>

#include <Random1.h>

#include <cmath>

// the basic math functions should be in namespace

// std but aren™t in VCPP6

#if !defined(_MSC_VER)

using namespace std;

#endif

double SimpleMonteCarlo4(const VanillaOption& TheOption,

double Spot,

const Parameters& Vol,

const Parameters& r,

unsigned long NumberOfPaths)

{

double Expiry = TheOption.GetExpiry();

double variance = Vol.IntegralSquare(0,Expiry);

double rootVariance = sqrt(variance);

64 Bridging with a virtual constructor

double itoCorrection = -0.5*variance;

double movedSpot = Spot*exp(r.Integral(0,Expiry) +

itoCorrection);

double thisSpot;

double runningSum=0;

for (unsigned long i=0; i < NumberOfPaths; i++)

{

double thisGaussian = GetOneGaussianByBoxMuller();

thisSpot = movedSpot*exp( rootVariance*thisGaussian);

double thisPayOff = TheOption.OptionPayOff(thisSpot);

runningSum += thisPayOff;

}

double mean = runningSum / NumberOfPaths;

mean *= exp(-r.Integral(0,Expiry));

return mean;

}

with the obvious corresponding changes in the header ¬le. The main difference

from the previous version is that instead of computing Vol*Vol*Expiry and

r*Expiry, the IntegralSquare and Integral methods of the Parameters

class are called. Note that by forcing us to go via the methods of the

Parameters class, the new design makes the code more comprehensible. We know

immediately from looking at it that the integral of r is used, as is the square-integral

of the volatility.

In VanillaMain4.cpp, we adapt our main program to use the Parameters

class. The main difference from the previous version is that we create a

ParametersConstant object from the inputs which is then passed into the Monte

Carlo routine. Note that we do not have to create the Parameters object explicitly;

the conversion is done implicitly by the compiler. The relevant portion of code

is

4.9 Exercises 65

ParametersConstant VolParam(Vol);

ParametersConstant rParam(r);

double result = SimpleMonteCarlo4(theOption,

Spot,

VolParam,

rParam,

NumberOfPaths);

4.8 Key points

• Cloning gives us a method of implementing a virtual copy constructor.

• The rule of three says that if we need any one of copy constructor, destructor and

assignment operator then we need all three.

• We can use a wrapper class to hide all the memory handling, allowing us to treat

a polymorphic object just like any other object.

• The bridge pattern allows us to separate interface and implementation, enabling

us to vary the two independently.

• The new command is slow.

• We have to be careful to ensure that self-assignment does not cause crashes.

4.9 Exercises

Exercise 4.1 Test how fast new is on your computer and compiler. Do this by using

it to allocate an array of doubles, ten thousand times. See how the speed varies with

array size. If you have more than one compiler see how they compare. Note that

you can time things using the clock() function.

Exercise 4.2 Find out about an auto ptr. Observe that it cannot be copied in the

usual sense of copying. Show that a class with an auto ptr data member will need

a copy constructor but not a destructor.

Exercise 4.3 Implement a piecewise-constant parameters class.

5

Strategies, decoration, and statistics

5.1 Differing outputs

Our Monte Carlo routine is now easily extendible to handle any pay-off and time-

dependent parameters. However, there are plenty of valid criticisms that could still

be made. One thing that is de¬nitely lacking is the absence of any indication of

the simulation™s convergence. We could make the routine return standard error, or

a convergence table, or simply have it return the value for every single path and

analyze the results elsewhere.

As we are trying to develop an object-oriented routine, we make the statistics

gatherer an input. Thus the Monte Carlo routine will take in a statistics gatherer

object, store the results in it and the statistics gatherer will then output the statistics

as required. This technique of using an auxiliary class to decide how part of an

algorithm is implemented is sometimes called the strategy pattern.

5.2 Designing a statistics gatherer

We want our statistics gatherer to be reusable; there are plenty of circumstances

where such a routine might be useful. For example, we might have many other

Monte Carlo routines such as an exotics pricer or a BGM pricer for interest-rate

derivatives. Also, if we are developing a risk system, we might be more interested

in the ninety-¬fth percentile, or in the conditional expected shortfall, than in the

mean or variance.

What should the routine do? It must have two principal methods. The ¬rst should

take in data for each path. The second must output the desired statistics.

Since we do not wish to specify what sort of statistics are being gathered in

advance, we proceed via an abstract base class using virtual methods, just as we

did for the PayOff and Parameters classes. However, as most of the time we will

not need to copy these objects we do not bother with the bridge but work with the

base class by reference.

66

5.2 Designing a statistics gatherer 67

We have to decide the precise interface for our two principal methods. They will

be pure virtual functions declared in the base class and de¬ned in the concrete in-

herited class. Our ¬rst method, DumpOneResult, we make quite simple: it takes

in a double and returns nothing. Note that it is not a const method, since by its

very nature it must update the statistics stored inside the object. Note that we have

not allowed the possibility of dumping more than one value per path, which could

be argued to be a defect. The object will store what it needs to in order to com-

pute the statistics desired and no more. So if statistics gatherer™s job is to compute

the mean, then it need only store the running sum and the number of paths. How-

ever, for a more complicated statistic we might need the value for every path to be

stored.

Our second method, which will indeed return the results, requires a little more

thought. We have to decide what sort of object to return the results in. Another

issue is whether it should be possible to ask for statistics en route or ought we be

able to call the method for returning results only once.

With regard to the form in which to return results, we opt for a vector of

vectors. This will allow us to easily return a table if we so desire. Whilst this

would not be a great way to implement a matrix-class if we were doing linear al-

gebra, the return statistics are not a matrix, just a table and this is suf¬cient for our

purposes.

We opt to allow the return statistics method to be called many times and therefore

name it GetResultsSoFar. This will cost us little (possibly nothing), and will be

more robust than an object that crashes if the get results method is called twice.

We make it a const method as it should not change the state of the object in any

substantive way: this enforces the rule that the method can be called many times.

Listing 5.1 (MCStatistics.h)

#ifndef STATISTICS_H

#define STATISTICS_H

#include <vector>

class StatisticsMC

{

public:

StatisticsMC(){}

virtual void DumpOneResult(double result)=0;

virtual std::vector<std::vector<double> >

GetResultsSoFar() const=0;

68 Strategies, decoration, and statistics

virtual StatisticsMC* clone() const=0;

virtual ˜StatisticsMC(){}

private:

};

class StatisticsMean : public StatisticsMC

{

public:

StatisticsMean();

virtual void DumpOneResult(double result);

virtual std::vector<std::vector<double> >

GetResultsSoFar() const;

virtual StatisticsMC* clone() const;

private:

double RunningSum;

unsigned long PathsDone;

};

#endif

Our abstract base class is StatisticsMC. It has the pure virtual functions Dump-

OneResult and GetResultsSoFar. We include the clone method to allow for the

possibility of virtual copy construction. We also make the destructor virtual, as any

cloned objects will likely need to be destroyed via pointers to the base class which

will not know their type, as usual. The base class does nothing but the important

task of de¬ning an interface.

We give a very simple inherited class StatisticsMean, which returns the mean

of the simulation, just as our routine previously did. The source code is included in

MCStatistics.cpp.

Listing 5.2 (MCStatistics.cpp)

#include<MCStatistics.h>

using namespace std;

StatisticsMean::StatisticsMean()

:

RunningSum(0.0), PathsDone(0UL)

{

5.3 Using the statistics gatherer 69

}

void StatisticsMean::DumpOneResult(double result)

{

PathsDone++;

RunningSum += result;

}

vector<vector<double> >

StatisticsMean::GetResultsSoFar() const

{

vector<vector<double> > Results(1);

Results[0].resize(1);

Results[0][0] = RunningSum / PathsDone;

return Results;

}

StatisticsMC* StatisticsMean::clone() const

{

return new StatisticsMean(*this);

}

Note that whilst we write the DumpOneResult method to be ef¬cient since it

will be called in every iteration of the loop, we do not worry about ef¬ciency for

GetResultsSoFar, as it will generally be called only once per simulation.

5.3 Using the statistics gatherer

Having written a class for gathering statistics, we now need to modify our Monte

Carlo routine to make use of it. We do this in the function SimpleMonteCarlo5

which is declared in SimpleMC7.h.

Listing 5.3 (SimpleMC7.h)

#ifndef SIMPLEMC7_H

#define SIMPLEMC7_H

#include <Vanilla3.h>

#include <Parameters.h>

70 Strategies, decoration, and statistics

#include <MCStatistics.h>

void SimpleMonteCarlo5(const VanillaOption& TheOption,

double Spot,

const Parameters& Vol,

const Parameters& r,

unsigned long NumberOfPaths,

StatisticsMC& gatherer);

#endif

Note that the StatisticsMC object is passed in by reference and is not const.

This is crucial as we want the object passed in to gather the information given to

it and for this information to be available after the function has returned. If we had

passed by value then the object outside would not change, and all the results would

disappear at the end of the function which we emphatically do not want. If the

object was const, then it would not be possible to put any new data into it which

would be useless. Previously, our routine returned a double; now it is void as all

the data to be returned is inside the object gatherer.

We de¬ne the function in SimpleMC7.cpp:

Listing 5.4 (SimpleMC7.cpp)

#include <SimpleMC7.h>

#include <Random1.h>

#include <cmath>

// the basic math functions should be

// in namespace std but aren™t in VCPP6

#if !defined(_MSC_VER)

using namespace std;

#endif

void SimpleMonteCarlo5(const VanillaOption& TheOption,

double Spot,

const Parameters& Vol,

const Parameters& r,

unsigned long NumberOfPaths,

StatisticsMC& gatherer)

{

double Expiry = TheOption.GetExpiry();

double variance = Vol.IntegralSquare(0,Expiry);

double rootVariance = sqrt(variance);

5.3 Using the statistics gatherer 71

double itoCorrection = -0.5*variance;

double movedSpot =

Spot*exp(r.Integral(0,Expiry)+itoCorrection);

double thisSpot;

double discounting = exp(-r.Integral(0,Expiry));

for (unsigned long i=0; i < NumberOfPaths; i++)

{

double thisGaussian = GetOneGaussianByBoxMuller();

thisSpot = movedSpot*exp( rootVariance*thisGaussian);

double thisPayOff = TheOption.OptionPayOff(thisSpot);

gatherer.DumpOneResult(thisPayOff*discounting);

}

return;

}

Our routine appears simpler than SimpleMonteCarlo4; all the work previously

spent accounting the results is now sublimated into the line on which we call

.DumpOneResult. Of course, the code has not disappeared; it has simply moved

into a different ¬le. Thus the strategy pattern gives us a readability bene¬t as well

as ¬‚exibility.

We illustrate how the gatherer might be called in StatsMain1.cpp:

Listing 5.5 (StatsMain1.cpp)

/*

uses source files

MCStatistics.cpp,

Parameters.cpp,

PayOff3.cpp,

PayOffBridge.cpp,

Random1.cpp,

SimpleMC7.cpp,

Vanilla3.cpp,

*/

#include<SimpleMC7.h>

#include<iostream>

using namespace std;

#include<Vanilla3.h>

#include<MCStatistics.h>

72 Strategies, decoration, and statistics

int main()

{

double Expiry;

double Strike;

double Spot;

double Vol;

double r;

unsigned long NumberOfPaths;

cout << "\nEnter expiry\n";

cin >> Expiry;

cout << "\nStrike\n";

cin >> Strike;

cout << "\nEnter spot\n";

cin >> Spot;

cout << "\nEnter vol\n";

cin >> Vol;

cout << "\nr\n";

cin >> r;

cout << "\nNumber of paths\n";

cin >> NumberOfPaths;

PayOffCall thePayOff(Strike);

VanillaOption theOption(thePayOff, Expiry);

ParametersConstant VolParam(Vol);

ParametersConstant rParam(r);

StatisticsMean gatherer;

SimpleMonteCarlo5(theOption,

Spot,

VolParam,

rParam,

5.4 Templates and wrappers 73

NumberOfPaths,

gatherer);

vector<vector<double> > results = gatherer.GetResultsSoFar();

cout <<"\nFor the call price the results are \n";

for (unsigned long i=0; i < results.size(); i++)

{

for (unsigned long j=0; j < results[i].size(); j++)

cout << results[i][j] << " ";

cout << "\n";

}

double tmp;

cin >> tmp;

return 0;

}

Our output of the results is now a bit more complicated in that we have to loop

over the vector of vectors in order to display the results. Of course, for the

particular class we have de¬ned, only one number is returned so this is not strictly

necessary. However, by writing it for the most general sort of return statement

from a statistics gatherer, we produce more robust code. After all, the object has

contracted to return a vector of vectors, and we should code in accordance with

this contract, and make no extra assumptions.

5.4 Templates and wrappers

We have created a class hierarchy for gathering statistics. This hierarchy includes

a virtual constructor, clone, so we can copy these objects without knowing their

type. However, if we wish to start copying and storing the objects then we have a

slight issue in that, just as with our PayOff class, this will have to be done manually

unless we write an extra wrapper class to do it for us. In the next section, we give an

example where this copying will be necessary, so we do need to provide a wrapper

class.

It becomes clear at this point that we will need to write these wrapper classes

repeatedly in a very similar way. We therefore present a templatized solution. The

name of the base class to be wrapped is taken as an argument to the wrapper class

and it can then be used for any class which provides a clone method.

74 Strategies, decoration, and statistics

Our Wrapper class provides various functionalities which are intended to make

it act like a pointer to a single object but with added responsibilities. The added

responsibilities are that the pointer is both responsible for and owns the object

pointed to at all times. Thus if we copy the Wrapper object, the pointed-to object

is also copied, so that each Wrapper object has its own copy of the pointed-to

object. When the Wrapper object ceases to exist because of going out of scope, or

being deleted, the pointed-to object is automatically deleted as well.

If we set one Wrapper object equal to another, then the object previously pointed

to must be deleted, and then a copy of the new object must be created so each

Wrapper still owns precisely one object.

In addition, it must be possible to dereference the Wrapper to obtain the under-

lying object. In other words, if you put *mywrapper then you should obtain the

object pointed to by mywrapper. We do this by overloading the operator*() and

just make it return the dereferenced inner pointer.

We also want to be able to access the methods of the inner object. Whilst this can

always be done by putting (*mywrapper).theMethod(), it is a lot less elegant

than being able to type myWrapper->theMethod(), which is the normal code for

an ordinary pointer. We provide this functionality by overloading operator->().

We provide the relevant code in Wrapper.h.

Listing 5.6 (Wrapper.h)

#ifndef WRAPPER_H

#define WRAPPER_H

template< class T>

class Wrapper

{

public:

Wrapper()

{ DataPtr =0;}

Wrapper(const T& inner)

{

DataPtr = inner.clone();

}

˜Wrapper()

{

if (DataPtr !=0)

5.4 Templates and wrappers 75

delete DataPtr;

}

Wrapper(const Wrapper<T>& original)

{

if (original.DataPtr !=0)

DataPtr = original.DataPtr->clone();

else

DataPtr=0;

}

Wrapper& operator=(const Wrapper<T>& original)

{

if (this != &original)

{

if (DataPtr!=0)

delete DataPtr;

DataPtr = (original.DataPtr !=0) ?

original.DataPtr->clone() : 0;

}

return *this;

}

T& operator*()

{

return *DataPtr;

}

const T& operator*() const

{

return *DataPtr;

}

const T* const operator->() const

{

return DataPtr;

}

76 Strategies, decoration, and statistics

T* operator->()

{

return DataPtr;

}

private:

T* DataPtr;

};

#endif

We start before each declaration with the command template<class T> to

let the compiler know we are writing template code. The class T will be spec-

i¬ed elsewhere. The compiler will produce one copy of the code for each dif-

ferent sort of T that is used. Thus if we declare Wrapper<MCStatistics>

TheStatsGatherer;, the compiler will then proceed to create the code by sub-

stituting MCStatistics for T everywhere, and then compile it. This has some side

effects: the ¬rst is that all the code for the Wrapper template is in the header ¬le “

there is no wrapper.cpp ¬le. The second is that if we use the Wrapper class many

times, we have to compile a lot more code than we might actually expect. Whilst

this is not really an issue for this class, it could be one for a complicated class,

where we might end up with rather slow compile times and a much larger than ex-

pected executable. There are some other effects and we will return to this topic in

Section 9.6.

We provide the Wrapper class with a default constructor. This means that it is

possible to have a Wrapper object which points to nothing. If we did not then a

declaration such as

std::vector<Wrapper<MCStatistics> > StatisticsGatherers(10);

would not compile: the constructor for vector would look for the default con-

structor for Wrapper<MCStatistics> in order to create the ten copies speci¬ed,

and not ¬nd it. Why would we want a vector of Wrappers? We saw in Section

3.4 that we can get into trouble if we try to copy inherited objects into base class

objects without a wrapper class. The same reasons apply here. We cannot declare

a vector of base class objects as they are abstract, and even if they were not, they

would be the wrong size. We therefore have to store pointers, references or wrap-

pers, and wrappers are the easiest option; they take care of all the memory handling

for us.

Given that a Wrapper object can point to nothing, we have to be able to take this

into account when writing the class™s methods. We indicate that the object points to

5.5 A convergence table 77

nothing by setting the pointer to zero. When carrying out copying and assignment,

we then have to take care of this special case.

We provide two different versions of the dereferencing operator *, as it should be

possible to dereference both const and non-const objects. As one would expect,

the const version returns a const object and the non-const version does not.

We have declared the two operators inline to ensure that there is no performance

overhead induced by going via a wrapper.

Similarly, we declare the operator-> to have both const and non-const ver-

sions. The syntax here is a little strange in that all the operator does is return the

pointer. However, there are special rules for overloading -> which ensure that any

method following -> is correctly invoked for the pointer returned.

5.5 A convergence table

If we use a statistics gatherer and run the simulation it will tell us the relevant

statistics for the entire simulation. However, it does not necessarily give us a feel

for how well the simulation has converged. One standard method of checking the

convergence is to examine the standard error of the simulation; that is measure the

sample standard deviation and divide by the square root of the number of paths. If

one is using low-discrepancy numbers this measure does not take account of their

special properties and, in fact, it predicts the same error as for a pseudo-random

simulation (see for example [10]). Which we can expect to be too large. (Else why

use low-discrepancy numbers?)

One alternative method is therefore to use a convergence table. Rather than re-

turning statistics for the entire simulation, we instead return them for every power

of two to get an idea of how the numbers are varying. We could just write a class

directly to return such a table for the mean, but since we might want to do this for

any statistic, we do it in a reusable fashion.

Our class must contain a statistics gatherer in order to decide for which statis-

tics to create a convergence table. On the other hand, it must implement the same

interface as all the other statistics gatherers so we can plug it into the same sim-

ulations. We therefore de¬ne a class ConvergenceTable which is inherited from

MCStatistics, and has a wrapper of an MCStatistics object as a data member.

The fact that the class is inherited from MCStatistics guarantees that from the

outside it looks just like any other statistics-gatherer object. The difference on the

inside is that we can make the data member refer to any kind of statistics gatherer

we like, and so we have a convergence table for any statistic for which a statistics

gatherer has been written. We give the implementation in ConvergenceTable.h

and ConvergenceTable.cpp.

78 Strategies, decoration, and statistics

Listing 5.7 (ConvergenceTable.h)

#ifndef CONVERGENCE_TABLE_H

#define CONVERGENCE_TABLE_H

#include <MCStatistics.h>

#include <wrapper.h>

class ConvergenceTable : public StatisticsMC

{

public:

ConvergenceTable(const Wrapper<StatisticsMC>& Inner_);

virtual StatisticsMC* clone() const;

virtual void DumpOneResult(double result);

virtual std::vector<std::vector<double> >

GetResultsSoFar() const;

private:

Wrapper<StatisticsMC> Inner;

std::vector<std::vector<double> > ResultsSoFar;

unsigned long StoppingPoint;

unsigned long PathsDone;

};

#endif

Listing 5.8 (ConvergenceTable.cpp)

#include<ConvergenceTable.h>

ConvergenceTable::ConvergenceTable(const

Wrapper<StatisticsMC>& Inner_)

: Inner(Inner_)

{

StoppingPoint=2;

PathsDone=0;

}

StatisticsMC* ConvergenceTable::clone() const

{

return new ConvergenceTable(*this);

}

5.5 A convergence table 79

void ConvergenceTable::DumpOneResult(double result)

{

Inner->DumpOneResult(result);

++PathsDone;

if (PathsDone == StoppingPoint)

{

StoppingPoint*=2;

std::vector<std::vector<double> >

thisResult(Inner->GetResultsSoFar());

for (unsigned long i=0; i < thisResult.size(); i++)

{

thisResult[i].push_back(PathsDone);

ResultsSoFar.push_back(thisResult[i]);

}

}

return;

}

std::vector<std::vector<double> >

ConvergenceTable::GetResultsSoFar() const

{

std::vector<std::vector<double> > tmp(ResultsSoFar);

if (PathsDone*2 != StoppingPoint)

{

std::vector<std::vector<double> >

thisResult(Inner->GetResultsSoFar());

for (unsigned long i=0; i < thisResult.size(); i++)

{

thisResult[i].push_back(PathsDone);

tmp.push_back(thisResult[i]);

}

}

return tmp;

}

Note that we do not write a copy constructor, destructor or assignment operator

as the class itself does no dynamic memory allocation. Dynamic memory allocation

80 Strategies, decoration, and statistics

does occur inside the class but it is all handled automatically by the Wrapper tem-

plate class.

The class does not do a huge amount; every result passed in is passed to the

inner class. When we reach a point where the number of paths done is a multiple of

two, the inner class™s GetResults() method is called, and the results stored with

the number of paths done so far added in. When the class™s own GetResults()

methods is called, it calls the inner class™s method one more time if necessary and

then spits out all the stored results.

In StatsMain2.cpp, we illustrate how the routine might be called:

Listing 5.9

StatisticsMean gatherer;

ConvergenceTable gathererTwo(gatherer);

SimpleMonteCarlo5(theOption,

Spot,

VolParam,

rParam,

NumberOfPaths,

gathererTwo);

vector<vector<double> > results =

gathererTwo.GetResultsSoFar();

First create a StatisticsMean object: then pass it into a ConvergenceTable

object, gatherTwo. Note the constructor takes a Wrapper<MCStatistics> object

but the compiler happily does this conversion for us. We then pass the new gatherer

into SimpleMonteCarlo5 which has not required any changes. We have also not

made any changes to either of the MCStatistics ¬les.

5.6 Decoration

The technique of the last section is an example of a standard design pattern called

the decorator pattern. We have added functionality to a class without changing the

interface. This process is called decoration. The most important point is that, since

the decorated class has the same interface as the undecorated class, any decoration

which can be applied to the original class can also be applied to the decorated

class.

We can therefore decorate as many times as we wish. It would be syntactically

legal (but not useful), for example, to have a convergence table of convergence

tables. We will more often wish to decorate several times but in differing manners.

5.8 Exercises 81

How else might we want to decorate? If we have a stream of numbers de¬ning

a time series, we often want the statistics of the successive increments instead of

the numbers themselves. A decorator class could therefore do this differencing and

pass the difference into the inner class.

We might want more than one statistic for a given set of numbers; rather than

writing one class to gather many statistics, we could write a decorator class which

contains a vector of statistics gatherers and passes the gathered value to each one

individually. The GetResults() method would then garner the results from each

of the inner gatherers and collate them.

We can also apply these decoration ideas to the Parameters class. We could

de¬ne a class that takes the linear multiple of an inner Parameters object for

example. This class would simple multiply the integral by a given constant, and

the square integral by its square.

5.7 Key points

In this chapter, we have seen that we can allow the user to specify aspects of how

an algorithm works by making part of the algorithm be carried out in an inputted

class. We have also examined the techniques of decoration and templatization.

• Routines can be made more ¬‚exible by using the strategy pattern.

• Making part of an algorithm be implemented by an inputted class is called the

strategy pattern.

• For code that is very similar across many different classes, we can use templates

to save time in rewriting.

• If we want containers of polymorphic objects, we must use wrappers or pointers.

• Decoration is the technique of adding functionality by placing a class around a

class which has the same interface; i.e. the outer class is inherited from the same

base class.

• A class can be decorated several times.

5.8 Exercises

Exercise 5.1 Write a statistics gathering class that computes the ¬rst four moments

of a sample.

Exercise 5.2 Write a statistics gathering class that computes the value at risk of a

sample.

Exercise 5.3 Write a statistics gathering class that allows the computation of sev-

eral statistics via inputted classes.

82 Strategies, decoration, and statistics

Exercise 5.4 Use the strategy pattern to allow the user to specify termination con-

ditions for the Monte Carlo, e.g., time spent or paths done.

Exercise 5.5 Write a terminator class that causes termination when either of two

inner terminator classes speci¬es termination.

Exercise 5.6 * Write a template class that implements a reference counted wrapper.

This will be similar to the wrapper class but instead of making a clone of the inner

object when the wrapper is copied, an internal counter is increased and the inner

object is shared. When a copy is destroyed, the inner counter is decremented. When

the inner counter reaches zero, the object is destroyed. Note that both the counter

and the inner object will be shared across copies of the object. (This exercise is

harder than most in this book.)

6

A random numbers class

6.1 Why?

So far, we have been using the inbuilt random number generator, rand. In this

chapter, we look at how we might implement a class to encapsulate random number

generation. There are a number of reasons we might wish to do this.

rand is implementation dependent. The standard speci¬es certain properties of

rand and gives an example of how it could be implemented but it does not actually

specify the details. This has important consequences for us. The ¬rst is simply that

we cannot expect any consistency across compilers. If we decide to test our code

by running it on multiple platforms, we can expect to obtain differing streams of

random numbers and whilst our Monte Carlo simulations should still converge to

the same number, this is a lot weaker than having every single random draw match-

ing. Thus our code becomes harder to test. A second issue is that we do not know

how good the compiler™s implementation is. Either we have to get hold of technical

documents for every compiler we use and make sure that the implementors have

done a good job, or we have to run a number of statistical tests to ensure that rand

is up to the job. Note that for most simulations we will actually need many random

draws for each path, and so it is not enough for us to check that single draws do a

good job of simulating the uniform distribution; instead we need a large number of

successive draws to do a good job of ¬lling out the unit hypercube, which is much

tougher.

rand is not predictable. A crucial aspect of running Monte Carlo simulations

is that they must be reproducible. If we run the same simulation twice we want

to obtain precisely the same random numbers. We can achieve this with rand by

using the srand command to set the seed which will guarantee the same number

stream from rand every time. The problem is that the seed is a global variable.

This means that calling rand in different parts of the program will cause totally

unrelated pieces of code to affect each other™s operation. We therefore want to be

83

84 A random numbers class

able to insulate the random number stream used by a particular simulation from the

rest of the program.

Another advantage of using a class is that we can decorate it. For example, sup-

pose we wish to use anti-thetic sampling. We could write a decorator class that

does anti-thetic sampling. This can then be combined with any random number

generator we have written, and plugged into the Monte Carlo simulator, with no

changes to the simulator class. If we used rand directly we would have to ¬ddle

with the guts of the simulator class. Similarly, if we wish to carry out moment

matching we could use a decorator class and then plug the decorated class into the

simulator.

A further reason is that we might decide not to use pseudo-random (i.e. random)

numbers but low-discrepancy numbers instead. Low-discrepancy numbers (some-

times called quasi-random numbers) are sequences of numbers designed to do a

good job of ¬lling out space. They are therefore anything but random. However,

they have the right statistical properties to guarantee that simulations converge to

the correct answer. Their space-¬lling properties mean they make simulations con-

verge faster. If we are using a random number class, we could replace this class

with a generator for low-discrepancy numbers without changing the interior of our

code.

6.2 Design considerations

As we want the possibility of having many random number generators and we want

to be able to add new ones later on without recoding, we use an abstract base class

to specify an interface. Each individual generator will then be inherited from it. In

order to specify the interface, we have to identify what we want from any random

number class.

Generally, when working with any Monte Carlo simulation, the simulation will

have a dimensionality which is the number of random draws needed to simulate

one path. This number is equal to the number of variables of integration in the

underlying integral which we are trying to approximate. It is generally cleaner

therefore to obtain all the draws necessary for a path in one go. This has the added

advantage that a random number generator can protest (i.e. throw an error) if it is

being used beyond its dimensional speci¬cation. Additionally, when working with

low-discrepancy numbers it is essential that the generator know the dimensionality

as the generator has to be set up speci¬cally for each choice of dimension.

This means that we need methods to set the dimensionality, and to obtain an

array of uniforms of size dimensionality from the generator. We also provide a

method that states the dimensionality.

6.2 Design considerations 85

For ¬nancial applications, we will want standard Gaussian draws more often

than uniforms so we will want a method of obtaining them instead. In fact, we

can separate out the creation of the uniforms and their conversion into Gaussians.

The conversion into Gaussians can therefore be done in a generator-independent

fashion and this means that it can be implemented as a method of the base class

which calls the virtual method that creates the uniform draws.

What else might we want? For many applications, it is necessary to generate

the same stream of random numbers twice. For example, if we wish to compute

Greeks by bumping parameters, the error is much smaller if the same numbers are

used twice. (See for example [11] or [13].) Or if we wish to carry out moment

matching, the reuse of the same random numbers stream twice enables us to avoid

storing all the numbers generated. Thus we include methods to reset the generator

to its initial state, and to set the seed of the generator.

Occasionally, we wish to be sure of having a different stream of random num-

bers. For example, when carrying out an optimization in order to estimate an exer-

cise strategy, we generally use one set of random numbers to optimize parameters

for the strategy, and then having chosen the strategy we run a second simulation

with different random numbers to estimate the price. This allows us to be sure

that the optimization has not exploited the micro-structure of the random number

stream. A simple way to achieve the differing streams of numbers is to make sure

the generator skips a number of paths equal to the number used for the ¬rst simu-

lation. We therefore include a method which allows us to skip paths.

Finally, we may wish to copy a random number generator for which we do not

know the type. We therefore include a clone method to enable virtual construc-

tion.

One extra issue we have to think about is in what range a uniform should lie. The

uniform distribution is generally de¬ned to be a density function on the interval

[0, 1] such that the probability that a draw X lies in an interval of length ± is ±.

The subtlety lies in whether we allow the values 0 and 1 to be taken. Since taking

either value is a probability zero event allowing or disallowing either value will

not effect the statistical properties of our simulation, but they can have practical

effects. For example, if we elect to convert the uniforms into Gaussians by using

the inverse cumulative normal function (which we will) then the numbers 0 and 1

cause us dif¬culties since the inverse cumulative normal function naturally maps

them to ’∞ and +∞. To avoid these dif¬culties, we therefore require that our

uniform variates never take these values and thus lie in the open interval (0, 1).

The main side effect of this choice is that if we use random generators written by

others then we need to check that they satisfy the same convention, and if not, adapt

them appropriately.

86 A random numbers class

6.3 The base class

We specify the interface via a base class as follows, Random2.h,

Listing 6.1 (Random2.h)

#ifndef RANDOM2_H

#define RANDOM2_H

#include <Arrays.h>

class RandomBase

{

public:

RandomBase(unsigned long Dimensionality);

inline unsigned long GetDimensionality() const;

virtual RandomBase* clone() const=0;

virtual void GetUniforms(MJArray& variates)=0;

virtual void Skip(unsigned long numberOfPaths)=0;

virtual void SetSeed(unsigned long Seed) =0;

virtual void Reset()=0;

virtual void GetGaussians(MJArray& variates);

virtual void ResetDimensionality(unsigned long

NewDimensionality);

private:

unsigned long Dimensionality;

};

unsigned long RandomBase::GetDimensionality() const

{

return Dimensionality;

}

#endif

Whilst most of the methods of RandomBase are pure virtual, three are not. The

method GetGaussians transforms uniforms obtained from the GetUniforms

method into standard Gaussian distributions. It does this via an approximation to

6.3 The base class 87

the inverse cumulative normal function due to Moro, [21]. As this method only

uses one uniform to produce a Gaussian and enacts precisely the de¬nition of the

Gaussian distribution it is very robust and works under all circumstances. Never-

theless, we make the method virtual to allow the possibility that for a particular

generator there is another preferred conversion method. Or even to allow the pos-

sibility that the generator provides normals which are then converted into uniforms

by the GetUniforms method.

The GetDimensionality method simply returns the dimensionality of the gen-

erator and there is no need for it to be virtual.

We also have the concrete virtual function ResetDimensionality. As the base

class stores dimensionality, it must be told when dimensionality changes: that is

the purpose of this function. However, the function is virtual because generally

if dimensionality changes, the random number generator will also need to know.

Suppose we have overriden this virtual function in an inherited class. Calling the

method thus only calls the inherited class method and the base class method is ig-

nored; however, we still need the base class method to be called; this has to be done

by the inherited class method. The syntax to do this is to pre¬x the method with

RandomBase::. The compiler then ignores the virtual function table and instead

knows to call the method associated to the base class.

Note that we de¬ne the interface for GetUniforms and GetGaussians via a ref-

erence to an array. The reason we do this is that we do not wish to waste time copy-

ing arrays. Also remember that arrays of dynamic size generally involve dynamic

memory allocation, i.e. new, and therefore are quite slow to create and to destroy.

We want to minimize unnecessary operations, and by passing the return values into

a pre-generated array we avoid all this. The array class used here is quite simple

and given in Appendix C. We assume that the array is of suf¬cient size. We could

check that it is big enough but that could result in substantial overhead. One solu-

tion would be to check the size only if a compiler ¬‚ag was set, e.g. in debug mode.

Note that one disadvantage of this approach is that we are now bound to this

array class. How could we overcome that disadvantage? One solution would be to

simply pass in a pointer, and write to the memory locations pointed to. However,

the use of raw pointers tends to lead to code that is hard to debug, and is therefore

best avoided. Another solution is to templatize so that the array class is a template

argument and the code will then work with any array class which has the requisite

methods. A related solution is to use iterators. An iterator is a generalization of

a pointer and we could templatize the code to work off any iterator. We do not

explore these options here but the reader should bear them in mind if he wishes to

adapt the code.

The source code for the base class is quite simple as it does not do very much:

88 A random numbers class

Listing 6.2 (Random2.cpp)

#include <Random2.h>

#include <Normals.h>

#include <cstdlib>

// the basic math functions should be in namespace

// std but aren™t in VCPP6

#if !defined(_MSC_VER)

using namespace std;

#endif

void RandomBase::GetGaussians(MJArray& variates)

{

GetUniforms(variates);

for (unsigned long i=0; i < Dimensionality; i++)

{

double x=variates[i];

variates[i] = InverseCumulativeNormal(x);

}

}

void RandomBase::ResetDimensionality(unsigned long

NewDimensionality)

{

Dimensionality = NewDimensionality;

}

RandomBase::RandomBase(unsigned long Dimensionality_)

: Dimensionality(Dimensionality_)

{

}

The inverse cumulative normal function is included in the ¬le Normals and is a

piece-wise rational approximation. See Appendix B.

6.4 A linear congruential generator and the adapter pattern

We now need to actually write a random number generator. A simple method

of generating random numbers is a linear congruential generator. We present a

6.4 A linear congruential generator and the adapter pattern 89

generator called by Park & Miller the minimal standard generator. In other words,

it is a generator that provides a minimum guaranteed level of statistical accuracy.

We refer the reader to [28] for further discussion of this and many other random

number generators.

We present the generator in two pieces. A small inner class that develops a ran-

dom generator that returns one integer (i.e., long) every time it is called, and a

larger class that turns the output into a vector of uniforms in the format desired. We

present the class de¬nition in ParkMiller.h.

Listing 6.3 (ParkMiller.h)

#ifndef PARK_MILLER_H

#define PARK_MILLER_H

#include <Random2.h>

class ParkMiller

{

public:

ParkMiller(long Seed = 1);

long GetOneRandomInteger();

void SetSeed(long Seed);

static unsigned long Max();

static unsigned long Min();

private:

long Seed;

};

class RandomParkMiller : public RandomBase

{

public:

RandomParkMiller(unsigned long Dimensionality,

unsigned long Seed=1);

virtual RandomBase* clone() const;

virtual void GetUniforms(MJArray& variates);

virtual void Skip(unsigned long numberOfPaths);

virtual void SetSeed(unsigned long Seed);

virtual void Reset();

90 A random numbers class

virtual void ResetDimensionality(unsigned long

NewDimensionality);

private:

ParkMiller InnerGenerator;

unsigned long InitialSeed;

double Reciprocal;

};

#endif

The inner class is quite simple. It develops a sequence of uncorrelated longs. The

seed can be set either in the constructor or via a set seed method. We give two extra

methods which indicate the minimum and maximum values that the generator can

give out. Such information is crucial to a user who wishes to convert the output into

uniforms, as they will need to subtract the minimum and divide by the maximum

minus the minimum to get a number in the interval [0, 1].

The bigger class is inherited from RandomBase. It has all the methods that it re-

quires. Its main data member is a ParkMiller generator object. It also remembers

the initial seed, and the reciprocal of the maximum value plus one, to save time

then turning the output of the inner generator into uniforms.

Our pattern here is an example of the adapter pattern. We have a random gen-

erator which works and is effective, however its interface is not what the rest of

the code expects. We therefore write a class around it which adapts its interface

into what we want. Whenever we use old code or import libraries, it is rare for the

interfaces to ¬t precisely with what we have been using, and the adapter pattern

is then necessary. To use the adapter pattern simply means to use an intermediary

class which transforms one interface into another. It is the coding equivalent of a

plug adapter.

The implementation of these classes is straightforward. The generator relies on

modular arithmetic. The basic idea is that if you repeatedly multiply a number by

a large number, and then take the modulus with respect to another number, then

the successive remainders are effectively random. We refer the reader to [28] for

discussion of the mathematics and the choice of the constants.

Listing 6.4 (ParkMiller.cpp)

#include <ParkMiller.h>

const long a = 16807;

const long m = 2147483647;

const long q = 127773;

const long r = 2836;

6.4 A linear congruential generator and the adapter pattern 91

ParkMiller::ParkMiller(long Seed_ ) : Seed(Seed_)

{

if (Seed ==0)

Seed=1;

}

void ParkMiller::SetSeed(long Seed_)

{

Seed=Seed_;

if (Seed ==0)

Seed=1;

}

unsigned long ParkMiller::Max()

{

return m-1;

}

unsigned long ParkMiller::Min()

{

return 1;

}

long ParkMiller::GetOneRandomInteger()

{

long k;

k=Seed/q;

Seed=a*(Seed-k*q)-r*k;

if (Seed < 0)

Seed += m;

return Seed;

}

RandomParkMiller::RandomParkMiller(unsigned long Dimensionality,

unsigned long Seed)

: RandomBase(Dimensionality),

InnerGenerator(Seed),

92 A random numbers class

InitialSeed(Seed)

{

Reciprocal = 1/(1.0+InnerGenerator.Max());

}

RandomBase* RandomParkMiller::clone() const

{

return new RandomParkMiller(*this);

}

void RandomParkMiller::GetUniforms(MJArray& variates)

{

for (unsigned long j=0; j < GetDimensionality(); j++)

variates[j] =

InnerGenerator.GetOneRandomInteger()*Reciprocal;

}

void RandomParkMiller::Skip(unsigned long numberOfPaths)

{

MJArray tmp(GetDimensionality());

for (unsigned long j=0; j < numberOfPaths; j++)

GetUniforms(tmp);

}

void RandomParkMiller::SetSeed(unsigned long Seed)

{

InitialSeed = Seed;

InnerGenerator.SetSeed(Seed);

}

void RandomParkMiller::Reset()

{

InnerGenerator.SetSeed(InitialSeed);

}

void RandomParkMiller::ResetDimensionality(unsigned long

NewDimensionality)

{

RandomBase::ResetDimensionality(NewDimensionality);

InnerGenerator.SetSeed(InitialSeed);

}

6.5 Anti-thetic sampling via decoration 93

Note that we check whether the seed is zero. If it is we change it to 1. The reason

is that a zero seed yields a chain of zeros. Note the advantage of a class-based

implementation here. The seed is only inputted in the constructor and the set seed

method, which are called only rarely, so we can put in extra tests to make sure

the seed is correct with no real overhead. If the seed had to be checked every time

the random number generator was called, then the overhead would be substantial

indeed.

The implementation of the adapter class is quite straightforward. Note that we

divide the outputs of the inner class by the maximum plus 1, and so ensure that we

obtain random numbers on the open interval (0, 1) rather than the closed one; this

means that we will have no trouble with the inverse cumulative normal function.

6.5 Anti-thetic sampling via decoration

A standard method of improving the convergence of Monte Carlo simulations is

anti-thetic sampling. The idea is very simple, if a X is a draw from a standard

Gaussian distribution so is ’X . This means that if we draw a vector (X 1 , . . . , X n )

for one path then instead of drawing a new vector for the next path we simply

use (’X 1 , . . . , ’X n ). This method guarantees that, for any even number of paths

drawn, all the odd moments of the sample of Gaussian variates drawn are zero,

and in particular the mean is correct. This generally, but not always, causes simula-

tions to converge faster. See [11] for discussion of the pros and cons of anti-thetic

sampling.

We wish to implement anti-thetic sampling in such a way that it can be used

with any random number generator and with any Monte Carlo simulation in such

a way that we only have to implement it once. The natural way to do this is the

decorator pattern. The decoration can be applied to any generator so it ful¬lls the

¬rst criterion, and the fact that the interface is unchanged means that we can plug

the decorated class into any socket which the original class ¬tted. We implement

such a decorator class in AntiThetic.h and AntiThetic.cpp.

Listing 6.5 (AntiThetic.h)

#ifndef ANTITHETIC_H

#define ANTITHETIC_H

#include <Random2.h>

#include <wrapper.h>

class AntiThetic : public RandomBase

{