Boost C++ Libraries and Applications to Computational Finance

By DANIEL J. DUFFY

Part I – Introduction and Initial Examples

In this article we give an overview of the Boost C++ libraries (www.boost.org)  and we discuss their applicability to certain problems in computational finance. There are approximately 90 libraries in Boost at the moment of writing and each one is focused on a particular aspect of software development. Part of the mission statement on www.boost .org reads:

Boost provides free peer-reviewed portable C++ source libraries. We emphasize libraries that work well with the C++ Standard Library. Boost libraries are intended to be widely useful, and usable across a broad spectrum of applications. The Boost license encourages both commercial and non-commercial use.

We aim to establish “existing practice” and provide reference implementations so that Boost libraries are suitable for eventual standardization. Ten Boost libraries are already included in the C++ Standards Committee’s Library Technical Report (TR1) and will be in the new C++0x Standard now being finalized. C++0x will also include several more Boost libraries in addition to those from TR1. More Boost libraries are proposed for TR2.

The Boost libraries make use of templates and in this sense the generic programming model (GP) plays a more important role than the more traditional object-oriented programming (OOP) model. I have conveniently grouped the individual libraries into the following categories and I give some libraries in each category that I find useful:

• Higher-Order Functions: define function types as objects. We can view these libraries as a generalization of and improvement on the use of function pointers in C.  Libraries are: Function, Bind, Phoenix, Signals
• Data Types: data and number types: Tuple, Variant, Any, Rational, Integer, Interval. These are data types that improve the reliability or extend the functionality of existing C/C++ data types
• Data Structures: Complex structures that model geometrical and mathematical entities such as n-dimensional tensors, graphs and matrices (and matrix algebra)
• Maths: Univariate statistical distributions (about 25), special functions (e.g. gamma functions, Legendre polynomials), random number generation (e.g. Mersenne Twister), Accumulators (framework for statistical accumulators)
• Text and string processing: libraries that process, analyse and parse strings: StringAlgo, Tokenizer, Regex, Xpressive
• Utilities, miscellaneous and other libraries: many other useful libraries. We mention Smart Pointer, Thread, Date-Time, Serialization and Asio (portable networking)

In the past developers would have created their own C++ libraries for the above functionality but now that Boost has done it it relieves them of having to maintain proprietary libraries. Furthermore, it also gives them more time focus on their core business.

In subsequent articles we shall discuss some specific Boost libraries and their applications in computational finance:

• Solving n-factor Black Scholes PDE with multi_array and finite differences
• Numerical linear algebra routines with uBLAS
• Why boost Function is so powerful
• Statistical distributions and special functions

We conclude this introduction with some code. In this case we give an example from the Random library that provides generators for producing random numbers. Most Boost libraries are header only so as developer you do not normally have to worry about linking libraries into your project.

The example we give is C++ code to calculate the value of (approximately 3.1415 for most quants). Here is the full and well-documented source code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

#include <boost/random.hpp> 	// Convenience header file
#include <iostream>
#include <ctime>			// std::time
using namespace std;
 
#include <boost/limits.hpp>
 
int main()
{
	// A. Define the lagged Fibonacci and Uniform objects
	boost::lagged_fibonacci607 rng;
	rng.seed(static_cast<boost::uint32_t> (std::time(0)));
	boost::uniform_real<> uni(0.0,1.0);
 
	// B. Produce Uniform (0, 1)
	boost::variate_generator<boost::lagged_fibonacci607&, boost::uniform_real<> >
         uniRng(rng, uni);
 
	// C. Choose the desired accuracy
	cout << "How many darts to throw? "; long N; cin >> N;
 
	// D. Thow the dart; does it fall in the circle or outside
	// Start throwing darts and count where they land
	long hits = 0;
	double x, y, distance;
	for (long n = 1; n <= N; ++n)
	{
		x = uniRng(); y = uniRng();
		distance = sqrt(x*x + y*y);
 
		if ( distance <=1)
		{
			hits++;
		}
	}
 
	// E.  Produce the answer
	cout << "PI is: " << hits << ", " << 4.0 * double(hits) / double (N);
 
	return 0;
 
}

As an exercise, you can modify the above code to compute the probability that a quadratic equation has two real roots.

Part 2 – Functions, Functions, Functions

The Boost.Function Library

In this article we give an introduction to one of most applicable (and at the same time, one of the easiest to understand) libraries in Boost. This is the Boost.Function library that allows developers to define generic function objects which can then be instantiated and subsequently invoked. The focus in this article is to discuss the library and how it works. Future articles will extend the scope of the discussion by showing how to apply Boost.Function to computational finance applications (for example, PDE, Monte Carlo and optimisation). In general, the library allows us to create another level of indirection between an client code and the specific classes and functions that it uses. Some typical applications are:

