<<

. 3
( 9)



>>


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
{

<<

. 3
( 9)



>>