Boost C++ Libraries and Applications to Computational Finance

[Daniel 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 pi (approximately 3.1415 for most quants). Here is the full and well-documented source code:

#include  // Convenience header file
#include
#include // std::time
using namespace std;

#include

int main()
{
// A. Define the lagged Fibonacci and Uniform objects
boost::lagged_fibonacci607 rng;
rng.seed(static_cast (std::time(0)));
boost::uniform_real<> uni(0.0,1.0);

// B. Produce Uniform (0, 1)
boost::variate_generator >
        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:

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

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:

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:

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

 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:

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:

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

We can even change f to refer to another function:

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

where the function AnotherFunc1d is:

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:

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:

 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:

template  struct PadeApproximation
{
double operator ()(double z) const;
}

Well-known specializations of this class are:

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:

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:

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

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

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

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

// Accessing elements using ()
typedef boost::array 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:

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

(\frac{\partial u}{\partial t} = \triangle u \equiv \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2} + F(x,y,z,t), \quad x,y,z \in (0,1), \quad t > 0)
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)​

(u(x,y,z,0) = \varphi (x,y,z), \quad x,y,z \in (0,1))
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)​

(u(0,y,z)=u(1,y,z)=0, \quad y,z \in (0,1))
(u(x,0,z)=u(x,1,z)=0, \quad x,z \in (0,1))
(u(x,y,0)=u(x,y,1)=0, \quad x,y \in (0,1))
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:
(x_i=ih_1, \quad 0 \leq i \leq I)
(y_j=jh_2, \quad 0 \leq j \leq J)
(z_k=kh_3, \quad 0 \leq k \leq K)
where the mesh sizes are defined by:
(h_1=1/I, \quad h_2=1/J, \quad h_3=1/K, I, J, K \geq 2) are integrers.
We also discretise the interval (0, T) into a number of equidistant mesh points in time as follows:
(\triangle t = T/NT, \quad t_n = n \triangle t, \quad 0 \leq n \leq NT)
Finally, we shall be working with mesh functions that are defined at the discrete mesh points in space and time:
(v^n_{i,j,k} \sim v(x_i,y_j,z_k,t_n))
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)​

(\frac{v^{n+1}_{i,j,k} - v^n_{i,j,k}}{\triangle t} = \frac{1}{h^2_1}(V^{n+1}_{i+1,j,k}-V^n_{i,j,k}-V^{n+1}_{i,j,k}+V^n_{i-1,j,k}))
(+\frac{1}{h^2_1}(V^{n+1}_{i,j+1,k}-V^n_{i,j,k}-V^{n+1}_{i,j,k}+V^n_{i,j-1,k}))
(+\frac{1}{h^2_3}(V^{n+1}_{i,j,k+1}-V^n_{i,j,k}-V^{n+1}_{i,j,k}+V^n_{i,j,k-1}))
(1 \leq i \leq I-1, \quad 1 \leq j \leq J-1, \quad 1 \leq k \leq K-1, \quad n \geq 0)
The second component computes a solution from the upper apex to the lower apex. We call it the downward sweep:

(5)​

(\frac{v^{n+1}_{i,j,k} - v^n_{i,j,k}}{\triangle t} = \frac{1}{h^2_1}(V^{n+1}_{i+1,j,k}-V^n_{i,j,k}-V^{n+1}_{i,j,k}+V^n_{i-1,j,k}))
(+\frac{1}{h^2_1}(V^{n+1}_{i,j+1,k}-V^n_{i,j,k}-V^{n+1}_{i,j,k}+V^n_{i,j-1,k}))
(+\frac{1}{h^2_3}(V^{n+1}_{i,j,k+1}-V^n_{i,j,k}-V^{n+1}_{i,j,k}+V^n_{i,j,k-1}))
(1 \leq i \leq I-1, \quad 1 \leq j \leq J-1, \quad 1 \leq k \leq K-1, \quad n \geq 0) where (h_{1}, h_{2}, h_{3}) are the mesh sizes in the x, y and z directions

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