• A Monte Carlo simulation that receives its random numbers from a variety of sources, for example Boost.Random (as discussed in Part I), the standard and unreliable rand() function or a propriety in-house random-number generator, or an open-source library (for example, Quantlib)
• An option pricing engine that interoperates with a range of numerical methods such as PDE/FDM, Monte Carlo and lattice models. The engine should have no knowledge of the internal workings of these numerical methods and it should use them solely by invoking their interfaces.

In general, real-life software is a combination of procedural, object-oriented and generic programming models and this raises a number of challenges when we develop software systems that need to interface with disparate code written using these models. In particular, the code we write should be able to work with C-style functions, function objects and member functions. To this end, Boost.Function provides the extra level of indirection between client code and the functions to which it delegates. As we shall see, the client code delegates to a Boost.Function object whose implementation or realization can be any one of a number of specific function types. In this sense Boost.Function is a Bridge (in the sense of Design Patterns) between client and implementation code.

The steps when using Boost.Function in applications are:

1. Specify signature the Boost.Function object (this entails its name, input arguments and return type).
2. Instantiate the object from Step 1 by assigning it to a function object, class member function or C-style function
3. Use the instantiated object in client code

Since many readers may not have hands-on experience with Boost we reduce the scope in this article by discussing a mini-application that serves as a model for more complex applications. In this case we create a class called MyClient that uses the exponential function exp(x), where x is a numeric type (typically a double, float or complex).
So, let’s get started. Since we are modelling all kinds of functions that accept a scalar value as input and return a scalar value we can specify this class of functions in Boost.Function as follows:

1
2

#include <boost/function.hpp>
typedef boost::function<double (double z)> ScalarFunction;

In this code we must use the Boost.Function library by including the appropriate files; second, as is common when working with templates we recommend using synonyms (using the typedef keyword). In this case we define a class called ScalarFunction that specifies a class of functions that have double as return type and double as input parameter. Similarly, the class of functions that accept two input parameters would read:

1

typedef boost::function<double (double z1, double z2)> TwoParameterFunction;

Thus, in order to understand the function signature, we read from left to right. We stress that the above specifications represent classes. We can now define instances (objects) of these classes as follows:

1
2

	ScalarFunction myFunc;
	TwoParameterFunction myFunc2;

We still cannot execute these function objects yet because they do not refer to any real code. In order to achieve this end we can assign to any kind of compatible function, that is one having the same signature as the corresponding class. In the above cases we assign the objects to standard C functions that are wrappers for the exponential function. The two functions are:

1
2
3
4
5
6
7
8
9
10
11

	// C function for exponential with 1 parameter
	double Exp1d(double z)
	{
		return exp(z);
	}
 
	// C function for exponential with 2 parameters
	double Exp2d(double z1, double z2)
	{
		return exp((z1+z2));
	}

Finally, we can assign the boost::function objects to these functions and we can subsequently execute the objects as it were because boost::function overloads the function call operator (). The following code shows how this is done:

1
2
3
4
5

	myFunc = Exp1d;
	cout << myFunc(1.0) << endl;
 
	myFunc2 = Exp2d;
	cout << myFunc2(1.0, 0.5) << endl;

This was our first encounter with Function. That was not too difficult!
We now discuss a more complicated example that shows how powerful Boost.Function is. To this end, we create a class that has an embedded algorithm, part of which is customizable because it has a boost::function as private member. The interface is:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

class MyClient
{ // A class with an embedded Boost.Function object. In
  // this sense it is policy-free as it does not depend on
  // a specific function form.
 
private:
 
	// Function specification: output/input
	ScalarFunction f;
 
public:
	MyClient(const ScalarFunction& myFunc) { f = myFunc; } 
 
	double calculate(double z) const
	{ // This simulates an algorithm
 
		return 2.0 * f(z) + sin(z*z); // The actual function called here
	}
 
	// Change the implementation of the calculator
	void SetFunction(const ScalarFunction& myFunc) { f = myFunc; }
};

In this case the precise form of the function f has not been specified; it is just a scalar function. This code is highly reusable, has no coupling with external base class references and such-like. In fact, you can create instances of this class by specializing f to be a C function, C++ function object (that is, one that implements the function call operator () ) or even a member function. This level of indirection is not possible with traditional object-oriented code. As a first example, we have:

1
2

	MyClient myC(Exp1d);
	cout <<  myC.calculate(x) << endl;

We can even change f to refer to another function:

1
2

	myC.SetFunction(AnotherFunc1d);
	cout << myC1.calculate(x) << endl;

where the function AnotherFunc1d is:

1
2
3
4

	double AnotherFunc1d(double z)
	{
		return exp(z) * sin(z);
	}

