double DiscountedFutureValue)

const;

virtual ˜TreeAmerican(){}

private:

PayOffBridge ThePayOff;

};

#endif

The header for the European option is very similar:

Listing 8.4 (TreeEuropean.h)

#ifndef TREE_EUROPEAN_H

#define TREE_EUROPEAN_H

#include <TreeProducts.h>

#include <PayOffBridge.h>

class TreeEuropean : public TreeProduct

{

public:

TreeEuropean(double FinalTime,

const PayOffBridge& ThePayOff_);

virtual TreeProduct* clone() const;

virtual double FinalPayOff(double Spot) const;

virtual double PreFinalValue(double Spot,

double Time,

double DiscountedFutureValue)

const;

virtual ˜TreeEuropean(){}

private:

PayOffBridge ThePayOff;

};

#endif

128 Trees

The source ¬les are also straightforward.

Listing 8.5 (TreeAmerican.cpp)

#include <TreeAmerican.h>

#include <minmax.h>

TreeAmerican::TreeAmerican(double FinalTime,

const PayOffBridge& ThePayOff_)

: TreeProduct(FinalTime),

ThePayOff(ThePayOff_)

{

}

TreeProduct* TreeAmerican::clone() const

{

return new TreeAmerican(*this);

}

double TreeAmerican::FinalPayOff(double Spot) const

{

return ThePayOff(Spot);

}

double TreeAmerican::PreFinalValue(double Spot,

double ,

// Borland compiler doesnt like unused named variables

double DiscountedFutureValue) const

{

return max(ThePayOff(Spot), DiscountedFutureValue);

}

and

Listing 8.6 (TreeEuropean.cpp)

#include <TreeEuropean.h>

#include <minmax.h>

TreeEuropean::TreeEuropean(double FinalTime,

const PayOffBridge& ThePayOff_)

: TreeProduct(FinalTime),

ThePayOff(ThePayOff_)

8.4 A tree class 129

{

}

double TreeEuropean::FinalPayOff(double Spot) const

{

return ThePayOff(Spot);

}

double TreeEuropean::PreFinalValue(double,

//Spot, Borland compiler

double,

//Time, doesnt like unused named variables

double DiscountedFutureValue) const

{

return DiscountedFutureValue;

}

TreeProduct* TreeEuropean::clone() const

{

return new TreeEuropean(*this);

}

The implementations of the methods are very simple “ all the pay-offs are sub-

contracted to the PayOffBridged class in any case. For the European option at

an interior node, the rule for computing the PreFinalValue is very simple: just

return the discounted future value that was input. For the American option, it is

only slightly harder; we take the maximum with the intrinsic value.

Note the slight oddity that we do not name the unused variables Spot and Time

in the PreFinalValue method. This is because some compilers issue a warning

if a variable is named but not used, to ensure that variables are not left unused

accidentally.

8.4 A tree class

We give a simple implementation of a binomial tree class in this section. We design

the tree to work in three pieces. The constructor does little except store the parame-

ters. The BuildTree method actually makes the tree. In particular, it computes the

locations of all the nodes, and the discounts needed to compute the expectations

backwards in the tree. As it is not intended that the BuildTree method should be

called from outside the class, we make it protected which allows the possibil-

ity that an inherited class may wish to use it without allowing any other external

access.

130 Trees

The method which does the actual pricing is GetThePrice. Note that it takes in

a TreeProduct by reference. As the argument is an abstract base class this means

that an object from an inherited class must be passed in. Note that since we not

need to store the object passed in, we do not need virtual constructors, wrappers or

bridges. This method builds the tree if necessary, checks that the product has the

right expiry time and then prices it. Our design is such that we can price multiple

products with the same expiry; we call the method multiple times and only build

the tree once. As we wish to be able to do this, we store the entire tree. Note

that this is not really necessary since for any given time slice, we only need the

next time slice and so we could easily save a lot of memory by only ever having

two arrays de¬ned. However, unless one is doing an awful lot of steps, memory

will not be an issue, and this approach has the added bene¬t that if one wishes to

analyze certain aspects of the product, such as where the exercise boundary lies

for an American option, it is better to have the entire tree. We present the class in

BinomialTree.h.

Listing 8.7 (BinomialTree.h)

#pragma warning( disable : 4786 )

#include <TreeProducts.h>

#include <vector>

#include <Parameters.h>

#include <Arrays.h>

class SimpleBinomialTree

{

public:

SimpleBinomialTree(double Spot_,

const Parameters& r_,

const Parameters& d_,

double Volatility_,

unsigned long Steps,

double Time);

double GetThePrice(const TreeProduct& TheProduct);

protected:

void BuildTree();

private:

8.4 A tree class 131

double Spot;

Parameters r;

Parameters d;

double Volatility;

unsigned long Steps;

double Time;

bool TreeBuilt;

std::vector<std::vector<std::pair<double, double> > >

TheTree;

MJArray Discounts;

};

Note that we store the tree as a vector of vectors of pairs of doubles. This is why

we have the #pragma at the start; without the pragma, we get a warning message

telling us that the debug info is too long (under Visual C++ 6.0).

A pair is a simple template class in the STL which simply gives a class with

two data members of the appropriate types. They are accessed as public members

as first and second. Note that an alternative implementation would be to have

two trees: one for the spot and another for option values. However, that would

require twice as much work when resizing, and more importantly one would then

have to be careful to ensure that spot and the associated option value were always

attached to the same indices. With the use of a pair, we get this for free.

We have allowed general Parameters for r and d: this is because variable in-

terest and dividend rates change little in the analysis or construction of the tree. If

we can add extra functionality with little cost, why not do so? We do not, however,

allow variable volatility as it greatly complicates node placement; we would lose

the property that the spot price is a simple function of the underlying Brownian