(u^n_{i,j,k} = \frac{1}{2}(U^n_{i,j,k}+V^n_{i,j,k}) \quad 0 \leq i \leq I, \quad 0 \leq j \leq J, \quad 0 \leq k \leq K, \quad n \ge 0)
In order to unambiguously define the solution (6) we define the discrete boundary conditions:

(7)​

(U^0_{i,j,k}= V^0_{i,j,k}=\varphi_{i,j,k} \equiv \varphi(x_i,y_j,z_k), \quad 1 \leq i \leq I-1,\quad 1 \leq j \leq J-1,\quad 1 \leq k \leq K-1)
and the discrete boundary conditions for the component (U^n_{i,j,k}) :

(8)​

(\left.\begin{aligned} U_{0,j,k}^n = U_{I,j,k}^n = 0, \; 0 \leq j \leq J, \; 0 \leq k \leq K\\ U_{i,0,k}^n = U_{i,J,k}^n = 0, \; 0 \leq i \leq I, \; 0 \leq k \leq K\\ U_{i,j,0}^n = U_{i,j,K}^n = 0, \; 0 \leq i \leq I, \; 0 \leq j \leq J \end{aligned}\right\rbrace n \geq 0)
with the same boundary conditions holding for (V^n_{i,j,k}). 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)​

(U^{n+1}_{i,j,k}(1+\sigma)=U^{n}_{i,j,k}(1-\sigma)+\lambda_1(U^{n+1}_{i+1,j,k}+U^{n}_{i-1,j,k})+\lambda_2(U^{n+1}_{i,j+1,k}+U^{n}_{i,j-1,k})+\lambda_3(U^{n+1}_{i,j,k+1}+U^{n}_{i,j,k+1})+\Delta t F^{n+1}_{i,j,k})
and

(10)​

(V^{n+1}_{i,j,k}(1+\sigma)=V^{n}_{i,j,k}(1-\sigma)+\lambda_1(V^{n+1}_{i+1,j,k}+V^{n}_{i-1,j,k}) +\lambda_2(V^{n+1}_{i,j+1,k}+V^{n}_{i,j-1,k})+\lambda_3(V^{n+1}_{i,j,k+1}+V^{n}_{i,j,k+1})+\Delta t F^{n+1}_{i,j,k})
where

(11)​

(\lambda_j = \triangle t/h^2_j, \quad j=1,2,3,\quad \sigma= \sum_{j=1}^3 \lambda_j)
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):

template  class ThreeFactorHeatPde
{ // The parameter 'T' is usually instantiated to double in client code

public:

// Coefficients of the elliptic part of PDE
boost::function a11;
boost::function a22;
boost::function a33;

boost::function 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):