We now turn our attention to showing how instances of ScalarFunction are assigned to function objects. A function object is a class that implements the function call operator. An example is:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

class ExponentialDecorator
{ // My own version of exponential with overload ()
  // This is a function object
 
private:
	double a;
public:
	ExponentialDecorator(double param) { a = param;}
 
	double operator () (double z) const
	{ // This is the function call operator
 
		return exp(a*z);
	}
};

In this case we go beyond procedural functions because function objects can have state in addition to input arguments. An example in which we create two instances of ExponentialDecorator is:

1
2
3
4
5
6
7

	ExponentialDecorator myExp(20.0);
	MyClient myC2(myExp);
	cout << myC2.calculate(x) << endl;
 
	ExponentialDecorator myExp2(1.0);
	myC2.SetFunction(Exp1d);
	cout << myC2.calculate(x) << endl;

In this case we see that we can customize the behaviour of the algorithm calculate(x) by providing with different instances of the function object.
The last example discusses the approximation of the exponential function by Padé rational functions which are quotients of two polynomials of degree q (the numerator) and p (the denominator). The corresponding template class has parameters representing the degrees q and p and is a function object:

1
2
3
4

	template <int p, int q> struct PadeApproximation
	{
		double operator ()(double z) const;
	}

Well-known specializations of this class are:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

template <> struct PadeApproximation<1,0>
{ // Explicit Euler
 
	double operator ()(double z) const
	{
		return 1.0 + z;
	}
};
 
template <> struct PadeApproximation<0,1>
{ // Implicit Euler
 
	double operator () (double z) const
	{
		return 1.0/(1.0 - z);
	}
};
 
template <> struct PadeApproximation<1,1>
{ // Crank Nicolson
 
	double operator () (double z) const
	{
		return (2.0 + z)/(2.0 - z);
	}
};

As before, we can use these in client code with abandon:

1
2
3
4
5
6
7
8
9
10
11
12

	double x = -0.09;
	PadeApproximation<1,0> approxExplicitEuler;
	PadeApproximation<1,0> approxImplicitEuler;
	PadeApproximation<1,1> approxCrankNicolson;
 
	MyClient myEE(approxExplicitEuler);
	MyClient myIE(approxImplicitEuler);
	MyClient myCN(approxCrankNicolson);
 
	cout << myEE.calculate(x) << endl;
	cout << myIE.calculate(x) << endl;
	cout << myCN.calculate(x) << endl;

We have now completed this short introduction. Boost.Function includes more functionality than what we have described here. On the other hand, understanding the above syntax is most important in the short term.
In the next installment we shall show how we have integrated Boost Function, multi_array and the ADE finite difference method to price two-factor and three-factor option PDE problems.
The full source code is on my web site http://www.datasimfinancial.com/forum/viewtopic.php?t=111

Part 3 – Multi-dimensional Data Stuctures in Boost

In this blog we introduce the Boost Multidimensional Array Library (MAL) that allows developers to define multidimensional arrays in C++ applications. The MultiArray concept defines an interface to hierarchically nested containers. In particular, we can define n-dimensional data structures using the template classes provided by the library and it provides operations to create multiarrays, access their elements, navigating in multiarrays as well as creating view of multiarray data. The underlying flexible memory model supports a number of layouts, thus making the library suitable for a wide range of applications. We postpone a discussion of applications until we have introduced MAL and given some examples. You can download the code for this blog from my site http://www.datasimfinancial.com/forum/viewtopic.php?t=111 and test the code for yourself.
We now give a first example of a three-dimensional hypercube structure. First, we assume that Boost has been installed on your system and that all necessary search paths have been defined (this aspect is system and compiler-dependent). Next, we need to include the following header file in our code:

	#include <boost/multi_array.hpp>

The basic class that allows us to create multi-dimensional data structures is called multi_array. We have a number of options when creating an instance of this class; the difference lies in when the array’s so-called extents (the number of elements in each dimension of n-space) is defined. In the first case we create the array object and afterwards we can define its extents:

1
2
3
4
5
6
7
8
9
10

	// Define structure of tensor
	const int dim = 3; // 3d matrix
	typedef boost::multi_array<double, dim> Tensor;
	Tensor tensor;
 
	// Define the extents of each separate dimension
	const int NT = 3;
	const int NSIM = 4;
	const int NDIM = 3;
	tensor.resize(extents[NT][NSIM][NDIM]);

The second option is to define the extents directly in the constructor (note that have used a typedef which promotes code readability):

	// Useful shorthand
	Tensor tensor2(extents [NT][NSIM][NDIM]);

Finally, we can create an array by defining its extents in a collection:

1
2
3

	// Create a tensor using a Collections, in this case in Boost
	array<Tensor::index, dim> extentsList = { { NT, NSIM, NDIM } };
	Tensor tensor3 (extentsList);