motion.

We use a bool to indicate whether the tree has been built yet. We store the

number of steps as an unsigned long and this is ¬xed in the constructor. One

could easily add an additional method to allow the number of steps to be changed;

however one would gain little over just instantiating a new object with a different

number of steps. We have a data member for the discount factors needed so that

they need only be computed once in BuildTree.

We present the source code in BinomialTree.cpp.

Listing 8.8 (BinomialTree.cpp)

#include <BinomialTree.h>

#include <Arrays.h>

#include <cmath>

132 Trees

// the basic math functions should be in namespace std

// but aren™t in VCPP6

#if !defined(_MSC_VER)

using namespace std;

#endif

SimpleBinomialTree::SimpleBinomialTree(double Spot_,

const Parameters& r_,

const Parameters& d_,

double Volatility_,

unsigned long Steps_,

double Time_)

: Spot(Spot_),

r(r_),

d(d_),

Volatility(Volatility_),

Steps(Steps_),

Time(Time_),

Discounts(Steps)

{

TreeBuilt=false;

}

void SimpleBinomialTree::BuildTree()

{

TreeBuilt = true;

TheTree.resize(Steps+1);

double InitialLogSpot = log(Spot);

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

{

TheTree[i].resize(i+1);

double thisTime = (i*Time)/Steps;

double movedLogSpot =

InitialLogSpot + r.Integral(0.0, thisTime)

- d.Integral(0.0, thisTime);

8.4 A tree class 133

movedLogSpot -=

0.5*Volatility*Volatility*thisTime;

double sd = Volatility*sqrt(Time/Steps);

for (long j = -static_cast<long>(i), k=0;

j <= static_cast<long>(i); j=j+2,k++)

TheTree[i][k].first = exp(movedLogSpot+ j*sd);

}

for (unsigned long l=0; l <Steps; l++)

{

Discounts[l] =

exp(- r.Integral(l*Time/Steps,(l+1)*Time/Steps));

}

}

double SimpleBinomialTree::GetThePrice(const TreeProduct&

TheProduct)

{

if (!TreeBuilt)

BuildTree();

if (TheProduct.GetFinalTime() != Time)

throw("mismatched product in SimpleBinomialTree");

for (long j = -static_cast<long>(Steps), k=0;

j <=static_cast<long>( Steps); j=j+2,k++)

TheTree[Steps][k].second =

TheProduct.FinalPayOff(TheTree[Steps][k].first);

for (unsigned long i=1; i <= Steps; i++)

{

unsigned long index = Steps-i;

double ThisTime = index*Time/Steps;

for (long j = -static_cast<long>(index), k=0;

j <= static_cast<long>(index); j=j+2,k++)

{

double Spot = TheTree[index][k].first;

134 Trees

double futureDiscountedValue = 0.5*Discounts[index]*

(TheTree[index+1][k].second +

TheTree[index+1][k+1].second);

TheTree[index][k].second =

TheProduct.PreFinalValue(Spot,ThisTime,

futureDiscountedValue);

}

}

return TheTree[0][0].second;

}

The code is fairly straightforward. The constructor initializes the basic class vari-

ables and does not do much else. The method BuildTree creates the tree. We

resize the vector describing all the layers ¬rst. We then resize each layer so that

it is of the correct size. Note that as the number of nodes in a layer grows with the

number of steps, these inner vectors are all of different sizes. We then compute a

basepoint for the layer in log-space which is just the zero Brownian motion point.

We then loop through the nodes in each layer writing in the appropriate value of

spot.

Note that as we are dealing with points in the tree corresponding to down moves

we count using a long. This long has to be compared to the unsigned long i

for the termination condition. We therefore have to be careful; if we simply com-

pare these numbers a likely effect is that the routine will conclude that ’1 is bigger

than 1. Why? The compiler will implicitly convert the long into a large unsigned

long. Clearly this is not the desired effect so we convert the unsigned long into

a long before the comparison.

After setting up the tree, we also set up the array of Discounts which are

product-independent, and therefore only need to be done once.

The main routine for actually doing the pricing, GetThePrice, is also straight-

forward. We have put in a test to make sure that the tree and product are compatible.

This throws an error which is a simple string if they are not. Note that this means

the program will terminate unless we have written a catch which is capable of

catching the string.

We simply iterate backwards through the tree. First we compute the ¬nal layer

using the FinalPayOff and write the values into the second element of each pair

in the ¬nal layer. After this we simply iterate backwards, computing as we go. The

¬nal value of the product is then simply the value of the option at the single node

in the ¬rst layer of the tree, and that is what we return.

8.5 Pricing on the tree 135

8.5 Pricing on the tree

Having written all the classes we need to actually put them together to price some-

thing. We give a simple interface in TreeMain.cpp.

Listing 8.9 (TreeMain.cpp)

/*

requires

Arrays.cpp

BinomialTree.cpp

BlackScholesFormulas.cpp

Normals.cpp

Parameters.cpp

PayOff3.cpp

PayOffBridge.cpp

PayOffForward.cpp

TreeAmerican.cpp

TreeEuropean.cpp

TreeProducts.cpp

*/

#include <BinomialTree.h>

#include <TreeAmerican.h>

#include <TreeEuropean.h>

#include <BlackScholesFormulas.h>

#include <PayOffForward.h>

#include <iostream>

using namespace std;

#include <cmath>

int main()