template  struct ThreeFactorPdeDomain
{
// 1. Domain
Range rx;
Range ry;
Range rz;
Range rt;

// 2. Boundary conditions, anticlockwise

// Specific boundaries in x
boost::function XLowerBC;
boost::function XUpperBC;

// Specific boundaries in y
boost::function YLowerBC;
boost::function YUpperBC;

// Specific boundaries in z
boost::function ZLowerBC;
boost::function ZUpperBC;

// 3. Initial condition
boost::function 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:

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& 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& pdeDomain)
{
pdeDomain.rx = Range(0.0, xMax);
pdeDomain.ry = Range(0.0, yMax);
pdeDomain.rz = Range(0.0, zMax);
pdeDomain.rt = Range(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:

class ThreeFactorADESolver
{ // Using ADE method for 3d heat equation PDE

private:
ThreeFactorHeatPde pde; // The pde coefficients
ThreeFactorPdeDomain 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 xmesh;
Vector ymesh;
Vector zmesh;
Vector tmesh;
};

where we have used the shorthand notation:

 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:

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

int main()
{

using namespace ThreeFactorHeatImpl;
xMax =1.0;
yMax = 1.0;
zMax = 1.0;
T = 2.0;

// Create pde and domain
ThreeFactorHeatPde pde;
createThreeFactorHeatPde(pde);
ThreeFactorPdeDomain 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:

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:

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

 

[sagrailo]

The first half of the code is showing exactly what I hate the most about Boost, and that is huge verbosity and over-complication: I guess when one is at point to actually understand what is going on in the first 5 statements of main() function, he will already know enough about RNGs that it's going to be rather tempting to write one on his own, and thus just got rid of dealing with Boost crap.

 

[Daniel Duffy]

I grant you the names of these classes get some time to get used to (I would have preferred LaggedFibonacci myself). There is very little I can do about the syntax.

And it takes some effort to understand the template code.

But the code is needed as it is, and a good knowledge of templates is needed.

And writing your own RNG is laudable but not realistic because it's all been done before.

The second half of the code was easier?

 

[sagrailo]

There exist some good part in Boost (shared pointers, regex handling, threads), and these are mostly going to be pulled into the next version of C++ standard library (actually, most modern compilers already support most of the std::tr1 stuff, so I see no reason to use Boost for this type of functionality). The rest of Boost is, in my experience, utter crap: for example, take a look into boost_program_options: anyone that ever used getopt() will plain tell you that boost_program_options is an insane over-complication, and the most of the Boost is just like that. Coupled with sparse/wrong documentation, and awful performance in cases when performance would be expected (I'm talking about Boost uBLAS here), it doesn't make convincing case for me. Worse-over, Boost is pushing programmers towards heavily use of templates and template meta-programming, and this is seemingly even some kind of norm among the people using it; now, imagine that you have to track a bug in Boost RNG code, and try to follow through all of various header files that you would have to read in order to understand what's going on in the mentioned first 5 statements of the code above - God help you if you find yourself in the situation like that (I did, whenever I was forced to use Boost, and for these memories I wouldn't touch it with a ten-foot pole any more).

 

[Daniel Duffy]

Fair enough. You do have a number of valid points. But you obviously have had a bad experience with Boost.

I have some small test cases to give an impression (101 examples) of what's in Boost. I created it precisely to try making it a bit more accessible.

http://www.datasimfinancial.com/forum/viewtopic.php?t=111

Some remarks:

1. As normal developer you need not worry about meta templates. It's a myth imho that MPL is being pushed.
2. String Algo does a lot of what Regex does, and more. And Xpressive is interesting as well.
3. I wrestled with the syntax of Random for a few days until I got the hang of it. Then it's easy. It's just math. The alternative is to build OO hierarchies for RNG with consequence of class explosion and less efficient etc. Random is particularly easy to fathom; it's just maths but the RNG names are cryptic to be sure
4. Yes, I spend time in header files, when necessary. On the other hand, I try to understand one case and take it from there.
5. Smart ptr and Thread are plumbing; the interesting ones are useful in scientific and finance applications

There is a possibility that I may revise the above list but at the moment I am on the whole content with the results till now.

Out of interest, what kinds of applications did you apply Boost to?


"imagine that you have to track a bug in Boost RNG code"
What kind of bug? accuracy, run-time, compile-time??

 

[sagrailo]

I was talking about using Boost in the context of implementing options pricing code, that was the project that I was involved with. With regards to mentioning bugs in RNG code: it's not that I've actually encountered any kind of bug there (haven't used RNG Boost library much - if you're serious with regards to this type of work, kind of RNG implementations as these offered by Boost are simply not performant enough, and you'll anyway have to implement RNG on your own, to fully utilize hardware your code is executing on), but I was more trying to point out that once when you deal with code involving several templates, parametrized say by several other classes from some class hierarchies, it gets very hard to find where the code executed actually comes from, even for such banal purpose as to know say where to put some printout statements for the purpose of debugging, or understanding better the code execution paths.
On the other side: while the overall impression certainly heavily depends on one's previous background, I'd say there exist other C++ "foundation" class libraries, with the purpose somewhat similar to Boost, that are demonstrating that it is possible to build this kind of library for C++ that is much easier to understand and use, performing better than Boost, and being much better documented - for example, I had great success with using Qt toolkit for much of this kind of stuff often needed in C++ code, but missing from C++ standard library.

 

[Daniel Duffy]

Just a few words on RNG in Random, in particular the popular Mersenne Twister. According to the documentation the original Matsumoto/Nishimura has been rewritten and checked to give the same output as before.
There are also performance indicators available; MT has relative speed of 70% compared to the fastest (which is lagged Fibonacci).
One advantage of use is that we do not have to use older and less reliable RNG such as rand() and Box Muller.

 

[missingpa]

I was wondering why you take the square root on line 29. That is an expensive and not needed in order to know if a dot is inside the circle.

 

[Daniel Duffy]

"I was wondering why you take the square root on line 29. That is an expensive and not needed in order to know if a dot is inside the circle."

In this case you can optimise this step, indeed. But the sqrt(..) represents the radius and not its square, so the code for radius R would be

sqrt() <= R

I would prefer then

x^2 + y^2 <= 1*1

 

[Daniel Duffy]

In my first article I mentioned the exercise of finding the probability of finding real roots of ax^2 + bx + c = 0 where a, b and c are standard normal variates (If I remember correctly, the value is .646... approx. which can be found by Edelman and Kostlan formula). You can run the code as a first check how accurate RNG in Boost is:

//
boost::normal_distribution<> nor(0,1);
boost::lagged_fibonacci607 rng;
rng.seed(static_cast<boost::uint32_t> (std::time(0)));
// Produce Normal (0, 1)
boost::variate_generator<boost::lagged_fibonacci607&, boost::normal_distribution<> >
norRng(rng, nor);
cout << "How many darts to throw? "; long N; cin >> N;
double a, b, c, factor;
long hits = 0;
// The dart board
for (long n = 1; n <= N; ++n)
{
// Generate 3 N(0,1) variates
a = norRng(); b = norRng(); c = norRng();
factor = b*b - 4.0*a*c;
if ( factor >= 0.0)
{
hits++;
}
}
cout << "Prob of positive roots: " << hits << ", " << double(hits) / double (N);

the rest of the code is the same as before.

 

[Daniel Duffy]

If you like a small exercise to see how to code with Boost.Function you could try implementing the Pade(2,2) approximation to exp(z):

exp(z) ~ (12 + 6z + z^2)/(12 - 6z + z^2)

and compare accuracy.

 

[Daniel Duffy]

Just in case you missed the linlk, the full code for multiarray is here
http://www.datasimfinancial.com/forum/viewtopic.php?t=111

 

[Daniel Duffy]

Some more upcoming developments in boost are found on www.datasimfinancial.com and in particular the forthcoming book

http://www.datasimfinancial.com/forum/viewtopic.php?t=338

 

[Daniel Duffy]

I am prepariing some new blogs on boost, Now that we have discussed the fundamentals I will progress to boost Maths (Distributions, Special Functions, numerics) and bond and option pricing using exact and numerical methods. These will be in Volume II of the forthcoming Volume I book with Robert Demmming.

We also discuss how to integrate boost with C# using C++/CLI and then use these as worksheet functions in Excel. In this way we can reuse existing code without having to rewrite a complete library in C#.

 

[Larry]

Hi Daniel,
That link just shows .pdf files.

  Figure 5.1 Essential Libraries.pdf
  Figure 5.2 Supporting Libraries.pdf

and the author is:

  Cuchulainnso I must be missing something.Could you please provide another link to the code?TIA.-regardsLarry

 

[Larry]

I should have been more explicit.  I'm interested in code for
the ThreeFactorADESolver.