Thursday, 2 October 2014

Problem needed for research on Bermudan Option Pricing Algorithms


NAG together with Prof. Oosterlee and an MSc student from TU Delft are investigating the recent Stochastic Grid Bundling Method (SGBM) [1,2]. The objective is to compare the performance of SGBM to the well-known Longstaff-Schwartz (least squares method or LSM) in a non-academic setting, i.e. on the pricing of a Bermudan option, with underlying asset(s) driven by a realistic process such as Heston or LMM. We are looking for an interesting case to test these two methods. This includes the type of option, the underlying processes and any other important features or details.


The well known LSM by Longstaff-Schwartz[3] is the industry standard for pricing multi-dimensional Bermudan options by simulation and regression. LSM is based on the regression now principle, whereas the Stochastic Grid Bundling Method (SGBM) by Jain and Oosterlee applies regression later in order to get more accurate approximations. However, this limits us to apply SGBM to processes where an analytical or approximate expression of the discounted moments are available.

Another advantage of SGBM is its regression on bundles instead of the whole data set in a further attempt to decrease the regression error, and SGBM allows the computation of the Greeks at almost no extra cost.

Numerical results show that, compared to LSM, a higher accuracy can be obtained at comparable computational time. However, these tests were performed on academic problems, using geometric Brownian motion for the underlying assets.

Very recent research (Feng and Oosterlee, Sept 2014) extended SGBM for stochastic volatility and stochastic interest rate dynamics (Heston-Hull-White). This led to the possibility of comparing the performance of SGBM in a non-academic setting and is the cause for our research. We therefore wish to compare LSM and SGBM on a problem which is interesting and relevant to industry. Realistically, handling all the complexities of traded products will probably require more time than we have (around two months), so ideally we seek a problem which captures all the salient features and allowing us to see whether SGBM outperforms LSM.

Industry Involvement

We would like some advice in defining the product to be valued, especially regarding dimensionality, process, correlation structure, payoff and exercise features. Any additional industry involvement will be light: perhaps a few emails to clarify details, a conf call or face to face meeting.

For correspondence or further details please email


[1] S. Jain and C. W. Oosterlee, "The Stochastic Grid Bundling Method: Efficient Pricing of Bermudan Options and their Greeks,'' papers, SSRN, Sept. 2013.
[2] S. Jain and C. W. Oosterlee, "Pricing high-dimensional Bermudan options using the stochastic grid method,'' International Journal of Computer Mathematics, 89(9):1186-1211, 2012.
[3] F. Longstaff and E. Schwartz,  "Valuing American Options by Simulation: a Simple Least-squares Approach,''  Review of Financial Studies, vol. 14, no. 1, pp. 113-147, 2001.

Thursday, 21 August 2014

Gaussian Mixture Model

With the release of Mark 24 of the NAG C Library comes a plethora of new functionality including matrix functions, pricing Heston options w/term structure, best subset selection, and element-wise weightings for the nearest correlation matrix.

Among the new routines I was excited to test out was the Gaussian mixture model (g03ga). This routine will take a set of data points and fit a mixture of Gaussians for a given (co)variance structure by maximizing the log-likelihood function. The user inputs the (co)variance structure, number of groups, and (optionally) the initial membership probabilities.

I decided to test out this new functionality, which is also in Mark 24 of the NAG Toolbox for MATLAB. Often I will use MATLAB with the NAG Toolbox before switching to C++ and the NAG C Library for my production code. So I generated some data and tried the routine to see if it could find the covariance structure. You can download the script and try it out for yourself here. The example will generate the test data, run the NAG Gaussian mixture routine and plot the results. An example of the output is given below:

The blue points are the generated data, while the red and yellow ovals show the covariance structure output from NAG Gaussian mixture model (the ovals are contours of ~0.60 density for their respective groups).

While running the example a couple times and re-sampling through the starting values for the initial membership probabilities, I noticed what I thought to be unusual behavior for the routine. Namely, the Gaussian mixture model algorithm isn't able to identify the Gaussian mixtures. The function would occasionally converge to the below structure (run the above script and click 'Resample' 3 times):