{

double Expiry;

double Strike;

double Spot;

double Vol;

double r;

double d;

unsigned long Steps;

136 Trees

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 << "\nNumber of steps\n";

cin >> Steps;

PayOffCall thePayOff(Strike);

ParametersConstant rParam(r);

ParametersConstant dParam(d);

TreeEuropean europeanOption(Expiry,thePayOff);

TreeAmerican americanOption(Expiry,thePayOff);

SimpleBinomialTree theTree(Spot,rParam,dParam,Vol,Steps,

Expiry);

double euroPrice = theTree.GetThePrice(europeanOption);

double americanPrice = theTree.GetThePrice(americanOption);

cout << "euro price" << euroPrice << "amer price"

<< americanPrice << "\n";

double BSPrice = BlackScholesCall(Spot,Strike,r,d,

Vol,Expiry);

cout << "BS formula euro price" << BSPrice << "\n";

8.5 Pricing on the tree 137

PayOffForward forwardPayOff(Strike);

TreeEuropean forward(Expiry,forwardPayOff);

double forwardPrice = theTree.GetThePrice(forward);

cout << "forward price by tree" << forwardPrice << "\n";

double actualForwardPrice =

exp(-r*Expiry)*(Spot*exp((r-d)*Expiry)-Strike);

cout << "forward price" << actualForwardPrice << "\n";

Steps++; // now redo the trees with one more step

SimpleBinomialTree theNewTree(Spot,rParam,dParam,Vol,

Steps,Expiry);

double euroNewPrice =

theNewTree.GetThePrice(europeanOption);

double americanNewPrice =

theNewTree.GetThePrice(americanOption);

cout << "euro new price" << euroNewPrice

<< "amer new price" << americanNewPrice << "\n";

double forwardNewPrice = theNewTree.GetThePrice(forward);

cout << "forward price by new tree" << forwardNewPrice

<< "\n";

double averageEuro = 0.5*(euroPrice+euroNewPrice);

double averageAmer = 0.5*(americanPrice+americanNewPrice);

double averageForward = 0.5*(forwardPrice+forwardNewPrice);

cout << "euro av price" << averageEuro << "amer av price"

<< averageAmer << "\n";

cout << "av forward" << averageForward << "\n";

double tmp;

cin >> tmp;

return 0;

}

138 Trees

To illustrate certain aspects of tree pricing, we price the European call, the Amer-

ican call and the forward. We then reprice them using one extra step. The rea-

son we do this is that often pricing on trees gives rise to a zig-zag behaviour as

nodes go above and below the money. The average of the price for two succes-

sive steps is therefore often a lot more accurate than a single price. We give the

comparison Black“Scholes price for the European option as an assessment of ac-

curacy. The code for the Black“Scholes price is in BlackScholesFormulas.h

and BlackScholesFormulas.cpp; we discuss it in Appendix A.

A standard way of improving the accuracy of American option pricing is to use

the European price as a control. As we know the true price of the European and

we know the tree price also, we can assume that the American option has the same

amount of error as the European and adjust its price accordingly. The principle here

is the same as that for a control variate in Monte Carlo simulation.

We also give the price of a forward. Note that if we had required the discounted

discretized process to be a martingale, then the forward would be priced absolutely

correctly; however as it is only an approximation to a martingale rather than actu-

ally being one, the forward need not be precise.

As we have not de¬ned a class for the forward™s pay-off before, we de¬ne one

in PayOffForward.h and PayOffForward.cpp.

Listing 8.10 (PayOffForward.h)

#ifndef PAY_OFF_FORWARD_H

#define PAY_OFF_FORWARD_H

#include <PayOff3.h>

class PayOffForward : public PayOff

{

public:

PayOffForward(double Strike_);

virtual double operator()(double Spot) const;

virtual ˜PayOffForward(){}

virtual PayOff* clone() const;

private:

double Strike;

};

#endif

8.7 Exercises 139

Listing 8.11 (PayOffForward.cpp)

#include <PayOffForward.h>

double PayOffForward::operator () (double Spot) const

{

return Spot-Strike;

}

PayOffForward::PayOffForward(double Strike_) : Strike(Strike_)

{

}

PayOff* PayOffForward::clone() const

{

return new PayOffForward(*this);

}

The class is straightforward and the only difference from the class de¬ned for

the call is that we take spot minus strike, instead of the call pay-off.

8.6 Key points

In this chapter we have used the patterns developed earlier in the book to develop

routines for pricing on trees.

• Tree pricing is based on the discretization of a Brownian motion.

• Trees are a natural way to price American options.

• On a tree, knowledge of discounted future values is natural but knowing about

the past is not.

• We can re-use the pay-off class when de¬ning products on trees.

• By having a separate class encapsulating the de¬nition of a derivative on a tree,

we can re-use the products for more general structures.

• European options can be used as controls for American options.

8.7 Exercises

We have developed a very simple tree and treated a couple of simple products.

The approach here can easily be extended to many more cases. We suggest a few

possibilities for the reader to try.

140 Trees

Exercise 8.1 Find a class that does barrier options in the same TreeProduct class

hierarchy. Try it out. How stable is the price? How might you improve the stability?

Exercise 8.2 Develop a binomial tree for which the memory requirements grow

linearly with the number of steps. How do the memory requirements grow for the

class here?

Exercise 8.3 Write a trinomial tree class.

Exercise 8.4 Modify the code so that it will work under variable volatility. The key

is to ensure that the integral of the square of the vol across each time step is the

same. This means that the time steps will be of unequal length.

Exercise 8.5 Modify the tree so the implied stock price process makes the dis-

counted price a martingale. Compare convergence for calls, puts and forwards.

Exercise 8.6 Implement an American knock-in option pricer on a tree. (Use an

additional auxiliary variable to indicate whether or not the option has knocked-in,

and compute the value at each node in both cases.)

9

Solvers, templates, and implied volatilities

9.1 The problem

Whatever model one is using to compute vanilla options™ prices, it is traditional to

quote prices in terms of the Black“Scholes implied volatility. The implied volatility

is by de¬nition the number to plug into the Black“Scholes formula to get the price

desired. Thus we have the problem that we must solve the problem of ¬nding the