Having created an array we obviously wish to access its data using indexing operators. We can choose between concatenated square brackets [] that accepts integers as index on the one hand and round brackets () that accept a collection as input on the other hand. We first define an alias that allows us to access the indices relating to an array:

1
2

	// Define a 3d loop to initialise the data
	typedef Tensor::index Index;

The following code shows how to initialize an array using square bracket notation:

1
2
3
4
5
6
7
8
9
10
11

	// Accessing elements using []
	for (Index i = 0; i != NT; ++i)
	{
		for (Index j = 0; j != NSIM; ++j)
		{
			for (Index k = 0; k != NDIM; ++k)
			{
				tensor[i][j][k] = i + j + k;
			}
		}
	}

We can also define an index collection and use round brackets to access the array elements (this option is more efficient than using square brackets in general):

1
2
3
4
5
6
7
8
9
10
11
12
13
14

	// Accessing elements using ()
	typedef boost::array<Tensor::index, dim> Indices;
 
	Indices indices;
	for (indices[2]=0; indices[2] < NZ; indices[2]++)
	{
		for (indices[0]=0; indices[0] < NX; indices[0]++)
		{
			for (indices[1]=0; indices[1] < NY; indices[1]++)
			{
				tensor(indices) = 0.0;
			}
		}
	}

This completes our first example. There is much more in the library in the way of whistles and bells but a discussion of these topics is outside the current scope.
As a last example, we wonder how a four-dimensional array can be created. In fact, it is not much more difficult than creating three-dimensional arrays. Here is an example:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

	// Typedef to make working with quads easier
	const int quadDim = 4;
	typedef boost::multi_array<double, quadDim> Quad;
 
	// The dimensions
	Quad::index dim1=3;
	Quad::index dim2=4;
	Quad::index dim3=2;
	Quad::index dim4=2;
 
	// Create a "quad" multi-array
	Quad quad(boost::extents[dim1][dim2][dim3][dim4]);
 
	// Fill the cube using the operator [index][index][index][index] syntax
	for (Quad::index d1=0; d1<dim1; d1++)
	{
		for (Quad::index d2=0; d2<dim2; d2++)
		{
			for (Quad::index d3=0; d3<dim3; d3++)
			{
				for (Quad::index d4=0; d4<dim4; d4++)
				{
					quad[d1][d2][d3][d4]= 0.0;
				}
			}
		}
	}

In the next blog I will give some applications of Boost.MultiArray to the numerical solution of partial differential equations (PDE) using the finite difference method (FDM). To this, the code in the current blog is sufficient to help us create finite solvers for a class of ADE schemes, which are unconditionally stable, second-order accurate and scale to higher dimensions.

Part 4 – Unconditionally stable and second-order accurate explicit Finite Difference Method for Three Factor PDE

The ADE Scheme meets the Boost Libraries

Summary and Goals

We now introduce the Alternating Direction Explicit ADE method to approximate the solution of the three-dimensional heat equation PDE in a bounded domain. ADE is a stable, second-order accurate finite difference scheme for general n-factor convection-diffusion equations (such as the Black Scholes PDE), it is easy to program and it is more efficient than conventional schemes such as Crank-Nicolson, ADI and Splitting.

The scope of this article is restricted to describing the ADE method, how we use it to solve a three-factor PDE and how we have implemented ADE using Boost. A future article will discuss the mathematical background to ADE and how we use it for approximating the solution of the n-factor Black-Scholes PDE.

The main goal of this note is to show how it is possible to solve high dimensional PDEs using a combination of the ADE method and Boost. In this sense we can model the problems using algebraic methods and we are closer to solving some problems associated with the curse of dimensionality.

The focus in this blog is on using Boost libraries with ADE; we exclude a discussion of the numerical analysis and this topic will be discussed elsewhere.

Background

One of the numerical techniques used in derivatives pricing is the Finite Difference Method (FDM) and it is used in Finance to approximate the solution of Partial Differential Equations (PDE) that describe the time evolution of a derivative that is contingent on one or more underlying variables (Duffy 2006). The theory of FDM for one and two-factor problems is reasonably well developed but out interest in this note is to discuss the application of the Alternating Direction Explicit (ADE) method to solving three factor derivatives PDE problems. In Duffy 2009 we discuss the applicability of ADE to one-factor option models where we motivate the method, discuss its mathematical and numerical foundations and we compare the scheme with other methods with regard to accuracy. We also give references to related work in Duffy 2009. Based on the results and successes with one-factor problems we have applied the method to multi-factor problems and in this blog we reduce the scope to the case of the three-dimensional heat equation as it contains many of the essential difficulties that we experience in multi-factor Black-Scholes PDEs.