It appears the routines has converged to local extrema of the likelihood function. This happens as a result of randomizing the initial membership probabilities.

Since we have the power of the NAG Library at our disposal, I've added a K-means clustering option in the above script to initialize the membership probabilities to a particular cluster before being input into the Gaussian mixture model.

My colleagues tell me that k-means can also get stuck in a local minima and exhibit this 'wrong' behavior as well, thus one should always be careful with initial allocations - luckily the NAG Library provides a generally acceptable default allocation as an option! Many thanks to Martyn Byng and Stephen Langdell for comments on this post.

Tuesday, 17 June 2014

Secrets of HPC Procurement

Liked my article today in HPC Wire "Secrets of the Supercomputers"? I firmly poke fun at various elements of an imaginary supercomputer procurement process. However, I'm sure many readers will also see familiar and painfully serious aspects in the fictional story.

As I mention at the bottom of that article, NAG can help make the process of buying HPC systems much better than the worrying case in the article.

For example, the tutorial I ran at SC13 with Terry Hewitt titled "Effective Procurement of HPC Systems" was very popular (~100 attendees). We have provided similar training as part of consulting engagements and we are now looking at running the tutorial again as an open short course.

Thursday, 10 April 2014

Testing Matrix Function Algorithms Using Identities

Edvin Deadman and Nick Higham (University of Manchester) write:

In a previous blog post we explained how testing new algorithms is difficult. We discussed the forward error (how far from the actual solution are we?) and the backward error (what problem have we actually solved?) and how we'd like the backward error to be close to the unit roundoff, u.

For matrix functions, we also mentioned the idea of using identities such as sin2A + cos2A = I to test algorithms. In practice, rather than I, we might find that we obtain a matrix R close to I, perhaps with ||R-I|| ≈ 10-13. What does this tell us about how the algorithms for sin A and cos A are performing? In particular, does it tell us anything about the backward errors? We've just written a paper which aims to answer these questions. This work is an output of NAG's Knowledge Transfer Partnership with the University of Manchester, so we thought we'd blog about it here.

Let's consider the identity exp(log A) - A = 0. Suppose that when we evaluate the left-hand side in floating point arithmetic we get a nonzero residual R rather than 0. We'll assume that this residual is  caused by some backward errors E1 and E2 so that exp(log(A + E1) + E2) = R. We'd like to investigate how big R can be when E1 and E2 are small, so we expand the left-hand side in a Taylor series to linear order. After a bit of algebra, the result is a linear operator relating R to E1 and E2R = L(E1E2). The operator is different for each identity considered, but it always involves the Fréchet derivatives of the matrix functions in the identity (the full gory details, including formulae for the linear operators associated with various identities, are in our paper).

We now have two options. The first option is to estimate the norm of the linear operator. This enables us to determine the maximum value of ||R|| consistent with backward errors of size u. The second option is to use the linear operator to explicitly estimate E1 and E2. This works just as well, but it is more expensive.

Our paper contains several examples demonstrating our approach in practice, so we'll just pick one here, involving the matrix exponential. The 10 x 10 Forsythe matrix is

where we set the parameter u =  2-53 ≈ 10-16 (the unit roundoff). It's known that the Schur-Parlett general-purpose matrix function algorithm struggles with this matrix, whereas algorithms based on scaling and squaring have no such problems. Does our approach reflect this behaviour? Using the identity ee-A = I, the maximum normwise residual consistent with backward stability is found to be 7.1 x 10-15. The Schur-Parlett approach gives a normwise residual of 7.7 x 10-11, but the scaling and squaring algorithm gives a residual of 2.2 x 10-15. Consistent with this, the backward error estimates are  2.5 x 10-11 and 4.5 x 10-16 respectively. So our approach has correctly identified the preferred algorithm.

There are some subleties which we haven't mentioned (for example certain matrices and identities can produce misleading cancellation effects) but we now have a useful method for testing matrix function algorithms when high-precision arithmetic isn't available. The relevant routines for computing matrix functions and Fréchet derivatives are available in chapter F01 of the NAG Library.

Friday, 21 March 2014

The Wilkinson Prize for Numerical Software