value σ such that

BS(S, K , r, d, T, σ ) = quoted price.

In other words, we must invert the map

σ ’ BS(S, K , r, d, T, σ )

with the other parameters ¬xed.

The Black“Scholes formula is suf¬ciently complicated that there is no analytic

inverse and this inversion must be carried out numerically. There are many algo-

rithms for implementing such inversions; we shall study two of the simplest: bisec-

tion and Newton“Raphson. Our objective, as usual, is to illustrate the programming

techniques for de¬ning the interfaces in a reusable fashion rather than to implement

the most ef¬cient algorithms available. Indeed, we hope that the reader will com-

bine the techniques here with algorithms found elsewhere, in for example [28], to

produce robust and ef¬cient reusable code.

Before proceeding to the coding and design issues, we recall the details of the

aforementioned algorithms. Given a function, f , of one variable we wish to solve

the equation

f (x) = y. (9.1)

In the above, f is the Black“Scholes formula, x is volatility and y is the price. If

141

142 Solvers, templates, and implied volatilities

the function f is continuous, and for some a and b we have

f (a) < y, (9.2)

f (b) > y, (9.3)

then there must exist some c in the interval (a, b) such that f (c) = x. Bisection is

one technique to ¬nd c. The idea is straightforward: we simply take the midpoint,

m, of the interval, then one of three things must occur:

• f (m) = y and we are done;

• f (m) < y in which case there must be a solution in (m, b);

• f (m) > y in which case there must be a solution in (a, m).

Thus by taking the midpoint, we either ¬nd the solution, or halve the size of the

interval in which the solution exists. Repeating we must narrow in on the solution.

In practice, we would terminate when we achieve

| f (m) ’ y| < , (9.4)

for some pre-decided tolerance, .

Bisection is robust but is not particularly fast. When we have a well-behaved

function with an analytic derivative then Newton“Raphson can be much faster. The

idea of Newton“Raphson is that we pretend the function is linear and look for the

solution where the linear function predicts it to be. Thus we take a starting point,

x0 , and approximate f by

g0 (x) = f (x0 ) + (x ’ x0 ) f (x0 ). (9.5)

We have that g0 (x) is equal to zero if and only if

y ’ f (x0 )

x= + x0 . (9.6)

f (x0 )

We therefore take this value as our new guess x1 . We now repeat until we ¬nd that

f (xn ) is within of y.

Newton“Raphson is much faster than bisection provided we have an easily eval-

uated derivative. This is certainly the case for the Black“Scholes function. Indeed,

for a call option the vega is easier to compute than the price is. However, as it in-

volves passing two functions rather than one to a solver routine, it requires more

sophisticated programming techniques to implement re-usably.

9.2 Function objects

We want to implement the bisection algorithm in a re-usable way; this means that

we will need a way to pass the function f into the bisection routine. Since f may

well be de¬ned, as it is in our case, in terms of the value of a more complicated

9.2 Function objects 143

function with many parameters ¬xed, we will also need somehow to pass in the

values of those auxiliary parameters.

There are, in fact, many different ways to tackle this problem. One method we

have already studied is the engine template. With this approach, we de¬ne a base

class for which the main method is to carry out the bisection. The main method

calls a pure virtual method to get the value of f (x). For any speci¬c problem,

we then de¬ne an inherited class which implements f appropriately. Whilst this

method can work effectively, there are a couple of disadvantages. The ¬rst is that

the function call is virtual which can lead to ef¬ciency problems. There are two

causes: the ¬rst is that to call a virtual function, the processor has to look up a

virtual function table each time the function is called, and then jump to a location

speci¬ed by the table. Clearly, it would be slightly faster not to have to look up

the table. A more subtle and serious speed issue is that it is not possible to inline

virtual functions. If the function is known beforehand, the compiler can inline it and

eliminate the mechanics of the function call altogether. In addition, the compiler

may be able to make additional optimizations as it sees all the code together at

once. Whilst these speed issues are not particularly important whilst designing a

solver, they are more critical when writing a class for other numerical routines such

as numeric integration where often a large of calls are made to a function which is

fast to evaluate. We therefore wish to develop a pattern which can be used in those

contexts too.

The second disadvantage of inheriting from a solver base class is that it inhibits

other inheritance. If we wish to inherit the class de¬ning our function from some

other class, we cannot inherit from the solver class as well without using multiple

inheritance. Of course, one could use multiple inheritance but it tends to be tricky

to get it to work in a bug-free fashion and I therefore tend to avoid it.

Having decided that we want to be able to input a function to our routine with-

out using virtual functions, what other options do we have? One solution would

be to use a function pointer but this would buy us little (if anything) over virtual

functions. Another approach is templatization. The crucial point is that with tem-

platization the type of the function being used in the optimization is decided at

compile time rather than at runtime. This means that the compiler can carry out

optimizations and inlining that depend on the type of the function since that infor-

mation is now available to it.

The approach we adopt for specifying the function we wish to optimize uses

the function object. We ¬rst encountered function objects in Section 2.1 when

de¬ning pay-offs. Recall that a function object is by de¬nition an object for which

operator() is de¬ned. So if we have an object f of a class T for which

const operator()( double x) const

144 Solvers, templates, and implied volatilities

has been de¬ned it is then legitimate to write f(y) for a double y, and this is

equivalent to

f.operator()(y).

Thus our object f can be used with function-like syntax. However, as f is an object

it can contain extra information. Thus if we want to solve for the implied volatility,

the function object will take the volatility as an argument, but will also have, as

extra parameters already stored, the values of r, d, T, S and K .

We thus obtain the class de¬ned in BSCallClass.h:

Listing 9.1

//

// BSCallClass.h

//

#ifndef BS_CALL_CLASS_H

