ńņš. 4 |

public:

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

NewDimensionality);

virtual void Reset();

private:

Wrapper<RandomBase> InnerGenerator;

bool OddEven;

MJArray NextVariates;

};

#endif

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),

InnerGenerator(innerGenerator)

{

6.5 Anti-thetic sampling via decoration 95

InnerGenerator->Reset();

OddEven =true;

NextVariates.resize(GetDimensionality());

}

RandomBase* AntiThetic::clone() const

{

return new AntiThetic(*this);

}

void AntiThetic::GetUniforms(MJArray& variates)

{

if (OddEven)

{

InnerGenerator->GetUniforms(variates);

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

NextVariates[i] = 1.0-variates[i];

OddEven = false;

}

else

{

variates = NextVariates;

OddEven = true;

}

}

void AntiThetic::SetSeed(unsigned long Seed)

{

InnerGenerator->SetSeed(Seed);

OddEven = true;

}

void AntiThetic::Skip(unsigned long numberOfPaths)

{

if (numberOfPaths ==0)

return;

if (OddEven)

96 A random numbers class

{

OddEven = false;

numberOfPaths--;

}

InnerGenerator->Skip(numberOfPaths / 2);

if (numberOfPaths % 2)

{

MJArray tmp(GetDimensionality());

GetUniforms(tmp);

}

}

void AntiThetic::ResetDimensionality(unsigned long

NewDimensionality)

{

RandomBase::ResetDimensionality(NewDimensionality);

NextVariates.resize(NewDimensionality);

InnerGenerator->ResetDimensionality(NewDimensionality);

}

void AntiThetic::Reset()

{

InnerGenerator->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-

tion.

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);

#endif

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<SimpleMC8.h>

#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;

#endif

void SimpleMonteCarlo6(const VanillaOption& TheOption,

double Spot,

const Parameters& Vol,

const Parameters& r,

unsigned long NumberOfPaths,

StatisticsMC& gatherer,

RandomBase& generator)

{

generator.ResetDimensionality(1);

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

generator.GetGaussians(VariateArray);

thisSpot = movedSpot*exp(rootVariance*VariateArray[0]);

double thisPayOff = TheOption.OptionPayOff(thisSpot);

gatherer.DumpOneResult(thisPayOff*discounting);

}

return;

}

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-

Main3.cpp.

Listing 6.9 (RandomMain3.cpp)

/*

uses source files

AntiThetic.cpp

Arrays.cpp,

ConvergenceTable.cpp,

MCStatistics.cpp

Normals.cpp

Parameters.cpp,

ParkMiller.cpp

PayOff3.cpp,

PayOffBridge.cpp,

Random2.cpp,

SimpleMC8.cpp

Vanilla3.cpp,

*/

#include<SimpleMC8.h>

#include<ParkMiller.h>

#include<iostream>

100 A random numbers class

using namespace std;

#include<Vanilla3.h>

#include<MCStatistics.h>

#include<ConvergenceTable.h>

#include<AntiThetic.h>

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);

SimpleMonteCarlo6(theOption,

Spot,

VolParam,

rParam,

NumberOfPaths,

gathererTwo,

GenTwo);

vector<vector<double> > results =

gathererTwo.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;

}

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

generator.

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-

bers.

ā¢ 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].)

7

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

12

1

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

otherwise.

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.

103

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

communication.

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,

std::vector<CashFlow>&

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

gatherer.

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

classes.

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)

#ifndef PATH_DEPENDENT_H

#define PATH_DEPENDENT_H

#include <Arrays.h>

#include <vector>

class CashFlow

{

public:

double Amount;

unsigned long TimeIndex;

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

: TimeIndex(TimeIndex_),

Amount(Amount_){};

};

class PathDependent

{

public:

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,

std::vector<CashFlow>&

GeneratedFlows) const=0;

virtual PathDependent* clone() const=0;

virtual ˜PathDependent(){}

private:

MJArray LookAtTimes;

};

#endif

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

size.

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)

#ifndef EXOTIC_ENGINE_H

#define EXOTIC_ENGINE_H

#include <wrapper.h>

#include <Parameters.h>

#include <PathDependent.h>

#include <MCStatistics.h>

#include <vector>

class ExoticEngine

{

public:

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;

private:

Wrapper<PathDependent> TheProduct;

Parameters r;

MJArray Discounts;

mutable std::vector<CashFlow> TheseCashFlows;

};

#endif

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>&

TheProduct_,

const Parameters& r_)

:

TheProduct(TheProduct_),

r(r_),

Discounts(TheProduct_->PossibleCashFlowTimes())

{

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

Discounts[i] = exp(-r.Integral(0.0, Discounts[i]));

TheseCashFlows.resize(TheProduct->MaxNumberOfCashFlows());

}

void ExoticEngine::DoSimulation(StatisticsMC& TheGatherer,

unsigned long NumberOfPaths)

{

MJArray SpotValues(TheProduct->GetLookAtTimes().size());

TheseCashFlows.resize(TheProduct->MaxNumberOfCashFlows());

double thisValue;

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

{

GetOnePath(SpotValues);

thisValue = DoOnePath(SpotValues);

TheGatherer.DumpOneResult(thisValue);

}

return;

}

double ExoticEngine::DoOnePath(const MJArray&

SpotValues) const