In honour of the outstanding contributions of James Hardy Wilkinson to the field of numerical software, Argonne National Laboratory, the National Physical Laboratory, and the Numerical Algorithms Group award the Wilkinson Prize for Numerical Software (US $3000).

The 2015 prize will be awarded at the International Conference in Industrial and Applied Mathematics (ICIAM) in Beijing, China, August 2015. Entries must be received by July 1, 2014. Additional details on the Wilkinson Prize for Numerical Software and the official rules can
be found at the URL:

Submissions can be sent by email to, contact this address for further information.

Previous prizes have also been awarded at ICIAM:

1991 - Linda Petzold for DASSL
1995 - Chris Bischof and Alan Carle for ADIFOR 2.0
1999 - Matteo Frigo and Steven Johnson for FFTW
2003 - Jonathan Shewchuk  for Triangle
2007 - Wolfgang Bangerth, Ralf Hartmann and Guido Kanschat for deal.II
2011 - Andreas Waechter and Carl D. Laird for Ipopt

Mike Dewar, Maurice Cox, Jorge Moré
Board of Trustees

Thursday, 6 February 2014

C++ wrappers for the NAG C Library


Occasionally, we receive requests to make the NAG C Library easier to call from C++. In the past, we found it difficult to build something that would work across all of the code our C++ users write. With the advent of the C++11 standard, many of the key features of the widely used Boost library have been incorporated into the STL, and finally provide a standardized way to address many of the difficulties we've encountered (the code we describe here works with Visual Studios 2010 and later, as well as several different versions of the Intel compiler and gcc).

We have created example wrappers that can serve as templates for creating C++ wrappers around NAG functions. Specifically, the examples now allow the user to:

1) Pass function pointers, functors, class member functions and lamda functions as callbacks to the NAG library.
2) Use raw pointers, smart pointers, STL containers or boost containers to store data and pass these to the NAG library.

 A note: these are NOT a C++ interface, but merely wrappers around the C Library.

Inside the Wrappers

Let's take a look at the protoype for one function. Here is the new signature for a NAG minimization routine:

void e04abcpp(const std::function<void(double,double*,NagComm*)> &callback, double e1, double e2, double *a, double *b, Integer max_fun, double *x, double *f, NagComm *user, NagError *error);

There are two differences from the NAG C routine; the callback function and the NagComm class (which we will discuss in a bit). Using the above signature, the NAG Library can be called via a function, a pointer to a function, a function object with a common interface, or a class member function. For example, the user could have the following NAG calls:

e04abcpp(MyFunctor(), e1, e2, &a, &b, max_fun, &x, &f, &comm, NULL);
e04abcpp([&](double x,double*fc,NagComm*comm) {
            return myclass.myFunc(x,fc,comm);
        }, e1, e2, &a, &b, max_fun, &x, &f, &comm, NULL);

Here, the MyFunctor and myclass.myFunc are the callbacks that return values of our objective function. In order to pass these to the C Library, we need a structure to evaluate and store them. The natural way to do this is to extend the NAG Communication class. Since we dont want a new NAG Comm class affecting the underlying Library, we create the first member to be the old comm structure (a pointer to a structure can be treated as as a pointer to the first member of the structure). We can then place our functional and new callback inside this structure:

class NagComm  // New Comm Class
    Nag_Comm comm;   // Old Comm Class
    friend void e04abcpp(const std::function<void(double,double*,NagComm*)> &callback, double e1, double e2, double *a, double *b, Integer max_fun, double *x, double *f, NagComm *user, NagError *error);

    void * puser;
    // A pointer to the C++11 std::function object
    const std::function<void(double,double*,NagComm*)> * callbackCPP11;

    static void NAG_CALL e04abDelegateCPP11(double xc, double *fc, Nag_Comm *comm) {
        // Safe cast since access to this routine is controlled
        NagComm * user = (NagComm*)comm;
        // Call the function polymorphically
        (* (user->callbackCPP11) )(xc, fc, user);

    NagComm() {
        puser = NULL;
        callbackCPP11 = NULL;

Once we have the new Comm class, we can make the cpp wrapper function:

void e04abcpp(const std::function<void(double,double*,NagComm*)> &callback, double e1, double e2, double *a, double *b, Integer max_fun, double *x, double *f, NagComm *user, NagError *error);
    // This is the only place where NagUser::callbackCPP11 is assigned
    user->callbackCPP11 = &callback;
    e04abc(NagComm::e04abDelegateCPP11, e1, e2, a, b, max_fun, x, f, (Nag_Comm*)user, error);

Smart Pointers/Containers

Want to pass different sets of pointers and containers into NAG routines? These are easily handled as well, for example in the follow wrapper to g01amc:

template<typename RV, typename Q, typename QV>
void g01amcpp(Integer n, RV && rv, Integer nq, Q && q, QV && qv, NagError *error)

Users can call the function using STL containers, boost types or smart pointers:

double rv[n]={.5, .729, .861, .44, .791, .001, .062, .912, .27, .141, .32, .133, .654,
     .285, .553, .438, .316, .696, .718, .293, .704, .029};
std::vector<double> q={0.0, .25, .73, .9, 1.0};
std::unique_ptr<double[]> qv(new double[nq]);
g01amcpp(n, rv, nq, q, qv, NULL);  

The g01amcpp wrapper is defined as follows:

template<typename RV, typename Q, typename QV>
void g01amcpp(Integer n, RV && rv, Integer nq, Q && q, QV && qv, NagError *error)
    // Get primitive pointers from input types.     
    double * p_rv = nagprivate::getPtr(rv);
    const double * p_q = nagprivate::getPtr(q);
    double * p_qv = nagprivate::getPtr(qv);

    g01amc(n, p_rv, nq, p_q, p_qv, error);

A different template parameter for each array argument means different container/smart pointer classes can be used for each argument. Inside the wrapper the function nagprivate::getPtr is called to return the address of the first element in each array. This function uses the auto/decltype feature which has a new meaning in C++11:

namespace nagprivate {
template<typename T> auto getPtr(T & p) -> decltype( &p[0] )
         return &(p[0]);

Note: the wrapper above assumes that the pointer/container classes implement operator[] and use contiguous arrays for internal storage. The wrapper will most likely not work for things such as linked lists.

Code Availability

We have created some example templates and more will available over time so key an eye on the NAG and C website and the bottom of this blog for updates. The examples are tested and work with VS10, g++ 4.7, and several versions of the Intel compilers. The wrappers are available in source form and can be compiled at the command line by:

$ cl -EHsc /MT -I"C:\Program Files\NAG\CL23\clw6i23dal"\include e04ab.cpp /link /LIBPATH:"C:\Program Files\NAG\CL23\clw6i23dal"\lib "C:\Program Files\NAG\CL23\clw6i23dal"\lib\nagc_nag_MT.lib user32.lib

 $ g++ -std=c++11 e04ab.cpp -I/opt/NAG/cll6a23dgl/include/ /opt/NAG/cll6a23dgl/lib/libnagc_nag.a -lpthread -lm -o e04ab.exe


We would be interested in hearing feedback from users. You can leave a comment below or email with the subject 'NAG and C++ wrappers'.

March 14 Update

We now have 14 examples including optimization, matrix functions, nearest correlation matrices, and quantile analysis that are available in the above links.

Thursday, 16 January 2014

Out and about this week – The London Thalesians Seminar

NAG’s Brian Spector gave a great talk to a packed audience of finance professionals in London this week. The Thalesians describe themselves as a “think tank of dedicated professionals with an interest in quantitative finance, economics, mathematics, physics and computer science”. Brian was delighted to present "Implied Volatility using Python's Pandas Library" at their recent London Seminar on Wednesday 15 January 2014.

Brian Spector presenting "Implied Volatility using Python's Pandas Library at the recent London Thalesians Seminar

You can learn more about the subject of Brian’s talk in one of his NAG and Python blogs below:

Additional NAG and Python information features on our website including how you can download the NAG Python Bindings that will enable use of the NAG C Library from Python.
You can follow The Thalesians on Twitter @thalesians and NAG @NAGTalk.