#define BS_CALL_CLASS_H

class BSCall

{

public:

BSCall(double r_, double d_,

double T, double Spot_,

double Strike_);

double operator()(double Vol) const;

private:

double r;

double d;

double T;

double Spot;

double Strike;

};

#endif

The source ¬le is simple:

9.3 Bisecting with a template 145

Listing 9.2

//

// BSCallClass.cpp

//

#include <BSCallClass.h>

#include <BlackScholesFormulas.h>

BSCall::BSCall(double r_, double d_,

double T_, double Spot_,

double Strike_)

:

r(r_),d(d_),

T(T_),Spot(Spot_),

Strike(Strike_)

{}

double BSCall::operator()(double Vol) const

{

return BlackScholesCall(Spot,Strike,r,d,Vol,T);

}

The constructor simply initializes the class data members, which are the parame-

ters needed to price a call option under Black“Scholes except the volatility. The

operator() takes in the volatility and then invokes the Black“Scholes formula.

This is the simplest possible implementation of the class. If we were truly wor-

ried about ef¬ciency considerations, we could code the formula directly and pre-

compute as much of it as possible in the constructor. We could have for example a

class data member, Moneyness, set to the log of Spot divided by Strike, and then

we would not have to compute it every time.

9.3 Bisecting with a template

In the previous section, we showed how we could de¬ne a class for which the syn-

tax f(x) makes sense when f was an object of the class, and x was a double. We

still need to get the object f into our solver routine, however. We do so via templa-

tization. The basic idea of templatization is that you can write code that works for

many classes simultaneously provided they are required to have certain operations

de¬ned with the same syntax. In this case, our requirement is that the class should

have

double operator()( double ) const

146 Solvers, templates, and implied volatilities

de¬ned, and thus that the syntax f(y) is well-de¬ned for class objects as we dis-

cussed above.

We present the Bisection function in Bisection.h:

Listing 9.3 (Bisection.h)

template<class T>

double Bisection(double Target,

double Low,

double High,

double Tolerance,

T TheFunction)

{

double x=0.5*(Low+High);

double y=TheFunction(x);

do

{

if (y < Target)

Low = x;

if (y > Target)

High = x;

x = 0.5*(Low+High);

y = TheFunction(x);

}

while

( (fabs(y-Target) > Tolerance) );

return x;

}

We only present a header ¬le, since for template code we cannot precompile in a

source ¬le “ we do not know the type of the object T. The function is quite simple.

We specify that it is templatized via the template<class T> at the top. If we

invoke the function with the template argument BSCall via Bisection<BSCall>

then every T will be converted into a BSCall before the function is compiled.

Once we have ¬xed the type of the template argument, there is really very little to

the function. We take the midpoint of the interval evaluate the function there, and

9.3 Bisecting with a template 147

switch to the left or right side of the interval by rede¬ning Low and High until the

value at the midpoint is close enough to the target, and we then return.

Note that we have de¬ned the type of the function object passed in as T The-

Function: we could equally well have put const T& TheFunction. The syntax

we have adopted involves copying the function object, and is therefore arguably

less good than the alternative. The reason I have done it this way is to highlight

the fact that the standard template library always uses the former syntax. A conse-

quence of this is that one needs to be careful when using function objects with the

STL not to de¬ne function objects which are expensive to copy (or, even worse,

impossible to copy.)

We now give a simple example of an implied volatility function:

Listing 9.4 (SolveMain1.cpp)

/*

Needs

BlackScholesFormulas.cpp

BSCallClass.cpp

Normals.cpp

*/

#include <Bisection.h>

#include <cmath>

#include <iostream>

#include <BSCallClass.h>

#include <BlackScholesFormulas.h>

using namespace std;

int main()

{

double Expiry;

double Strike;

double Spot;

double r;

double d;

double Price;

cout << "\nEnter expiry\n";

cin >> Expiry;

cout << "\nStrike\n";

148 Solvers, templates, and implied volatilities

cin >> Strike;

cout << "\nEnter spot\n";

cin >> Spot;

cout << "\nEnter price\n";

cin >> Price;

cout << "\nr\n";

cin >> r;

cout << "\nd\n";

cin >> d;

double low,high;

cout << "\nlower guess\n";

cin >> low;

cout << "\nhigh guess\n";

cin >> high;

double tolerance;

cout << "\nTolerance\n";

cin >> tolerance;

BSCall theCall(r,d,Expiry,Spot,Strike);

double vol =

Bisection(Price,low,high,tolerance,theCall);

double PriceTwo =

BlackScholesCall(Spot,Strike,r,d,vol,Expiry);

cout << "\n vol " << vol << " pricetwo "

<< PriceTwo << "\n";

double tmp;

cin >> tmp;

return 0;

}

9.4 Newton“Raphson and function template arguments 149

As usual, we input all the necessary parameters. We then put them together to

create a BSCall object. We then call the Bisection function to ¬nd the volatility.

Note that we have not put Bisection<BSCall>. The compiler deduces the type of

the template argument from the fact that our ¬nal argument in the function call is

theCall. We only have to specify the template argument when the compiler does

not have a suf¬cient amount of other information to deduce it.

Our function ¬nishes by getting the price via the Black“Scholes functions for

the implied volatility that was found. If everything is working correctly then this

will give the original price inputted.

9.4 Newton“Raphson and function template arguments

We now want to adapt the pattern we have presented to work for Newton“Raphson

as well as for bisection. The fundamental difference from a design viewpoint is that

Newton“Raphson involves two functions, the value and the derivative, whereas bi-

section involves just one. One solution would be simply to pass in two function

objects: one for the value and another for the derivative. This is unappealing, how-

ever, in that we would then need to initialize a set of parameters for each object and

we would have to be careful to make sure they are the same. More fundamentally,

