QuantLib : Least squares method implementation

A couple of years ago I published a blog post on using Microsoft Solver Foundation for curve fitting. In this post we go through that example again, but this time using corresponding Quantlib (QL) tools. In order to make everything a bit more concrete, the fitting scheme (as solved by Frontline Solver in Microsoft Excel) is presented in the screenshot below.

Second degree polynomial fitting has been performed, in order to find a set of coefficients which are minimizing sum of squared differences between observed rates and estimated rates. The same minimization scheme has also been used in my example program presented later below.

LeastSquares class

After some initial woodshedding with QL optimization tools, I started to dream if it could be possible to create kind of a generic class, which could then handle all kinds of different curve fitting schemes.

When using LeastSquares class, client is creating its instance by setting desired optimization method (Quantlib::OptimizationMethod) and all parameters needed for defining desired end criteria (Quantlib::EndCriteria) for optimization procedure. The actual fitting procedure will be started by calling Fit method of LeastSquares class. For this method, client is giving independent values (Quantlib::Array), vector of function pointers (boost::function), initial parameters (Quantlib::Array) and parameter constraints (Quantlib::Constraint).

Objective function implementation (Quantlib::CostFunction) used by QL optimizers is hidden as a nested class inside LeastSquares class. Since the actual task performed by least squares method is always to minimize a sum of squared errors, this class is not expected to change. Objective function is then using a set of given independent values and function pointers (boost::function) to calculate sum of squared errors between independent (observed) values and estimated values. By using function pointers (boost::function), algorithm will not be hard-coded into class method.

Finally, for function pointers used in objective function inside LeastSquares class, implementation class with correct method signature is expected to be available. In the example program, there is CurvePoint class, which is modelling a rate approximation for a given t (maturity) and a set of given external parameters.

Example program

Create a new C++ console project and copyPaste the following into corresponding header and implementation files. Help is available here if needed. In documentations page, there are three excellent presentation slides on Boost and Quantlib libraries by Dimitri Reiswich. "Hello world" examples of QuantLib optimizers can be found within these files.

Finally, thanks again for reading my blog.
-Mike

// CurvePoint.h
#pragma once
#include <ql/quantlib.hpp>
using namespace QuantLib;
//
// class modeling rate approximation for a single 
// curve point using a set of given parameters
class CurvePoint
{
private:
 Real t;
public:
 CurvePoint(Real t);
 Real operator()(const Array& coefficients);
};
//
//
//
// CurvePoint.cpp
#include "CurvePoint.h"
//
CurvePoint::CurvePoint(Real t) : t(t) { }
Real CurvePoint::operator()(const Array& coefficients)
{
 Real approximation = 0.0;
 for (unsigned int i = 0; i < coefficients.size(); i++)
 {
  approximation += coefficients[i] * pow(t, i);
 }
 return approximation;
}
//
//
//
// LeastSquares.h
#pragma once
#include <ql/quantlib.hpp>
using namespace QuantLib;
//
class LeastSquares
{
public:
 LeastSquares(boost::shared_ptr<OptimizationMethod> solver, Size maxIterations, Size maxStationaryStateIterations, 
  Real rootEpsilon, Real functionEpsilon, Real gradientNormEpsilon);
 //
 Array Fit(const Array& independentValues, std::vector<boost::function<Real(const Array&)>> dependentFunctions, 
  const Array& initialParameters, Constraint parametersConstraint);
private:
 boost::shared_ptr<OptimizationMethod> solver;
 Size maxIterations; 
 Size maxStationaryStateIterations;
 Real rootEpsilon; 
 Real functionEpsilon;
 Real gradientNormEpsilon;
 //
 // cost function implementation as private nested class 
 class ObjectiveFunction : public CostFunction
 {
 private:
  std::vector<boost::function<Real(const Array&)>> dependentFunctions;
  Array independentValues;
 public:
  ObjectiveFunction(std::vector<boost::function<Real(const Array&)>> dependentFunctions, const Array& independentValues);
  Real value(const Array& x) const;
  Disposable<Array> values(const Array& x) const;
 };
};
//
//
//
// LeastSquares.cpp
#pragma once
#include "LeastSquares.h"
//
LeastSquares::LeastSquares(boost::shared_ptr<OptimizationMethod> solver, Size maxIterations, 
 Size maxStationaryStateIterations, Real rootEpsilon, Real functionEpsilon, Real gradientNormEpsilon)
 : solver(solver), maxIterations(maxIterations), maxStationaryStateIterations(maxStationaryStateIterations),
 rootEpsilon(rootEpsilon), functionEpsilon(functionEpsilon), gradientNormEpsilon(gradientNormEpsilon) { }