{

unsigned long NumberFlows =

TheProduct->CashFlows(SpotValues,

TheseCashFlows);

double Value=0.0;

7.5 A Blackā“Scholes path generation engine 111

for (unsigned i =0; i < NumberFlows; ++i)

Value += TheseCashFlows[i].Amount *

Discounts[TheseCashFlows[i].TimeIndex];

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

process

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

1

log St0 = log S0 + r (s) ā’ d(s) ā’ Ļ (s)2 ds + Ļ (s)2 dsW0 , (7.2)

2

0 0

and put

tj tj

1

log St j = log St jā’1 + r (s) ā’ d(s) ā’ Ļ (s)2 ds + Ļ (s)2 dsW j . (7.3)

2

t jā’1 t jā’1

112 An exotics engine and the template pattern

We implement this procedure in ExoticBSEngine.h and ExoticBS-

Engine.cpp.

Listing 7.5 (ExoticBSEngine.h)

#ifndef EXOTIC_BS_ENGINE_H

#define EXOTIC_BS_ENGINE_H

#include <ExoticEngine.h>

#include <Random2.h>

class ExoticBSEngine : public ExoticEngine

{

public:

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(){}

private:

Wrapper<RandomBase> TheGenerator;

MJArray Drifts;

MJArray StandardDeviations;

double LogSpot;

unsigned long NumberOfTimes;

MJArray Variates;

};

#endif

Listing 7.6 (ExoticBSEngine.cpp)

#include <ExoticBSEngine.h>

#include <cmath>

void ExoticBSEngine::GetOnePath(MJArray& SpotValues)

{

TheGenerator->GetGaussians(Variates);

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);

}

return;

}

ExoticBSEngine::ExoticBSEngine(const Wrapper<PathDependent>&

TheProduct_,

const Parameters& R_,

const Parameters& D_,

const Parameters& Vol_,

const Wrapper<RandomBase>&

TheGenerator_,

double Spot_)

:

ExoticEngine(TheProduct_,R_),

TheGenerator(TheGenerator_)

{

MJArray Times(TheProduct_->GetLookAtTimes());

NumberOfTimes = Times.size();

TheGenerator->ResetDimensionality(NumberOfTimes);

Drifts.resize(NumberOfTimes);

StandardDeviations.resize(NumberOfTimes);

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 =

Vol_.IntegralSquare(Times[j-1],Times[j]);

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_);

Variates.resize(NumberOfTimes);

}

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)

#ifndef PATH_DEPENDENT_ASIAN_H

#define PATH_DEPENDENT_ASIAN_H

#include <PathDependent.h>

#include <PayOffBridge.h>

class PathDependentAsian : public PathDependent

{

public:

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;

private:

double DeliveryTime;

PayOffBridge ThePayOff;

unsigned long NumberOfTimes;

};

#endif

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

destruction.

The source ļ¬le is fairly simple too.

Listing 7.8 (PathDependentAsian.cpp)

#include <PathDependentAsian.h>

PathDependentAsian::PathDependentAsian(const MJArray&

LookAtTimes_,

double DeliveryTime_,

const PayOffBridge&ThePayOff_)

:

PathDependent(LookAtTimes_),

DeliveryTime(DeliveryTime_),

ThePayOff(ThePayOff_),

NumberOfTimes(LookAtTimes_.size())

{

}

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&

SpotValues,

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

AntiThetic.cpp

Arrays.cpp,

ConvergenceTable.cpp,

ExoticBSEngine.cpp

ExoticEngine.cpp

MCStatistics.cpp

Normals.cpp

Parameters.cpp,

ParkMiller.cpp,

PathDependent.cpp

PathDependentAsian.cpp

PayOff3.cpp,

PayOffBridge.cpp,

Random2.cpp,

*/

118 An exotics engine and the template pattern

#include<ParkMiller.h>

#include<iostream>

using namespace std;

#include<MCStatistics.h>

#include<ConvergenceTable.h>

#include<AntiThetic.h>

#include<PathDependentAsian.h>

#include<ExoticBSEngine.h>

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 =

gathererTwo.GetResultsSoFar();

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

Trees

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

1

. (8.3)

t

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

T

W(l+1)T /N ā’ WlT /N = N (0, 1). (8.4)

N

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

121

122 Trees

therefore take a set of N independent random variables X j with this distribution,

and approximate WlT /N by

l

T

Yl = X j. (8.5)

N j=1

The approximation for SlT /N is then

S0 e(r ā’dā’ 2 Ļ

2 )lT /N +Ļ Y

1

.

l

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)

2

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

1

,

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

up:

ā¢ 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

pay-off.

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

TreeProducts.h.

Listing 8.1 (TreeProducts.h)

#ifndef TREE_PRODUCTS_H

#define TREE_PRODUCTS_H

class TreeProduct

{

public:

TreeProduct(double FinalTime_);

virtual double FinalPayOff(double Spot) const=0;

virtual double PreFinalValue(double Spot,

double Time,

double DiscountedFutureValue)

const=0;

virtual ˜TreeProduct(){}

virtual TreeProduct* clone() const=0;

double GetFinalTime() const;

private:

double FinalTime;

};

#endif

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)

#ifndef TREE_AMERICAN_H

#define TREE_AMERICAN_H

#include <TreeProducts.h>

#include <PayOffBridge.h>

class TreeAmerican : public TreeProduct

{

public:

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 |