the value and the derivative are really two aspects of the same thing rather than two

separate functions and so having two objects does not express well our conceptual

model of them.

A second solution is to assume a name for the derivative function. After all, that

is essentially what we did for the value function; it was a special name with special

syntax but ultimately it was just assuming a name. Thus we could assume that the

class had a method

double Derivative(double ) const

de¬ned and at appropriate points in our function we would then put

TheFunction.Derivative(x).

This would certainly work. However, it is a little ugly and if our class already had

a derivative de¬ned under a different name, it would be annoying.

Fortunately, there is a way of specifying which class member function to call

at compile time using templatization. The key to this is a pointer to a member

function. A pointer to a member function is similar in syntax and idea to a function

pointer, but it is restricted to methods of a single class. The difference in syntax is

that the class name with a :: must be attached to the * when it is declared. Thus

to declare a function pointer called Derivative which must point to a method of

150 Solvers, templates, and implied volatilities

the class T, we have

double (T::*Derivative)(double) const

The function Derivative is a const member function which takes in a double

as argument and outputs a double as return value. If we have an object of class T

called TheObject and y is a double, then the function pointed to can be invoked

by

TheObject.*Derivative(y)

Whilst the syntax is a little cumbersome, the key is to realise that it is the same for

ordinary function pointers except that T:: must be added to the * for declarations

and TheObject. must be added for invocations.

We can now use a function pointer to specify both the derivative and the value.

As we would like to avoid the time spent on evaluating the pointers, we can make

them template parameters rather than arguments to our function. This means that

the compiler can treat them just like any other function call, and as their types are

decided at compile time they can be inlined.

Our Newton“Raphson routine is therefore as follows:

Listing 9.5 (NewtonRaphson.h)

template<class T, double (T::*Value)(double) const,

double (T::*Derivative)(double) const >

double NewtonRaphson(double Target,

double Start,

double Tolerance,

const T& TheObject)

{

double y = (TheObject.*Value)(Start);

double x=Start;

while ( fabs(y - Target) > Tolerance )

{

double d = (TheObject.*Derivative)(x);

x+= (Target-y)/d;

y = (TheObject.*Value)(x);

}

return x;

}

9.5 Using Newton“Raphson to do implied volatilities 151

We have three template parameters: the class, the pointer to the value function for

that class, and the pointer to the derivative function for that class. The routine is

short and simple, now that we have the right structure. As usual, we keep repeating

until close enough to the root. We have not included the checks required to ensure

that the loop does not repeat endlessly if the sequence fails to converge to a root,

but such additions are easily put in.

9.5 Using Newton“Raphson to do implied volatilities

Now that we have developed a Newton“Raphson routine, we want to use it to

compute implied volatilities. Our class will therefore have to support pricing as a

function of volatility and the vega as a function of volatility. As before, the other

parameters will be class data members which are not inputted in the constructor

rather than via these methods. We present a suitable class in BSCallTwo.h and

BSCallTwo.cpp.

Listing 9.6 (BSCallTwo.h)

#ifndef BS_CALL_TWO_H

#define BS_CALL_TWO_H

class BSCallTwo

{

public:

BSCallTwo(double r_, double d_,

double T, double Spot_,

double Strike_);

double Price(double Vol) const;

double Vega(double Vol) const;

private:

double r;

double d;

double T;

double Spot;

double Strike;

};

#endif

152 Solvers, templates, and implied volatilities

The methods just call the relevant functions. As before, we could optimize by pre-

computing as much as possible, and inlining the methods.

Listing 9.7 (BSCallTwo.cpp)

#include <BSCallTwo.h>

#include <BlackScholesFormulas.h>

BSCallTwo::BSCallTwo(double r_, double d_,

double T_, double Spot_,

double Strike_)

:

r(r_),d(d_),

T(T_),Spot(Spot_),

Strike(Strike_)

{}

double BSCallTwo::Price(double Vol) const

{

return BlackScholesCall(Spot,Strike,r,d,Vol,T);

}

double BSCallTwo::Vega(double Vol) const

{

return BlackScholesCallVega(Spot,Strike,r,d,Vol,T);

}

We present an example of using the combination of NewtonRaphson and

BSCallTwo in SolveMain2.cpp.

Listing 9.8 (SolveMain2.cpp)

/*

Needs

BlackScholesFormulas.cpp

BSCallTwo.cpp

Normals.cpp

*/

#include <NewtonRaphson.h>

#include <cmath>

#include <iostream>

#include <BSCallTwo.h>

#include <BlackScholesFormulas.h>

9.5 Using Newton“Raphson to do implied volatilities 153

using namespace std;

int main()

{

double Expiry;

double Strike;

double Spot;

double r;

double d;

double Price;

cout << "\nEnter expiry\n";

cin >> Expiry;

cout << "\nStrike\n";

cin >> Strike;

cout << "\nEnter spot\n";

cin >> Spot;

cout << "\nEnter price\n";

cin >> Price;

cout << "\nr\n";

cin >> r;

cout << "\nd\n";

cin >> d;

double start;

cout << "\nstart guess\n";

cin >> start;

double tolerance;

cout << "\nTolerance\n";

cin >> tolerance;

BSCallTwo theCall(r,d,Expiry,Spot,Strike);

154 Solvers, templates, and implied volatilities

double vol=NewtonRaphson<BSCallTwo, &BSCallTwo::Price,

&BSCallTwo::Vega>(Price, start,

tolerance, theCall);

double PriceTwo =

BlackScholesCall(Spot,Strike,r,d,vol,Expiry);

cout << "\n vol " << vol << " \nprice two:" << PriceTwo

<< "\n";

double tmp;

cin >> tmp;

return 0;

}

Our new main program is very similar to the one we had before. The main change