Some of the advantages of using ADE for this class of PDE are:

• It is an unconditionally stable, second-order accurate and explicit finite difference scheme and this means that the method will give accurate results irrespective of the mesh sizes in the time and underlying (space) dimensions.
• ADE is efficient; in contrast to ADI and Splitting method we do not have to break a multi-factor PDE into a sequence of one-dimensional problems and in contrast to these latter methods we do not have to solve a tridiagonal matrix system at each time level and for each space dimension. The scheme is also very easy to design and to implement.
• We have applied ADE to multi-factor Black-Scholes PDEs using previously developed methods such as exponential fitting, the Janenko scheme for problems with mixed derivatives and the ability to transform the PDE in semi-infinite space to a PDE that is defined on a unit hypercube (Duffy 2009).
• It is easy to parallelise ADE by inserting compiler directives from the OpenMP library into the C++ code that implements the scheme, for example by using the parallel sections construct.
• ADE resolves the curse of dimensionality: this refers to the problems that we encounter when we try to design and implement finite difference schemes in three or more dimensions. Traditional finite difference schemes break down due to the complexity. The ADE scheme – in combination with Boost.MultAarray – resolves this curse to a large extent.

The Model Problem

In this section we discuss how to compute the solution of three-dimensional heat equation in a unit cube. The heat equation is one of the most important equations in mathematical physics and it has applications in many domains.

For this reason it is important to know how to find the solution (either in exact form or by numerical means). In this case we choose the finite difference method (FDM) to approximate the solution of the heat equation at equidistant mesh points in three-dimensional space. We assume that the reader knows what partial derivatives, partial differential equations (PDE) and divided differences are.

The PDE for the heat equation in three dimensions is given by:

(1)

There are an infinite number of solutions to (1). In order to specify a unique solution to (1) we must define some auxiliary conditions. First, we define the initial condition at t = 0:

(2)

In this case we state that the unknown solution u is known at t = 0. Second, we define zero Dirichlet boundary conditions when the independent variables x, y and z are 0 and 1:

(3)

We have now defined the initial boundary-value problem (1), (2), (3) defined on the unit cube and for all time t in the interval (0, T), where T > 0.

Turning to the finite difference solution of (1), (2), (3) we partition cube [0,1]^3 into three-dimensional mesh points.
To this end, we partition the interval [0,1] in each of the directions x, y and z into equal subintervals and we create three arrays as follows:

where the mesh sizes are defined by:

.

We also discretise the interval (0, T) into a number of equidistant mesh points in time as follows:

Finally, we shall be working with mesh functions that are defined at the discrete mesh points in space and time:

In all further discussions we assume that all functions are mesh functions of this form.

The ADE Scheme

We are now ready to introduce the Alternating Direction Explicit (ADE) method (Duffy 2009). The scheme is stable, has second-order accuracy and it is explicit which means that we compute a solution without having to solve a matrix system at each time level. In general, we compute the final solution in a series of steps. First, we compute one component from the lower apex (0,0,0) to the upper apex (1,1,1) of the cube as follows. We call it the upward sweep:

(4)

The second component computes a solution from the upper apex to the lower apex. We call it the downward sweep:

where are the mesh sizes in the x, y and z directions

(5)

Having computed these solutions (and this can be done in parallel) we arrive at the final solution as the average of the two components:

(6)

In order to unambiguously define the solution (6) we define the discrete boundary conditions:

(7)

and the discrete boundary conditions for the component :

(8)

with the same boundary conditions holding for . Summarising, the scheme defined by equations (4) to (8) is an unambiguous specification of the finite difference method for the heat equation. We shall map this scheme using multi_array. To this end, we use some algebra to cast equations (4) and (5) to computable form:

(9)

and

(10)

where

(11)

Design and Implementation in C++

We now discuss how to set up the software system for solving the model PDE using the ADE scheme. The code that we show presently is based on a larger ongoing software project for the three-factor Black-Scholes equation. We have modified it for pedagogical reasons to make it readable.

The design approach is based on the strategy of decomposing the problem at hand into loosely coupled and autonomous (cohesive) subsystems and classes. The main classes that we model are for:

• The coefficients of PDE (class ThreeFactorPde).
• The PDE domain in conjunction with its boundary and initial conditions (class ThreeFactorPdeDomain).
• The class that implements the ADE scheme (ThreeFactorADESolver).

We now discuss each class in detail. The class ThreeFactorPde models the PDE (1). In this case the class is slightly more general because it supports general diffusion equations and coefficients. The class interface uses boost::function to model all the coefficients in equation (1):

