. 5
( 9)


double Time,
double DiscountedFutureValue)
virtual ˜TreeAmerican(){}

PayOffBridge ThePayOff;

The header for the European option is very similar:
Listing 8.4 (TreeEuropean.h)

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

class TreeEuropean : public TreeProduct

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

PayOffBridge ThePayOff;
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),

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


Listing 8.6 (TreeEuropean.cpp)

#include <TreeEuropean.h>
#include <minmax.h>

TreeEuropean::TreeEuropean(double FinalTime,
const PayOffBridge& ThePayOff_)
: TreeProduct(FinalTime),
8.4 A tree class 129


double TreeEuropean::FinalPayOff(double Spot) const
return ThePayOff(Spot);

double TreeEuropean::PreFinalValue(double,
//Spot, Borland compiler
//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

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

Listing 8.7 (BinomialTree.h)
#pragma warning( disable : 4786 )

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

class SimpleBinomialTree
SimpleBinomialTree(double Spot_,
const Parameters& r_,
const Parameters& d_,
double Volatility_,
unsigned long Steps,
double Time);

double GetThePrice(const TreeProduct& TheProduct);

void BuildTree();

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

SimpleBinomialTree::SimpleBinomialTree(double Spot_,
const Parameters& r_,
const Parameters& d_,
double Volatility_,
unsigned long Steps_,
double Time_)
: Spot(Spot_),

void SimpleBinomialTree::BuildTree()
TreeBuilt = true;

double InitialLogSpot = log(Spot);

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


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

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&
if (!TreeBuilt)

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 =

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][k].second =

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

#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,
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,
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 =
cout << "forward price" << actualForwardPrice << "\n";

Steps++; // now redo the trees with one more step
SimpleBinomialTree theNewTree(Spot,rParam,dParam,Vol,

double euroNewPrice =
double americanNewPrice =

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)


#include <PayOff3.h>
class PayOffForward : public PayOff
PayOffForward(double Strike_);

virtual double operator()(double Spot) const;
virtual ˜PayOffForward(){}
virtual PayOff* clone() const;

double Strike;
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.)

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

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


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

class BSCall


BSCall(double r_, double d_,
double T, double Spot_,
double Strike_);

double operator()(double Vol) const;


double r;
double d;
double T;
double Spot;
double Strike;

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

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

if (y < Target)
Low = x;

if (y > Target)
High = x;

x = 0.5*(Low+High);

y = TheFunction(x);
( (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)
#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 =
double PriceTwo =

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


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

Listing 9.6 (BSCallTwo.h)

#ifndef BS_CALL_TWO_H
#define BS_CALL_TWO_H

class BSCallTwo

BSCallTwo(double r_, double d_,
double T, double Spot_,
double Strike_);

double Price(double Vol) const;
double Vega(double Vol) const;

double r;
double d;
double T;
double Spot;
double Strike;
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_)

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)
#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 =

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

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.

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

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

#include <PayOff3.h>

#if defined(_MSC_VER)
#pragma warning( disable : 4786)

#include <map>
#include <string>

class PayOffFactory
typedef PayOff* (*CreatePayOffFunction)(double );

static PayOffFactory& Instance();
void RegisterPayOff(std::string, CreatePayOffFunction);
PayOff* CreatePayOff(std::string PayOffId,double Strike);

std::map<std::string, CreatePayOffFunction>
PayOffFactory(const PayOffFactory&){}
PayOffFactory& operator=
(const PayOffFactory&){ return *this;}

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

#include <PayOffFactory.h>
#include <iostream>
using namespace std;

void PayOffFactory::RegisterPayOff(string PayOffId,
CreatePayOffFunction 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;


. 5
( 9)