is that this time we specify the template parameters for our NewtonRaphson func-

tion, whereas for the Bisection function we did not bother. The reason for the

change is that there is not suf¬cient information for the compiler to deduce the

types of the parameters. There is nothing to indicate which member functions of

the class are to be used. Even for our class which only has two member functions,

these two functions could equally well be the other way round as far the compiler

knows.

9.6 The pros and cons of templatization

In this chapter, we have used template arguments to achieve re-usability whereas

in other chapters, we have used virtual functions and polymorphism. There are

advantages and disadvantages to each approach. The principal difference is that

for templates argument types are decided at the time of compilation, whereas for

virtual functions the type is not determined until runtime.

What consequences does this difference have? The ¬rst is speed. No time is

spent on deciding which code to run when the code is actually running. In addition,

the fact that the compiler knows which code will be run allows it to make extra

optimizations which would be hard if not impossible when the decision is made at

run time.

A second consequence is size. As the code is compiled for each template argu-

ment used separately, we have multiple copies of very similar code. For a simple

9.6 The pros and cons of templatization 155

routine such as a solver this is not really an issue but for a complicated routine, this

could result in a very large executable. Another aspect of this is slower compiler

times since much more code would have to be compiled. If we had several tem-

plate parameters the size could multiply out of control. For example, suppose when

designing our Monte Carlo path-dependent exotic pricer, we had templatized both

the random number generator and the product. If we had six random number gen-

erators and ten products, and we wished to allow any combination then we would

have sixty times as much code.

A third consequence is that it becomes harder for the user of the code to make

choices. In the example of the exotics pricer, if the user was allowed to choose the

number generator and the product via outside input, we would have to write code

that branched into each of the sixty cases and within each branch called the engine

and gathered results.

As well as the run time versus compile time decision, there are other issues with

using templatized code. A simple one is that it is harder to debug. Some debuggers

get confused by template code and, for example, refuse to set breakpoints within

templates (or else set and ignore them.) Related to this is the fact that compilers

will often not actually compile lines of template code that are not used. Thus if a

templatized class has a method that is not called anywhere, the code will compile

even if it has syntax errors. Only when a line is added that calls the particular

method will the compiler errors appear. This can be infuriating as the error may

show up a long time afterwards.

One way of avoiding these problems is ¬rst to write non-template code for a

particular choice of the template parameter. This code can be thoroughly tested and

debugged, and then afterwards the code can be rewritten by changing the particular

parameter into a template parameter.

So when should one use templates and when use virtual functions? My prefer-

ence is not to use templates unless certain conditions are met. These are that the

routine should be short, and potentially re-usable in totally unrelated contexts. So

for example, I would use templates for a numerical integration routine and a solver.

I would also use templates for a container class; in fact I would generally use the

templates given in the standard template library. I would not, however, use tem-

plates for an option pricing engine since the code would probably be long and is

only relevant to a quite speci¬c context.

The general trend in C++ is towards templates. The principal reason is that

they are the way to achieve the same speed as lower-level languages. In general,

languages exhibit a trade-off between abstraction and ef¬ciency; C++ has always

striven to achieve both. Templates are ultimately a way of achieving abstraction

without sacri¬cing ef¬ciency.

156 Solvers, templates, and implied volatilities

9.7 Key points

In this chapter we have looked at how to implement solvers using template code.

• Templates are an alternative to inheritance for coding without knowing an ob-

ject™s precise type.

• Template code can be faster as function calls are determined at compile time.

• Extensive use of template code can lead to very large executables.

• Pointers to member functions can be a useful way of obtaining generic behaviour.

• Implied volatility can only be computed numerically.

9.8 Exercises

Exercise 9.1 Modify the Newton“Raphson routine so that it does not endlessly

loop if a root is not found.

Exercise 9.2 Take your favourite numerical integration routine, e.g. the trapezium

rule, and write a template routine to carry it out.

Exercise 9.3 Write a routine to price a vanilla option by Monte Carlo or trees where

the pay-off is passed in as a template parameter expressed via a function object.

10

The factory

10.1 The problem

Suppose we wish to design an interface which is a little more sophisticated than

those we have used so far. The user will input the name of a pay-off and a strike,

and the program will then price a vanilla option with that pay-off. We therefore

need a conversion routine from strings and strikes to pay-offs. How might we write

this?

One simple solution is to write a function that takes in the string and the strike,

checks against all known types of pay-offs and when it comes across the right one,

creates a pay-off of the right type. We would probably implement this via a switch

statement. Our conversion routine would then have to include the header ¬les for

all possible forms of pay-off, and every time we added a new pay-off we would

have to modify the switch statement. Clearly, this solution violates the open-closed

principle as any addition involves modi¬cation.

In this chapter, we present a solution that allows us to add new pay-offs without

changing any of the existing ¬les. We simply add new ¬les to the project. Our

solution is a design pattern known as the factory pattern. It is so called because it

can be used to manufacture objects. Whilst we restrict attention to a simple factory

which manufactures pay-offs, the basic pattern can be used in much wider contexts.

10.2 The basic idea

Our solution requires each type of pay-off to tell the factory that it exists, and to

give the factory a blueprint for its manufacture. In this context a blueprint means

an identity string to distinguish that class and a pointer to a function that will create

objects of that class.

How can we get the class to communicate with the factory, without explicitly

calling anything from the main routine? The key lies in global variables. Every

global variable is initialized when the program commences before anything else

157

158 The factory

happens. If we de¬ne a class in such a way that initializing a global variable of

that class registers a pay-off class with the factory, then we have achieved what we

wanted. This is possible because the initialization involves a call to a constructor,

and we can make the constructor do whatever we want.

So for each pay-off class, we write an auxiliary class whose constructor registers

the pay-off class with our factory, and we declare a global variable of the auxiliary