1
2
3
4
5
6
7
8
9
10
11
12
13
14

	template <typename T> class ThreeFactorHeatPde
	{ // The parameter 'T' is usually instantiated to double in client code
 
	public:
 
		// Coefficients of the elliptic part of PDE
		boost::function<T (T,T,T,T)> a11;
		boost::function<T (T,T,T,T)> a22;
		boost::function<T (T,T,T,T)> a33;
 
		boost::function<T (T,T,T,T)> F;	// RHS
 
		ThreeFactorHeatPde() {}
	};

We now need to define the region in which the PDE (1) is defined and to this end we define it as three ranges (a range is an interval in the mathematical sense) objects. We also define the Dirichlet boundary conditions on the six faces of the unit hypercube and the associated initial condition, examples of which are seen in equations (2) and (3):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

	template <typename T> struct ThreeFactorPdeDomain
	{
		//	1. Domain
		Range<T> rx;
		Range<T> ry;
		Range<T> rz;
		Range<T> rt;
 
		//	2. Boundary conditions, anticlockwise
 
		// Specific boundaries in x
		boost::function<T (T y, T z, T t)> XLowerBC;
		boost::function<T (T y, T z, T t)> XUpperBC;
 
		// Specific boundaries in y
		boost::function<T (T x, T z, T t)> YLowerBC;
		boost::function<T (T x, T z, T t)> YUpperBC;
 
		// Specific boundaries in z
		boost::function<T (T x, T y, T t)> ZLowerBC;
		boost::function<T (T x, T y, T t)> ZUpperBC;
 
		//	3. Initial condition
		boost::function<T (T x,T y, T z)> IC;
	};

We need to instantiate these two classes because we are solving the three-dimensional heat equation with specific boundary and initial conditions. To this end, we define a dedicated namespace called ThreeFactorHeatImpl that defines all the parameters and functions needed to define a particular concrete instantiation of an initial boundary problem for the heat equation. Its interface is:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103

namespace ThreeFactorHeatImpl
{
 
	// Domain information
	double T;
	double xMax, yMax, zMax;	// Truncated domain
 
	double diffusionX(double x, double y, double z, double t)
	{
		return 1.0;
	}
 
	double diffusionY(double x, double y, double z, double t)
	{
		return 1.0;
	}
 
	double diffusionZ(double x, double y, double z, double t)
	{
		return 1.0;
	}
 
	double H(double v)
	{
		return v*(1.0 - v);
	}
 
	double phi(double x, double y, double z)
	{
		return exp(-x*x) * exp(-y*y) * exp(-z*z);
	}
 
	double termInhomogeneous(double x, double y, double z, double t)
	{
		return (+6.0 - 4.0*(x*x+y*y+z*z)) * phi(x,y,z);
	}
 
	// Factories
	void createThreeFactorHeatPde(ThreeFactorHeatPde<double>& pde)
	{
		// Initialise all functions
		pde.a11 = &ThreeFactorHeatImpl::diffusionX;
		pde.a22 = &ThreeFactorHeatImpl::diffusionY;
		pde.a33 = &ThreeFactorHeatImpl::diffusionZ;
		pde.F = &ThreeFactorHeatImpl::termInhomogeneous;
	}
 
	double IC(double x, double y, double z)
	{
		return phi(x,y,z);
	}
 
	double BCXLower(double y, double z, double t)
	{
		return  exp(-y*y) * exp(-z*z);
	}
	double BCXUpper(double y, double z, double t)
	{
		return  exp(-xMax*xMax) * exp(-y*y) * exp(-z*z);
	}
 
	double BCYLower(double x, double z, double t)
	{
		return  exp(-x*x) * exp(-z*z);
	}
 
	double BCYUpper(double x, double z, double t)
	{
		return  exp(-x*x) * exp(-yMax*yMax) * exp(-z*z);
	}
 
	double BCZLower(double x, double y, double t)
	{
		return  exp(-x*x) * exp(-y*y);
	}
 
	double BCZUpper(double x, double y, double t)
	{
		return  exp(-x*x) * exp(-y*y) * exp(-zMax*zMax);
	}
 
	// Factory
	void createThreeFactorPdeDomain(
			ThreeFactorPdeDomain<double>& pdeDomain)
	{
		pdeDomain.rx = Range<double>(0.0, xMax);
		pdeDomain.ry = Range<double>(0.0, yMax);
		pdeDomain.rz = Range<double>(0.0, zMax);
		pdeDomain.rt = Range<double>(0.0, T);
 
		pdeDomain.XLowerBC = &ThreeFactorHeatImpl::BCXLower;
		pdeDomain.XUpperBC = &ThreeFactorHeatImpl::BCXUpper;
 
		pdeDomain.YLowerBC = &ThreeFactorHeatImpl::BCYLower;
		pdeDomain.YUpperBC = &ThreeFactorHeatImpl::BCYUpper;
 
		pdeDomain.ZLowerBC = &ThreeFactorHeatImpl::BCZLower;
		pdeDomain.ZUpperBC = &ThreeFactorHeatImpl::BCZUpper;
 
		pdeDomain.IC = &ThreeFactorHeatImpl::IC;
	}
 
} // End of namespace

