. 4
( 9)


94 A random numbers class

AntiThetic(const Wrapper<RandomBase>& innerGenerator );

virtual RandomBase* clone() const;

virtual void GetUniforms(MJArray& variates);

virtual void Skip(unsigned long numberOfPaths);

virtual void SetSeed(unsigned long Seed);

virtual void ResetDimensionality(unsigned long

virtual void Reset();
Wrapper<RandomBase> InnerGenerator;

bool OddEven;

MJArray NextVariates;

The decorator class is quite simple. It has an array as a data member to store the
last vector drawn, and a boolean to indicate whether the next draw should be drawn
from the inner generator, or be the anti-thetic of the last draw. A copy of the gen-
erator we are using is stored using the Wrapper template class and cloning, as
usual. Note that we are actually taking a copy of the generator here so that the se-
quence of draws from the original generator will not be affected by drawing from
the anti-thetic generator.

Listing 6.6 (AntiThetic.cpp)

#include <AntiThetic.h>

AntiThetic::AntiThetic(const Wrapper<RandomBase>&
innerGenerator )
: RandomBase(*innerGenerator),
6.5 Anti-thetic sampling via decoration 95

OddEven =true;

RandomBase* AntiThetic::clone() const
return new AntiThetic(*this);

void AntiThetic::GetUniforms(MJArray& variates)
if (OddEven)

for (unsigned long i =0; i < GetDimensionality(); i++)
NextVariates[i] = 1.0-variates[i];

OddEven = false;
variates = NextVariates;

OddEven = true;

void AntiThetic::SetSeed(unsigned long Seed)
OddEven = true;

void AntiThetic::Skip(unsigned long numberOfPaths)
if (numberOfPaths ==0)

if (OddEven)
96 A random numbers class

OddEven = false;

InnerGenerator->Skip(numberOfPaths / 2);

if (numberOfPaths % 2)
MJArray tmp(GetDimensionality());



void AntiThetic::ResetDimensionality(unsigned long



void AntiThetic::Reset()
OddEven =true;

The implementation of the class is quite straightforward. Most of the methods
consist of simply forwarding the request to the inner class, together with book-
keeping for odd and even paths. The main GetUniforms method, gets uniforms
from the inner generator for the odd draws, stores the results, X j , and returns
(1 ’ X 1 , . . . , 1 ’ X n ) for the even draws. Note that

N ’1 (1 ’ x) = ’N ’1 (x), (6.1)

so this will yield the negative of the Gaussian variates if the GetGaussians method
is used, as we wanted.
6.6 Using the random number generator class 97

Note the syntax for initialization in the constructor. We have RandomBase
(*innerGenerator). As innerGenerator is a wrapped pointer, * gives us the
value of the inner object which is a member of some inherited class. However,
we can always treat any inherited class object as a base class object so the call
to RandomBase invokes the base class copy constructor, copying the base class
data in innerGenerator, and thus ensuring that the new object has the correct
dimensionality stored.

6.6 Using the random number generator class
Now that we have a random number generator class, we need to adapt our Monte
Carlo code to work with it. We give an adapted vanilla option pricer in
SimpleMC8.h and SimpleMC8.cpp. The header ¬le declares the new func-

Listing 6.7 (SimpleMC8.h)
#ifndef SIMPLEMC8_H
#define SIMPLEMC8_H

#include <Vanilla3.h>
#include <Parameters.h>
#include <Random2.h>
#include <MCStatistics.h>

void SimpleMonteCarlo6(const VanillaOption& TheOption,
double Spot,
const Parameters& Vol,
const Parameters& r,
unsigned long NumberOfPaths,
StatisticsMC& gatherer,
RandomBase& generator);
We have chosen to take the random number generator in as a non-const reference.
It cannot be a const reference as the act of drawing a random number changes the
generator and is therefore implemented by a non-const method. The effect of this
is that any random numbers drawn inside the function will not be produced outside
the function, but instead the generator will continue where the function left off.
If we wanted the generator to be totally unaffected by what happened inside the
function, we would change the function to take in the object by value instead. Or
alternatively, we could copy the object and pass in the copy to the function, which
98 A random numbers class

would have the same net effect. As usual, we use a reference to the base class in
order to allow the caller to decide how to implement the generator.
The implementation is as follows:

Listing 6.8 (SimpleMC8.cpp)
#include <cmath>
#include <Arrays.h>

// the basic math functions should be in
// namespace std but aren™t in VCPP6
#if !defined(_MSC_VER)
using namespace std;

void SimpleMonteCarlo6(const VanillaOption& TheOption,
double Spot,
const Parameters& Vol,
const Parameters& r,
unsigned long NumberOfPaths,
StatisticsMC& gatherer,
RandomBase& generator)

double Expiry = TheOption.GetExpiry();
double variance = Vol.IntegralSquare(0,Expiry);
double rootVariance = sqrt(variance);
double itoCorrection = -0.5*variance;
double movedSpot = Spot*exp(r.Integral(0,Expiry)
+ itoCorrection);

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

MJArray VariateArray(1);

for (unsigned long i=0; i < NumberOfPaths; i++)
6.6 Using the random number generator class 99

thisSpot = movedSpot*exp(rootVariance*VariateArray[0]);
double thisPayOff = TheOption.OptionPayOff(thisSpot);


We only comment on the new aspects of the routine. We ¬rst reset the generator™s
dimensionality to 1 as pricing a vanilla option is a one-dimensional integral “ we
just need the location of the ¬nal value of spot.
We set up the array in which to store the variate before we set up the main
loop, once and for all. This avoids any dif¬culties with speed in the allocation of
dynamically sized arrays. The GetGaussians method of the generator is used
to write the variates (in this case just one variate, of course) into the array. This
variate is then used as before to compute the ¬nal value of spot.
We give an example of using this routine with anti-thetic sampling in Random-

Listing 6.9 (RandomMain3.cpp)

uses source files
100 A random numbers class

using namespace std;

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);
6.6 Using the random number generator class 101

StatisticsMean gatherer;
ConvergenceTable gathererTwo(gatherer);

RandomParkMiller generator(1);

AntiThetic GenTwo(generator);


vector<vector<double> > results =

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;

We create a Park“Miller random number generator object and then wrap it with
an anti-thetic decorator. This decorated object is then passed into the new Monte
Carlo routine. As usual, the routine is not aware of the fact that the passed-in object
has been decorated but simply uses it in the same way as any other random number
The big difference between our new program and the old ones is that the results
are now compiler-independent. The numbers returned are now precisely the same
under Borland 5.5, Visual C++ 6.0 and MingW 2.95, since we have removed
102 A random numbers class

the dependency on the inbuilt rand() function which previously made our results
compiler-dependent. This gives us an extra robustness test; if our results are not
now compiler-independent we should be worried and ¬nd out why!

6.7 Key points
In this chapter, we developed a random number generator class and saw how anti-
thetic sampling could be implemented via decoration.
• rand is implementation-dependent.
• Results from rand() are not easily reproducible.
• We have to be sure that a random generator is capable of the dimensionality
necessary for a simulation.
• Using a random number class allows us to use decoration.
• The inverse cumulative normal function is the most robust way to turn uniform
variates from the open interval, (0, 1), into Gaussian variates.
• Using a random number class makes it easier to plug in low-discrepancy num-
• Anti-thetic sampling can be implemented via decoration.

6.8 Exercises
Exercise 6.1 For various cases compare convergence of Monte Carlo simulations
with and without anti-thetic sampling.

Exercise 6.2 Obtain another random number generator and ¬t it into the class hi-
erarchy given here. (See [28] or www.boost.org for other generators.)

Exercise 6.3 Code up a low-discrepancy number generator and integrate it into the
classes here. (See [28] or [11].)

An exotics engine and the template pattern

7.1 Introduction
We have now developed quite a few sets of components: random number gener-
ators, parameters classes, pay-off classes, statistics gatherers and a wrapper tem-
plate. Having developed all these components with the objective of reusability, we
examine in this chapter how to put them together to price path-dependent exotic op-
tions. Our objective is to develop a ¬‚exible Monte Carlo pricer for exotic options
which pay off at some future date according to the value of spot on a ¬nite number
of dates. We will work within a deterministic interest rate world, and assume the
Black“Scholes model of stock price evolution.
We assume that our derivative is discrete, i.e. that it depends upon the value of
spot on a discrete set of times. Thus our derivative is associated to a set of times,
t1 , t2 , . . . , tn , and pays at some time T a function f (St1 , . . . , Stn ) of the value of
spot at those times. For example, a one-year Asian call option struck at K with
monthly resets would pay
St j ’ K ,
12 j=1 +

where t j = j/12, at time 1.
More generally, the derivative could possibly pay cash-¬‚ows at more than
one time. For example, a discrete barrier knock-out option could pay an ordi-
nary vanilla pay-off at the time of expiry, and a rebate at the time of knock-out
We do not consider American or Bermudan options here as the techniques in-
volved are quite different. Note, however, that once an exercise strategy has been
chosen the option is really just a path-dependent derivative and so the option can
be evaluated by these techniques for any given ¬xed exercise strategy.

104 An exotics engine and the template pattern

7.2 Identifying components
Before designing our engine let™s identify what it will have to do. For each path,
we generate a discounted pay-off which is then averaged over all paths to obtain a
price. To generate a pay-off for a path, we have to know the path of stock prices at
the relevant times, plug this path into the pay-off function of the option to obtain
the cash-¬‚ows, and then discount these cash-¬‚ows back to the start to obtain the
price for that path.
We can therefore identify four distinct actions:
(i) the generation of the stock price path;
(ii) the generation of cash-¬‚ows given a stock price path;
(iii) the discounting and summing of cash-¬‚ows for a given path;
(iv) the averaging of the prices over all the paths.
We already have a suitable component for the last of these actions: the statistics
gatherer. We will just plug in the class we have already written at the appropriate
point. For the second action, we are purely using the de¬nition of the derivative to
determine what its pay-off is, given a set of stock prices. An obvious component for
our model is therefore a path-dependent exotic option class which will encapsulate
the information which would be written in the term-sheet for the option.
Note that by de¬ning the concept we model by the term-sheet, we make it clear
that the class will not involve interest rates nor knowledge of volatility nor any
aspect of the stock price process. A consequence of this is that the option class
can only ever say what cash-¬‚ows occur and when; it cannot say anything about
their discounted values because that would require knowledge of interest rates.
Note the general point here that, by de¬ning the concept in real-world terms, it
becomes very easy to decide what should and should not be contained in the class.
Another consequence is that the component is more easily reusable; if we decide
to do jump-diffusion pricing or stochastic interest rates, this class will be reusable
without modi¬cation, and that would not be the case if we had included information
about interest rates or volatility.
There is more than one way to handle the two remaining tasks: path generation
and cash-¬‚ow accounting. The latter will be the same for any deterministic interest-
rate model and it is therefore natural to include it as part of our main engine class.
We can require path generation to be an input to our main class, and therefore
de¬ne it in terms of its own class hierarchy, or via a virtual method of the base
class. A third option, which we do not explore, would be to make it a template
parameter, which would avoid the small overhead of a virtual function call.
The option we will pursue here is to make path generation a virtual method of
the base class. This is an example of the template design pattern, which should not
be confused with templatized code. The idea here is that the base class sets up a
structure with methods that control everything “ in this case it would be methods
7.3 Communication between the components 105

to run the simulation and account for each path “ though the main work of actually
generating the path is not de¬ned in the base class, but is instead called via a pure
virtual function. This pure virtual function must therefore be de¬ned in an inherited
class. Thus as in templatized code, the code is set up more to de¬ne a structure than
to do the real work which is coded elsewhere.

7.3 Communication between the components
Having identi¬ed what our components will be, we still need to assess what infor-
mation has to be passed between them, and we need to decide how to carry out the
The product will take in a vector of spot values for its relevant times and spit out
the cash-¬‚ows generated. A couple of immediate consequences of this are that there
has to be a mechanism for the product to tell the path generator for which times it
needs the value of spot, and that we need to decide how to de¬ne a cash-¬‚ow object.
To deal with the ¬rst of these, we include a method GetLookAtTimes; this
passes back an array of times that are relevant to the pay-off function of the prod-
uct. For the de¬nition of cash-¬‚ow objects, there are a couple of options. The ¬rst
obvious approach is to simply make a cash-¬‚ow a pair of doubles which de¬ne
the amount and the time of the cash-¬‚ow. This approach has the disadvantage that
if we have a complicated term structure of interest rates, the action of computing
the discount factor for the cash-¬‚ow time may be slow, and with the product being
allowed to pass back an arbitrary time, this discounting will have to be done on
every path. Whilst one could cache already-computed times, there is then still the
problem that searching the cache for the already-computed times will take time.
In practice, many products can only pay off at one time. This means that it would
be better to pre-compute the discount factor for that time. However, we would still
need to know in advance what that time is. We therefore require our product to
have another method, PossibleCashFlowTimes, which returns an array de¬ning
the possible times. As the engine will know all the possible times in advance we
can return a cash-¬‚ow as a pair: an index and an amount. The index is now an
unsigned long instead of a double. The engine will now precompute all the
discount factors and then simply use the index to look up an array to get the dis-
counting for any given cash-¬‚ow.
We still have to decide the syntax for the main method CashFlows. The method
takes in an array de¬ning spot values, and returns cash-¬‚ows. As we allow the
possibility of more than one cash-¬‚ow, we must use a container to pass them back.
We use the STL vector class. Whilst it would be tempting to make the return type
of the class vector<CashFlow>, this would have timing disadvantages. We would
have to create a new vector every time the method was called, and this could be time
consuming because any dynamically sized container class must involve memory
106 An exotics engine and the template pattern

allocation. We therefore take an argument of type vector<CashFlow>& into which
we write the cash-¬‚ows.
We still have the issue that the vector will need to be of the correct size. One
solution is for the method to resize it as necessary but this could be time consuming.
First, resizing can involve memory allocation though this is not a huge issue since
the memory allocated for an STL vector never shrinks so if the same vector is
used every time it will rapidly grab enough memory and then will need no more.
Second, some implementations of the STL explicitly destroy all the objects in the
vector during a resize, which means that every resize involves looping, and is
therefore unnecessarily slow even when no memory allocation is necessary.
The solution we adopt is to tell the outside engine how big the vector has to be,
and then each time the method is called, to return an unsigned long saying how
many cash-¬‚ows have been generated. Thus we have two pure virtual methods:
virtual unsigned long MaxNumberOfCashFlows() const=0;
virtual unsigned long CashFlows(const MJArray& SpotValues,
GeneratedFlows) const=0;
So in summary our objects will communicate as follows:
(i) The path generator asks the product what times it needs spot for, and it passes
back an array.
(ii) The accounting part of the engine asks the product what cash-¬‚ow times are
possible, and it passes back an array. The engine then computes all the possi-
ble discount factors.
(iii) The accounting part of the engine asks the product the maximum number of
cash ¬‚ows it can generate, and sets up a vector of that size.
(iv) For each path, the engine gets an array of spot values from the path generator.
(v) The array of spot values is passed into the product, which passes back the
number of cash-¬‚ows, and puts their values into the vector.
(vi) The cash-¬‚ows are discounted appropriately and summed, and the total value
is passed into the statistics gatherer.
(vii) After all the looping is done, the ¬nal results are obtained from the statistics

7.4 The base classes
Having discussed in previous sections what classes will be needed and how they
should communicate, in this section we give the implementations of the base
In PathDependent.h, we de¬ne the CashFlow and the PathDependent classes
which give our path-dependent exotic option.
7.4 The base classes 107

Listing 7.1 (PathDependent.h)

#include <Arrays.h>
#include <vector>

class CashFlow
double Amount;
unsigned long TimeIndex;

CashFlow(unsigned long TimeIndex_=0UL, double Amount_=0.0)
: TimeIndex(TimeIndex_),

class PathDependent
PathDependent(const MJArray& LookAtTimes);

const MJArray& GetLookAtTimes() const;

virtual unsigned long MaxNumberOfCashFlows() const=0;
virtual MJArray PossibleCashFlowTimes() const=0;
virtual unsigned long CashFlows(const MJArray& SpotValues,
GeneratedFlows) const=0;
virtual PathDependent* clone() const=0;

virtual ˜PathDependent(){}
MJArray LookAtTimes;

The CashFlow class is really just a struct as it has public data members. Note
that we ensure that it has a default constructor by giving the constructor default
arguments, this is necessary in order to use it with STL containers which need
108 An exotics engine and the template pattern

a default constructor for certain operations such as creating a vector of arbitrary
The base class for PathDependent really does not do much except de¬ne the
interface. We have made the base class store the LookAtTimes as every possi-
ble product will need these times, and provided the method GetLookAtTimes to
obtain them. As usual, we include a clone method for virtual copy construction,
and a virtual destructor to make sure that there are no memory leaks arising from
destroying base class objects instead of inherited ones.
The source code is suitably short:

Listing 7.2 (PathDependent.cpp)
#include <PathDependent.h>

PathDependent::PathDependent(const MJArray& LookAtTimes_)
: LookAtTimes(LookAtTimes_)

const MJArray& PathDependent::GetLookAtTimes() const
return LookAtTimes;
There is a bit more to the base class for the engine as it will actually handle the
accounting for the cash-¬‚ows.

Listing 7.3 (ExoticEngine.h)
#include <wrapper.h>
#include <Parameters.h>
#include <PathDependent.h>
#include <MCStatistics.h>
#include <vector>
class ExoticEngine
ExoticEngine(const Wrapper<PathDependent>&
The Product_, const Parameters& r_);

virtual void GetOnePath(MJArray& SpotValues)=0;
7.4 The base classes 109

void DoSimulation(StatisticsMC& TheGatherer,
unsigned long NumberOfPaths);
virtual ˜ExoticEngine(){}
double DoOnePath(const MJArray& SpotValues) const;

Wrapper<PathDependent> TheProduct;
Parameters r;
MJArray Discounts;
mutable std::vector<CashFlow> TheseCashFlows;

The engine has four data members. The product is stored using the Wrapper tem-
plate as we do not know its type. The interest rates are stored using the Parameters
class which will allow us variable ones if we so desire. We also delegate computa-
tion of integrals to the Parameters class, and not have to worry about them here.
We have an array Discounts, which will be used to store the discount factors
in order for the possible cash-¬‚ow times. Finally we have a mutable data mem-
ber TheseCashFlows. This means that it can change value inside const member
functions. The idea is that the variable is not really a data member, but instead it is
a workspace variable: this it is faster to declare once and for all in the class de¬ni-
tion. Remember that creating and destroying containers can be time-consuming so
we design the class so that the vector is created once and for all at the beginning.
Note that we split our main method; it has two auxiliary methods, DoOnePath
and GetOnePath. The second of these is pure virtual and therefore will be de¬ned
in an inherited class which will involve a choice of stochastic process and model.
Note that this method is not constant as we will want a different set of spot values
every time, and so it will necessarily change something about the state of the object.
The other of the methods does everything necessary for one path given the spot
values. This is const as turning spot values into prices is a purely functional action
with no underlying changes. Both these methods pass arrays by reference in order
to avoid any memory allocation. Note the implicit assumption that the array passed
into GetOnePath is of the correct size.
The source code for implementing the base class is fairly simple and straightfor-
ward as all the hard work has been hived off into auxiliary classes.

Listing 7.4 (ExoticEngine.cpp)

#include <ExoticEngine.h>
#include <cmath>
110 An exotics engine and the template pattern

ExoticEngine::ExoticEngine(const Wrapper<PathDependent>&
const Parameters& r_)
for (unsigned long i=0; i < Discounts.size(); i++)
Discounts[i] = exp(-r.Integral(0.0, Discounts[i]));


void ExoticEngine::DoSimulation(StatisticsMC& TheGatherer,
unsigned long NumberOfPaths)
MJArray SpotValues(TheProduct->GetLookAtTimes().size());

double thisValue;

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


double ExoticEngine::DoOnePath(const MJArray&
SpotValues) const
unsigned long NumberFlows =
double Value=0.0;
7.5 A Black“Scholes path generation engine 111

for (unsigned i =0; i < NumberFlows; ++i)
Value += TheseCashFlows[i].Amount *

return Value;
The constructor stores the inputs, computes the discount factors necessary, and
makes sure the cash-¬‚ows vector is of the correct size. The DoSimulation method
loops through all the paths, calling GetOnePath to get the array of spot value and
then passes them into DoOnePath to get the value for that set of spot values. This
value is then dumped into the statistics gatherer.
DoOnePath is only slightly more complicated. The array of spot values is passed
into the product to get the cash-¬‚ows. These cash-¬‚ows are then looped over and
discounted appropriately. The discounting is simpli¬ed by using the precomputed
discount factors.
We have now set up the structure for pricing path-dependent exotic derivatives
but we still have to actually de¬ne the classes which will do the path generation
and de¬ne the products.

7.5 A Black“Scholes path generation engine
The Black“Scholes engine will produce paths from the risk-neutral Black“Scholes
process. The paths will be an array of spot values at the times speci¬ed by the
product. We allow the possibility of variable interest rates and dividend rates, as
well as variable but deterministic volatility. The stock price therefore follows the
d St = (r (t) ’ d(t))St dt + σ (t)St dWt , (7.1)
with S0 given. To simulate this process at times t0 , t1 , . . . , tn’1 , we need n
independent N (0, 1) variates W j and we set
t0 t0
log St0 = log S0 + r (s) ’ d(s) ’ σ (s)2 ds + σ (s)2 dsW0 , (7.2)
0 0

and put

tj tj
log St j = log St j’1 + r (s) ’ d(s) ’ σ (s)2 ds + σ (s)2 dsW j . (7.3)
t j’1 t j’1
112 An exotics engine and the template pattern

We implement this procedure in ExoticBSEngine.h and ExoticBS-

Listing 7.5 (ExoticBSEngine.h)
#include <ExoticEngine.h>
#include <Random2.h>

class ExoticBSEngine : public ExoticEngine
ExoticBSEngine(const Wrapper<PathDependent>& TheProduct_,
const Parameters& R_,
const Parameters& D_,
const Parameters& Vol_,
const Wrapper<RandomBase>& TheGenerator_,
double Spot_);

virtual void GetOnePath(MJArray& SpotValues);
virtual ˜ExoticBSEngine(){}

Wrapper<RandomBase> TheGenerator;
MJArray Drifts;
MJArray StandardDeviations;
double LogSpot;
unsigned long NumberOfTimes;
MJArray Variates;

Listing 7.6 (ExoticBSEngine.cpp)
#include <ExoticBSEngine.h>
#include <cmath>

void ExoticBSEngine::GetOnePath(MJArray& SpotValues)
7.5 A Black“Scholes path generation engine 113

double CurrentLogSpot = LogSpot;

for (unsigned long j=0; j < NumberOfTimes; j++)
CurrentLogSpot += Drifts[j];
CurrentLogSpot += StandardDeviations[j]*Variates[j];
SpotValues[j] = exp(CurrentLogSpot);


ExoticBSEngine::ExoticBSEngine(const Wrapper<PathDependent>&
const Parameters& R_,
const Parameters& D_,
const Parameters& Vol_,
const Wrapper<RandomBase>&
double Spot_)
MJArray Times(TheProduct_->GetLookAtTimes());
NumberOfTimes = Times.size();


double Variance = Vol_.IntegralSquare(0,Times[0]);

Drifts[0] = R_.Integral(0.0,Times[0])
- D_.Integral(0.0,Times[0]) - 0.5 * Variance;
StandardDeviations[0] = sqrt(Variance);

for (unsigned long j=1; j < NumberOfTimes; ++j)
114 An exotics engine and the template pattern

double thisVariance =
Drifts[j] = R_.Integral(Times[j-1],Times[j])
- D_.Integral(Times[j-1],Times[j])
- 0.5 * thisVariance;
StandardDeviations[j] = sqrt(thisVariance);

LogSpot = log(Spot_);

The integrals and square-roots are the same for every path and so can be precom-
puted. The constructor therefore gets the times from the product, and uses them to
compute the integrals of the drifts and the standard deviations which are stored as
data members. Note that the class does not bother to store the times as it is only the
constructor which needs to know what they are. In any case, the product is passed
up to the base class and it could be retrieved from there if it were necessary.
The generation will of course require a random number generator and we pass
in a wrapped RandomBase object to allow us to plug in any one we want without
having to do any explicit memory handling. We have a data member Variates
so that the array can be de¬ned once and for all at the beginning: once again this
is with the objective of avoiding unnecessary creation and deletion of objects. We
store the log of the initial value of spot as this is the most convenient for carrying
out the path generation.
As we have done a lot of precomputation in the constructor, the routine to actu-
ally generate a path is fairly simple. We simply get the variates from the generator
and loop through the times. For each time, we add the integrated drift to the log,
and then add the product of the random number and the standard deviation. To
minimize the number of calls to log and exp, we keep track of the log of the spot at
all times, and convert into spot values as necessary. We thus have NumberOfTimes
calls to exp each path and no calls to log. As we will have to exponentiate to change
our Gaussian into a log-normal variate at some point, this appears to be optimal for
this design. If we were really worried that too much time was being spent on com-
puting exponentials, one solution would be to change the design and pass the log of
the values of spot back, and then pass these log values into the product. The prod-
uct would then have the obligation to exponentiate them if necessary. For certain
products such as a geometric Asian option this might well be faster as it would only
involve one exponentiation instead of many. The main downside would be that for
7.6 An arithmetic Asian option 115

certain processes, such as a normal process or displaced diffusion, one might end
up having to take unnecessary logs.

7.6 An arithmetic Asian option
Before we can run our engine, we need one last thing, namely a concrete product
to put in it. One simple example is an arithmetic Asian option. Rather than de¬ne a
different class for each sort of pay-off, we use the already developed PayOff class
as a data member.
The header ¬le for the class is quite simple:

Listing 7.7 (PathDependentAsian.h)

#include <PathDependent.h>
#include <PayOffBridge.h>

class PathDependentAsian : public PathDependent
PathDependentAsian(const MJArray& LookAtTimes_,
double DeliveryTime_,
const PayOffBridge& ThePayOff_);

virtual unsigned long MaxNumberOfCashFlows() const;
virtual MJArray PossibleCashFlowTimes() const;
virtual unsigned long CashFlows(const MJArray& SpotValues,
std::vector<CashFlow>& GeneratedFlows) const;
virtual ˜PathDependentAsian(){}
virtual PathDependent* clone() const;
double DeliveryTime;
PayOffBridge ThePayOff;
unsigned long NumberOfTimes;
The methods de¬ned are just the ones required by the base class. We pass in the
averaging times as an array and we provide a separate delivery time to allow for the
possibility that the pay-off occurs at some time after the last averaging date. Note
116 An exotics engine and the template pattern

that the use of PayOffBridge class means that the memory handling is handled
internally, and this class does not need to worry about assignment, copying and
The source ¬le is fairly simple too.

Listing 7.8 (PathDependentAsian.cpp)
#include <PathDependentAsian.h>

PathDependentAsian::PathDependentAsian(const MJArray&
double DeliveryTime_,
const PayOffBridge&ThePayOff_)

unsigned long PathDependentAsian::MaxNumberOfCashFlows() const
return 1UL;

MJArray PathDependentAsian::PossibleCashFlowTimes() const
MJArray tmp(1UL);
tmp[0] = DeliveryTime;
return tmp;

unsigned long PathDependentAsian::CashFlows(const MJArray&
std::vector<CashFlow>& GeneratedFlows) const
double sum = SpotValues.sum();
double mean = sum/NumberOfTimes;

GeneratedFlows[0].TimeIndex = 0UL;
7.7 Putting it all together 117

GeneratedFlows[0].Amount = ThePayOff(mean);

return 1UL;

PathDependent* PathDependentAsian::clone() const
return new PathDependentAsian(*this);

Note that our option only ever returns one cash-¬‚ow so the maximum number
of cash-¬‚ows is 1. It only ever generates cash-¬‚ows at the delivery time so the
PossibleCashFlowTimes method is straightforward too. The CashFlows method
takes the spot values, sums them, divides by the number of them and calls ThePay-
Off to ¬nd out what the pay-off is. The answer is then written into the Generated-
Flows array and we are done.

7.7 Putting it all together
We now have everything we need to price an Asian option. We give an example of
a simple interface program in EquityFXMain.cpp.

Listing 7.9 (EquityFXMain.cpp)

uses source files
118 An exotics engine and the template pattern

using namespace std;
int main()
double Expiry;
double Strike;
double Spot;
double Vol;
double r;
double d;
unsigned long NumberOfPaths;
unsigned NumberOfDates;

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 << "\nd\n";
cin >> d;

cout << "Number of dates\n";
cin >> NumberOfDates;
7.7 Putting it all together 119

cout << "\nNumber of paths\n";
cin >> NumberOfPaths;

PayOffCall thePayOff(Strike);

MJArray times(NumberOfDates);

for (unsigned long i=0; i < NumberOfDates; i++)
times[i] = (i+1.0)*Expiry/NumberOfDates;

ParametersConstant VolParam(Vol);
ParametersConstant rParam(r);
ParametersConstant dParam(d);

PathDependentAsian theOption(times, Expiry, thePayOff);

StatisticsMean gatherer;
ConvergenceTable gathererTwo(gatherer);

RandomParkMiller generator(NumberOfDates);

AntiThetic GenTwo(generator);

ExoticBSEngine theEngine(theOption, rParam, dParam,
VolParam, GenTwo, Spot);

theEngine.DoSimulation(gathererTwo, NumberOfPaths);

vector<vector<double> > results =

cout <<"\nFor the Asian 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";
120 An exotics engine and the template pattern


double tmp;
cin >> tmp;

return 0;

7.8 Key points
In this chapter, we saw how we can put the ideas developed in the previous chapters
together to build a pricer for exotic options.
• An important part of the design process is identifying the necessary components
and specifying how they talk to each other.
• The template pattern involves deferring the implementation of an important part
of an algorithm to an inherited class.
• If an option class knows nothing that is not speci¬ed in the term-sheet then it is
much easier to reuse.
• We can reuse the PayOff class to simplify the coding of our more complicated
path-dependent derivatives.

7.9 Exercises
Exercise 7.1 Write a class to do geometric Asian options.

Exercise 7.2 Write a class to do discrete knock-out options that pay a rebate at the
time of rebate.

Exercise 7.3 Rewrite the classes here so that they pass the logs of spot values
around instead of the spot values. Show that the discrete barrier option and the
geometric Asian need fewer exponentiations.

Exercise 7.4 Implement an engine for pricing when the spot price is normal instead
of log-normal.

Exercise 7.5 Write a class that pays the difference in pay-offs of two arbitrary
path-dependent derivatives.


8.1 Introduction
We have studied Monte Carlo code in detail: now we examine how we can apply
similar techniques to pricing on trees. Before we can start designing the code, we
need to ¬x the underlying mathematics. The point of view we adopt is that a tree is
a method of approximating the risk-neutral expectation. In particular, we assume
that we are pricing a derivative on a stock following geometric Brownian motion
with a constant volatility σ . We let the continuously compounding interest rate be
r and the continuous dividend rate be d. The dynamics for the stock in the risk-
neutral measure are therefore
d S = (r ’ d)Sdt + σ Sd Wt . (8.1)
The value of a European option with expiry T is then
e’r T E(C(S, T )), (8.2)
where C(S, T ) is the ¬nal pay-off.
When we price on a binomial tree, we divide time into steps and across each step
we assume that the underlying Brownian motion can only move a ¬xed amount up
or down. The dynamics of the stock price under geometric Brownian motion are
such that
St = S0 e(r ’d’ 2 σ
2 )t+σ W
. (8.3)

We wish to discretize Wt . We have N steps to get from 0 to T . Each time step is
therefore of length T /N . Across step l, we need to approximate
W(l+1)T /N ’ WlT /N = N (0, 1). (8.4)
There is only one random variable taking precisely two values which has the same
mean and variance as N (0, 1), and this variable takes ±1 with probability 1/2. We

122 Trees

therefore take a set of N independent random variables X j with this distribution,
and approximate WlT /N by
Yl = X j. (8.5)
N j=1

The approximation for SlT /N is then

S0 e(r ’d’ 2 σ
2 )lT /N +σ Y

Note the crucial point here that since the value of St does not depend on the path
of Wt but solely upon its value at time t, it is only the value of Yl that matters not
the value of each individual X j . This is crucial because it means that our tree is
recombining; it does not matter whether we go down then up, or up then down.
This is not the case if we allow variable volatility, which is why we have assumed
its constancy.
The nature of martingale pricing means that the value at a given time-step is
equal to the discounted value at the next time-step, thus if we let Sl,k be the value
of the stock at time T l/N if Yl is k, we have that

C(Sl,k , T l/N ) = e’r T /N E(Sl+1 (Yl+1 |Yl = k)),

1 ’r T /N (r ’d’ 1 σ 2 )T /N +σ T /N
=e , (l + 1)T /N
C Sl,k e 2
2 √
(r ’d’ 1 σ 2 )T /N ’σ T /N
+ C Sl,k e , (l + 1)T /N . (8.6)

Note that we are not doing true martingale pricing in the discrete world in that we
are not adjusting probabilities to make sure the assets are discrete martingales; we
are instead approximating the continuous martingales with discrete random vari-
ables which are almost, but not quite, martingales.
What does all this buy us? The value of C(S N ,k ) is easy to compute for any value
of k: we just plug the stock price into the ¬nal pay-off. Since we have a formula
that expresses the value at step l into that at step l + 1, we can now just backwards
iterate through the tree in order to get the price at time zero.
However, the purpose of a tree is not to price a European option; there are lots of
better ways of doing that, including analytical solutions and numerical integration.
The reason trees were introduced was that they give an effective method for pricing
American options. The analysis for an American option is similar except that the
value at a point in the tree is the maximum of the exercise value at that point and
the discounted expected value at the next time. This corresponds to the optimal
strategy of exercise if and only if exercise gives more money than not exercising.
Our algorithm for pricing an American option is therefore as follows:
8.2 The design 123

(i) Create an array of ¬nal spot values which are of the form

(r ’d’ 1 σ 2 )T +σ T /N j
S0 e 2

where j ranges from ’N to N .
(ii) For each of these spot values evaluate the pay-off and store it.
(iii) At the previous time-slice compute the possible values of spot: these will be
of the form

S0 e(r ’d’ 2 σ
2 )(N ’1)T /N +σ T /N j

where j ranges from ’(N ’ 1) to N ’ 1.
(iv) For each of these values of spot, compute the pay-off and take the maximum
with the discounted pay-off of the two possible values of spot at the next time.
(v) Repeat 3,4 until time zero is reached.
What else could we price on a tree? We could do a barrier option or an Ameri-
can option that could only be exercised within certain time periods. For a knock-out
barrier option, the procedure would be the same as for the European, except that the
value at a point in the tree would be zero if it lay behind the barrier. For an Amer-
ican option with limited early exercise the procedure would be the same again,
except that we would only take the maximum at times at which early exercise was
allowed. So in each case, when we change the option, all that alters is the rule for
updating the value at a point in the tree.
Note that in our formulation, we have not used any no-arbitrage arguments. The
reason is that we have implicitly assumed that the no-arbitrage arguments have
been done in the continuous case before any discretization has been carried out.
This means that the completeness property of a binomial tree, i.e. that it generates
a unique no-arbitrage price, is not relevant. In particular, we could replace X j by
any random variable with mean 0 and variance 1. If we use X j which takes values
√ √
’ 2, 0, 2 with probabilities 0.25, 0.5 and 0.25 respectively, then we obtain a
trinomial tree on which we could carry out a very similar analysis and algorithm.
We could in fact go much further and have as many values as we wanted as long as
we took care to make sure that the tree still recombined.

8.2 The design
Having re-examined the mathematics and the algorithms, we are now in a position
to think about the design. Here are some concepts that our discussion has thrown
• the discretization;
• the ¬nal pay-off of the option;
124 Trees

• the rule for deciding the value of an option at a point in the tree given spot and
the discounted future value of the option.
The ¬rst of these concepts decides the shape of the tree, whereas the second
and third are properties of the option. There is thus an obvious orthogonalization:
we have a tree class which handles the discretization, and a second class which
deals with the ¬nal pay-off and the rule at previous times. In fact, we have already
developed a class, PayOff, to encapsulate vanilla option pay-offs and it is ideal for
reuse here.
There are a number of ways we could approach the second class. We could in-
herit from PayOffBridged since we could view our class as adding structure to an
existing class. Whilst this would work in code, I personally dislike it as an option
being priced on a tree is not a type of pay-off, and so the inheritance is not express-
ing an is a relationship. Another approach might be simply to de¬ne a second class
to do the calculation rule, and plug both the ¬nal pay-off and the calculation rule
into the tree. Since for American options the ¬nal pay-off is generally relevant at
all times, such an approach seems sub-optimal as it might require two copies of the
Ultimately, the pay-off is an aspect of the option, and it therefore makes more
sense to de¬ne it as data member of the class which can referenced via a ¬nal pay-
off method. Thus we de¬ne options on trees via an abstract base class which has
three de¬ning methods:
• FinalTime indicates the expiry time of the option;
• FinalPayOff gives the ¬nal pay-off as a function of spot;
• PreFinalValue gives the value at a point in the tree as a function of spot, time
and the discounted future value of the option.
Note that by de¬ning the option class in this fashion, we have not allowed it to
know anything about interest rates nor the process of the underlying. This means it
can be used in any tree-like structure provided the structure is always in a position
to let it know its own discounted future value. Note the difference here between the
option classes we are de¬ning here and those we de¬ned for Monte Carlo: whilst
we are trying to encapsulate similar concepts, the difference is in the information
we are able to feed in at a given time. For Monte Carlo, the entire history of spot
is easy but the discounted future value of an option is hard, whereas on a tree the
discounted future value is easy but the history is hard. However, both classes have
easy access to the pay-off which means we are able to share the pay-off class.
Our other concept is the tree itself. The tree really has two aspects: the placement
of the nodes of the tree and the computing of the option value at each of the nodes.
Whilst we could further orthogonalize and de¬ne separate classes for each of these,
we write a single class to do the binomial tree which takes in an option as an
argument. An important point to note is that as the placement of the nodes does
8.3 The TreeProduct class 125

not depend upon the option, we can save ourselves time if we want to price several
options by placing the nodes once and then pricing all the options on the same tree.
Given this fact, we design our tree class in such a way that the tree is built once,
and then any option can be valued on the tree via a separate method.

8.3 The TreeProduct class
As we have decided to model the tree and the product separately, we develop a
class hierarchy for the products we can value on trees. As usual, we use an ab-
stract base class to de¬ne an interface. We de¬ne the class, TreeProduct, in

Listing 8.1 (TreeProducts.h)


class TreeProduct
TreeProduct(double FinalTime_);
virtual double FinalPayOff(double Spot) const=0;
virtual double PreFinalValue(double Spot,
double Time,
double DiscountedFutureValue)
virtual ˜TreeProduct(){}
virtual TreeProduct* clone() const=0;
double GetFinalTime() const;

double FinalTime;

The only data member for the base class is FinalTime, and we provide a
GetFinalTime to allow its value to be read. Note that this places the constraint
on our products that they actually have a time of expiry! Thus we are implicitly
disallowing perpetual options, however this is a good thing as it is not clear how
to go about valuing such an option using a tree. Ultimately, one would have to
approximate using a product with a ¬nite expiry and it is probably better to do so
explicitly than implicitly.
126 Trees

We provide the usual clone method and a virtual destructor to allow virtual
copying, and to ensure the absence of memory leaks after virtual copying. The
remaining methods are pure virtual and specify the value of the product at expiry
and at previous times as we discussed above.
As most of the base class is abstract the source ¬le is short:

Listing 8.2 (TreeProducts.cpp)
#include <TreeProducts.h>

TreeProduct::TreeProduct(double FinalTime_)
: FinalTime(FinalTime_)

double TreeProduct::GetFinalTime() const
return FinalTime;
We provide below two concrete implementations of tree products: the European
option and the American option. As we specify the pay-off using the PayOff-
Bridged class, we do not need to write separate classes for calls and puts. The
header ¬les are straightforward, the only changes from the base class being the
addition of a data member to specify the pay-off and the fact that virtual methods
are now concrete instead of abstract. Note that everything to do with the expiry
time is taken care of by the base class; the constructors need only pass it on to the
base class.

Listing 8.3 (TreeAmerican.h)

#include <TreeProducts.h>
#include <PayOffBridge.h>

class TreeAmerican : public TreeProduct

TreeAmerican(double FinalTime,
const PayOffBridge& ThePayOff_);
8.3 The TreeProduct class 127

virtual TreeProduct* clone() const;
virtual double FinalPayOff(double Spot) const;
virtual double PreFinalValue(double Spot,


. 4
( 9)