class. In fact, as these auxiliary classes will all be very similar to each other, we

adopt a template solution for de¬ning these classes.

We also need a factory object for these auxiliary classes to talk to. We cannot

make this factory object a global variable, as we have no control over the order in

which the global variables are initialized. We need it to exist before the other glob-

als are de¬ned as they refer to it. Fortunately, there is a type of variable guaranteed

to come into existence at the moment it is ¬rst referred to: the static variable.

Thus if the registration function contains a static variable which is the factory,

then on the ¬rst call to the registration function the factory comes into existence.

Recall that a static variable de¬ned in a function persists from one call to the

next, and only disappears when the program exits. So all the registration function

calls will register the pay-off blueprints with the same factory.

However, we are not yet done because the creator function will need to have

access to the same factory object as the registration function, and if the factory is

hidden inside the registration function this will not be possible. The solution to this

problem is known as the singleton pattern.

10.3 The singleton pattern

We saw in the last section that we need to de¬ne a factory via a static variable

since it must come into existence as soon as it is referred to when registering the

blueprints. We also saw that the same factory must be referred to by every regis-

tration, and that the same factory will be needed when creating the pay-offs from

strings.

So what we need is a factory class for which an object exists as soon as it is

required, and for this object to exist until the end of the program. We also do not

want any other factory objects to exist as they will just confuse matters; everything

must be registered with and built by the same factory.

The singleton pattern gives a way of creating a class with these properties. The

¬rst thing is that all constructors and assignment operators are made private. This

means that factories can only be created from inside methods of the class this

gives us ¬rm control over the existence of factory objects. In order to get the one

class object that we need, we de¬ne a very simple method that de¬nes a class ob-

ject as a static variable. Thus if our class is called PayOffFactory, we de¬ne a

class method Instance as follows:

10.4 Coding the factory 159

PayOffFactory& PayOffFactory::Instance()

{

static PayOffFactory theFactory;

return theFactory;

}

The ¬rst time that Instance is called, it creates the static data member the-

Factory. As it is a member function, it can do this by using the private default

constructor. Every subsequent time the Instance is called, the address of the

already-existing static variable theFactory is returned. Thus Instance cre-

ates precisely one PayOffFactory object which can be accessed from anywhere

by calling PayOffFactory::Instance().

Note that Instance will have to be a static method of PayOffFactory, as

the whole point is that it provides you with a PayOffFactory object, and it would

be useless if you had to access it from an existing object. Note also that the mean-

ing of static here for a function is quite different from the meaning above for

a variable; for a function it means that the function can be called directly with-

out any attachment to an object. One still has to pre¬x with the name of class,

however.

We have now achieved what we needed: we have a way of creating a single

factory which can be referenced from anywhere at any time in the program. Note

that the name singleton pattern was chosen because precisely one object from the

class exists.

10.4 Coding the factory

In the last section, we saw how the singleton pattern could be used to create a single

factory accessible in any place at any time. We now use this pattern to implement

the factory. As well as the instance method discussed above, we will need a method

for registering pay-off classes and a method for creating then.

How will registration work? Upon registration, we need to know the string iden-

ti¬er for the speci¬c pay-off class and the pointer to the function which actually

creates the object in question. These will therefore be the arguments for the regis-

tration method. The factory will need to store this information for when the create

pay-off method is called. This will require a container class.

Fortunately, there is a container in the standard template library which is de-

signed for associating identi¬ers to objects. This container is called the map class.

We therefore need a data member which is a map with template arguments std::

string and pointers to create functions.

Finally, we need a method which turns a string plus a strike into a PayOff object.

Our header ¬le is therefore as follows:

160 The factory

Listing 10.1 (PayOffFactory.h)

#ifndef PAYOFF_FACTORY_H

#define PAYOFF_FACTORY_H

#include <PayOff3.h>

#if defined(_MSC_VER)

#pragma warning( disable : 4786)

#endif

#include <map>

#include <string>

class PayOffFactory

{

public:

typedef PayOff* (*CreatePayOffFunction)(double );

static PayOffFactory& Instance();

void RegisterPayOff(std::string, CreatePayOffFunction);

PayOff* CreatePayOff(std::string PayOffId,double Strike);

˜PayOffFactory(){};

private:

std::map<std::string, CreatePayOffFunction>

TheCreatorFunctions;

PayOffFactory(){}

PayOffFactory(const PayOffFactory&){}

PayOffFactory& operator=

(const PayOffFactory&){ return *this;}

};

#endif

Note the typedef: this allows us to refer to pointers to functions which take

in a double and spit out a PayOff* as CreatePayOffFunction. Without this

typedef, the syntax would quickly become unmanageable. Note also that we make

the CreatePayOff method return a pointer to a PayOff. The reason for this is that

it allows the possibility of returning a null pointer if the identity string was not

found; otherwise we would have to throw an error or return a default sort of pay-

off.

We present the source code in PayOffFactory.cpp:

10.4 Coding the factory 161

Listing 10.2 (PayOffFactory.cpp)

#if defined(_MSC_VER)

#pragma warning( disable : 4786)

#endif

#include <PayOffFactory.h>

#include <iostream>

using namespace std;

void PayOffFactory::RegisterPayOff(string PayOffId,

CreatePayOffFunction CreatorFunction)

{

TheCreatorFunctions.insert(pair<string,CreatePayOffFunction>

(PayOffId,CreatorFunction));

}

PayOff* PayOffFactory::CreatePayOff(string PayOffId,

double Strike)

{

map<string, CreatePayOffFunction>::const_iterator

i = TheCreatorFunctions.find(PayOffId);

if (i == TheCreatorFunctions.end())

{

std::cout << PayOffId

<< " is an unknown payoff" << std::endl;

return NULL;

}