We now describe the class that models the ADE finite difference scheme. The scheme is a one-step time-marching one and this means that the information at time level n is used to compute the information at time level n+1. In particular, the upsweep and downsweep vectors defined by equations (9) and (10), respectively and their averaged value defined by equation (6) are represented in C++ as Boost multi-arrays:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

class ThreeFactorADESolver
{ // Using ADE method for 3d heat equation PDE
 
private:
	ThreeFactorHeatPde<double> pde;	// The pde coefficients
	ThreeFactorPdeDomain<double> pdeDomain;	// Domain, BC, IC
 
	// Data structures
	BoostTensor* U;	// upper sweep, n+1
	BoostTensor* UOld;	// upper sweep, n
	BoostTensor* V;	// lower sweep, n+1
	BoostTensor* VOld;	// lower sweep, n
 
public:
	BoostTensor* MatNew;	// averaged solution, level n+1
private:
	// Mesh-related data
	double hx, hy, hz, delta_k, 	// Step lengths in each direction
		 hx2, hy2, hz2,	// 1/h^2
		 Sum,	A11, A22, A33,	// delta_k * (hx2 + hy2 + hz2)
		 G;	// ~F	// Right side of FDM for u_t = Lu + F
 
	// Other variables
	double tprev, tnow, T;	// Time levels
	long NX, NY,NZ,NT;	// Number of subdivisions
	double t1, t2, t3;	// mesh values in each (x,y,z) direction
 
	double coeffI, coeffII;
 
public:
	// Mesh arrays in each direction
	Vector<double, long> xmesh;
	Vector<double, long> ymesh;
	Vector<double, long> zmesh;
	Vector<double, long> tmesh;
};

where we have used the shorthand notation:

1
2
3

	int const Dim = 3;
	typedef boost::multi_array<double, Dim> BoostTensor;
	typedef boost::array<BoostTensor::index, 3> BoostTensorIndices;

The code implementing the equations (10), (11) and (6) is:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

for (BoostTensor::index i = 1; i <= NX-1; i++)
{
	for (BoostTensor::index j = 1; j <= NY-1; j++)
	{
		for (BoostTensor::index k = 1; k <= NZ-1; k++)
		{
 
			t1 = xmesh[i]; t2 = ymesh[j]; t3 = zmesh[k];
			G = pde.F(t1,t2,t3,tnow) * delta_k;
 
			(*U)[i][j][k] =  (
					(*UOld)[i][j][k]*coeffI
					+ A11*((*UOld)[i+1][j][k] + (*U)[i-1][j][k])
					+ A22*((*UOld)[i][j+1][k] + (*U)[i][j-1][k])
					+ A33*((*UOld)[i][j][k+1] + (*U)[i][j][k-1])
					+ G
					  )/coeffII;
		}
	}
}
for (BoostTensor::index i = NX-1; i >= 1; i--)
{
	for (BoostTensor::index j = NY-1; j >= 1; j--)
	{
		for (BoostTensor::index k = NZ-1; k >= 1; k--)
		{
 
			t1 = xmesh[i]; t2 = ymesh[j]; t3 = zmesh[k];
			G = pde.F(t1,t2,t3,tnow) * delta_k;
 
			(*V)[i][j][k] = (
					(*VOld)[i][j][k]*coeffI
					+ A11*((*V)[i+1][j][k] + (*VOld)[i-1][j][k])
					+ A22*((*V)[i][j+1][k] + (*VOld)[i][j-1][k])
					+ A33*((*V)[i][j][k+1] +  (*VOld)[i][j][k-1])
					+ G
					)/coeffII;
		}
	}
}
for (BoostTensor::index i = 0; i <= NX; i++)
{
 
	for (BoostTensor::index j = 0; j <= NY; j++)
	{
		for (BoostTensor::index k = 0; k <= NZ; k++)
		{
			(*MatNew)[i][j][k] = 0.5*((*U)[i][j][k] + (*V)[i][j][k]);	// == u
 
			(*UOld)[i][j][k] = (*U)[i][j][k];
			(*VOld)[i][j][k] = (*V)[i][j][k];
		}
	}
}

Finally, we show how to create a program to show the application of ADE in C++:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