//
Array LeastSquares::Fit(const Array& independentValues, std::vector<boost::function<Real(const Array&)>> dependentFunctions, 
 const Array& initialParameters, Constraint parametersConstraint)
{
 ObjectiveFunction objectiveFunction(dependentFunctions, independentValues);
 EndCriteria endCriteria(maxIterations, maxStationaryStateIterations, rootEpsilon, functionEpsilon, gradientNormEpsilon);
 Problem optimizationProblem(objectiveFunction, parametersConstraint, initialParameters);
 //LevenbergMarquardt solver;
 EndCriteria::Type solution = solver->minimize(optimizationProblem, endCriteria);
 return optimizationProblem.currentValue();
}
LeastSquares::ObjectiveFunction::ObjectiveFunction(std::vector<boost::function<Real(const Array&)>> dependentFunctions, 
 const Array& independentValues) : dependentFunctions(dependentFunctions), independentValues(independentValues) { }
//
Real LeastSquares::ObjectiveFunction::value(const Array& x) const
{
 // calculate squared sum of differences between
 // observed and estimated values
 Array differences = values(x);
 Real sumOfSquaredDifferences = 0.0;
 for (unsigned int i = 0; i < differences.size(); i++)
 {
  sumOfSquaredDifferences += differences[i] * differences[i];
 }
 return sumOfSquaredDifferences;
}
Disposable<Array> LeastSquares::ObjectiveFunction::values(const Array& x) const
{
 // calculate differences between all observed and estimated values using 
 // function pointers to calculate estimated value using a set of given parameters
 Array differences(dependentFunctions.size());
 for (unsigned int i = 0; i < dependentFunctions.size(); i++)
 {
  differences[i] = dependentFunctions[i](x) - independentValues[i];
 }
 return differences;
}
//
//
//
// tester.cpp
#include "LeastSquares.h"
#include "CurvePoint.h"
//
int main()
{
 // 2nd degree polynomial least squares fitting for a curve
 // create observed market rates for a curve
 Array independentValues(9);
 independentValues[0] = 0.19; independentValues[1] = 0.27; independentValues[2] = 0.29; 
 independentValues[3] = 0.35; independentValues[4] = 0.51; independentValues[5] = 1.38; 
 independentValues[6] = 2.46; independentValues[7] = 3.17; independentValues[8] = 3.32;
 //
 // create corresponding curve points to be approximated
 std::vector<boost::shared_ptr<CurvePoint>> curvePoints;
 curvePoints.push_back(boost::shared_ptr<CurvePoint>(new CurvePoint(0.083)));
 curvePoints.push_back(boost::shared_ptr<CurvePoint>(new CurvePoint(0.25)));
 curvePoints.push_back(boost::shared_ptr<CurvePoint>(new CurvePoint(0.5)));
 curvePoints.push_back(boost::shared_ptr<CurvePoint>(new CurvePoint(1.0)));
 curvePoints.push_back(boost::shared_ptr<CurvePoint>(new CurvePoint(2.0)));
 curvePoints.push_back(boost::shared_ptr<CurvePoint>(new CurvePoint(5.0)));
 curvePoints.push_back(boost::shared_ptr<CurvePoint>(new CurvePoint(10.0)));
 curvePoints.push_back(boost::shared_ptr<CurvePoint>(new CurvePoint(20.0)));
 curvePoints.push_back(boost::shared_ptr<CurvePoint>(new CurvePoint(30.0)));
 //
 // create container for function pointers for calculating rate approximations
 std::vector<boost::function<Real(const Array&)>> dependentFunctions;
 //
 // for each curve point object, bind function pointer to operator() and add it into container
 for (unsigned int i = 0; i < curvePoints.size(); i++)
 {
  dependentFunctions.push_back(boost::bind(&CurvePoint::operator(), curvePoints[i], _1));
 }
 // perform least squares fitting and print optimized coefficients
 LeastSquares leastSquares(boost::shared_ptr<OptimizationMethod>(new LevenbergMarquardt), 10000, 1000, 1E-09, 1E-09, 1E-09);
 NoConstraint parametersConstraint;
 Array initialParameters(3, 0.0);
 Array coefficients = leastSquares.Fit(independentValues, dependentFunctions, initialParameters, parametersConstraint);
 for (unsigned int i = 0; i < coefficients.size(); i++)
 {
  std::cout << coefficients[i] << std::endl;
 }
 return 0;
}