int main()
{
 
	using namespace ThreeFactorHeatImpl;
	xMax =1.0;
	yMax = 1.0;
	zMax = 1.0;
	T = 2.0;
 
	// Create pde and domain
	ThreeFactorHeatPde<double> pde;
	createThreeFactorHeatPde(pde);
	ThreeFactorPdeDomain<double> pdeDomain;
	createThreeFactorPdeDomain(pdeDomain);
	cout << "Now creating solvern";
 
	// FD scheme
	long V = 100;
	long NXX = V; long NYY = V; long NZZ = V; long NTT = V;
	DatasimClock myClock;
	myClock.start();
 
	ThreeFactorADESolver solver(pde, pdeDomain, NXX, NYY, NZZ, NTT);
	cout << "Now starting computationn";
	solver.result();
 
	myClock.stop();
	cout << endl << "Running time " << myClock.duration() << endl;
 
	BoostTensor exactTensor = CreateDiscreteFunction(&ExactSolution,
					solver.xmesh, solver.ymesh, solver.zmesh,
					T);
 
	HotSpot maxError = MaxNorm(exactTensor, *(solver.MatNew));
 
	// Get the spot where the largest error occurs
	cout << endl << "Max error " << maxError.get<0>() << ", "
			<< maxError.get<1>() << ", " << maxError.get<2>() << ", " << maxError.get<3>() << endl;
 
	return 0;
}

In this example, we have tested the scheme against the following exact solution:

1
2
3
4

double ExactSolution(double x, double y, double z, double t)
{
	return exp(-x*x) * exp(-y*y) * exp(-z*z);
}

and in order to compare the exact and approximate values we have created a utility function called MaxNorm()that returns a 3-tuple with fields for the values of the biggest error and where this error occurs:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

typedef tuple<double, BoostTensor::index, BoostTensor::index, BoostTensor::index> HotSpot;
 
HotSpot MaxNorm(const BoostTensor& tensor, const BoostTensor& tensor2)
{ // Compute largest element of m1 - m2 and state where it occurs
 
	// For readability, define start and end indices for each dimension
	BoostTensor::size_type size0 = tensor.shape()[0];
	BoostTensor::index start0 = tensor.index_bases()[0];
	BoostTensor::index end0 = start0 + size0;
 
	BoostTensor::size_type size1 = tensor.shape()[1];
	BoostTensor::index start1 = tensor.index_bases()[1];
	BoostTensor::index end1 = start1 + size1;
 
	BoostTensor::size_type size2 = tensor.shape()[2];
	BoostTensor::index start2 = tensor.index_bases()[2];
	BoostTensor::index end2 = start2 + size2;
 
	double val = fabs(tensor[start0][start1][start2]
			- tensor2[start0][start1][start2]);
	double tmp;
 
	HotSpot result;
 
	for (BoostTensor::index row=start0; row<end0; row++)
	{
		for (BoostTensor::index column=start1; column<end1; column++)
		{
			for (BoostTensor::index layer=start2; layer<end2; layer++)
			{
				tmp = fabs(tensor[row][column][layer]
					- tensor2[row][column][layer]);
				if (tmp > val)
				{
					val = tmp;
 
					// Update the hotspot
					result.get<0>() = val;
					result.get<1>() = row;
					result.get<2>() = column;
					result.get<3>() = layer;
				}
			}
		}
	}
 
	return result;
}

We conclude this note with some observations on the accuracy and efficiency of the ADE scheme as applied to the current problem. We have carried out a number of tests with a range of boundary and initial conditions. In general, the scheme is very accurate and fast. Even with NX = NY = NZ = NT = 50 we get O(e-5) accuracy and the total computation time is approximately 8 seconds on a laptop. We report on our findings for the Black Scholes PDE in the next blog.

Summary

This note builds on my previous two previous Quantnet.com blogs on the Boost libraries for higher-order functions and multi-dimensional arrays. We were particularly interested in showing how the ADE (Alternating Direction Explicit) method is applied to approximating the solution of the three-dimensional heat equation and implementing this scheme using Boost libraries.
The results from this feasibility study suggest our continuation of research into linear and nonlinear multi-factor models for equity and interest rate derivatives problem. We hope to report on these issues at a later stage.

References

• Duffy, D.J. 2006 Finite Difference methods in financial engineering John Wiley and Sons Chichester UK.
• Duffy, D.J. 2009 Unconditionally stable and second-order accurate explicit Finite Difference Schemes using Domain Transformation Part I. One-Factor Equity Problems. SSRN http://papers.ssrn.com/sol3/papers.cfm?abstract_id=1552926

About the Author:

Daniel Duffy is an author and trainer. His company Datasim specializes in methods and techniques for solving problems in quantitative finance. He is the author of Monte Carlo Frameworks: Building Customisable High-performance C++ Applications and Introduction to C++ for Financial Engineers: An Object-Oriented Approach. For more information on the author, see Quant Network’s interview with Daniel Duffy