Python-SciPy: Optimizing Smooth Libor Forward Curve

In my previous post, I presented how to use Python SciPy Optimization package for solving zero-coupon rate term structure from a given set of zero-coupon bond prices numerically. In this post, the same approach will be used in order to solve smooth Libor forward curve from a given set of vanilla swaps. Example of applying this scheme using Solver and Excel/VBA can be found in here.

In order to value fixed income derivatives cash flows, relevant forward rates and discount factors have to be defined from bootstrapped zero-coupon curve. In a Libor world, we use cash and FRA contracts (or futures contracts) in a short-end of the curve, while in a long-end of the curve we use swap contracts. Let us assume for a moment, that we bootstrap USD 3M zero-coupon curve, in order to get 3M forward rates on a quarterly basis. While bootstrapping approach usually gives smooth curve within short-end of the curve, it will usually fail to do this within long-end of the curve. This happens, because there is no swap contracts available in a long-end of the curve and hereby, we need to do a lot of interpolations for required swap rates. The resulting forward curve then looks like a chain of waterfalls and for the long-end of the curve, smoothness will be lost. This topic has been thoroughly covered in chapter three within this excellent book by Richard Flavell.

In a nutshell, we first define the shortest end of the curve (using spot cash rate) and then we just guess the rest of the forward rates. After this, we minimize sum of squared differences between adjacent forward rates, subject to constraints. In this case, constraints are present values of swaps, which by definition will be zero. When this optimization task has been performed, the resulting forward curve will successfully re-price the current swap market.

Bloomberg forward curve screen for the same set of data is available in here. Finally, thanks again for reading my blog.
-Mike

import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as pl

# sum of squared errors of all decision variables
# args: 0 = array of given initial rates, 1 = scaling factor
def ObjectiveFunction(x, args):
    # concatenate given initial rates and given 'guesses' for forward rates
    x = np.concatenate([args[0], x])
    return np.sum(np.power(np.diff(x), 2) * args[1])

# x = array of Libor forward rates
# args: 0 = swap rate, 1 = years to maturity, 2 = floating leg payments per year, 
# 3 = notional, 4 = array of given initial rates
def VanillaSwapPV(x, args):
    # PV (fixed payer swap) = -swapRate * Q(T) + (1 - DF(T))
    # where Q(T) is sumproduct of fixed leg discount factors and corresponding time accrual factors
    # assumption : fixed leg is always paid annually
    # since fixed leg is paid annually and maturity is always integer, 
    # time accrual factor will always be exactly 1.0
    
    # concatenate given initial rates and given 'guesses' for forward rates
    x = np.concatenate([args[4], x])
    DF = 0.0
    Q  = 0.0
    floatingLegTenor = 1 / args[2]
    nextFixedLegCouponDate = 1
    currentTimePoint = 0.0
    nCashFlows = int(args[1] * args[2])
    
    for i in range(nCashFlows):
        currentTimePoint += floatingLegTenor
        # first rate is always spot rate: calculate first spot discount factor
        if (i == 0):
            DF = 1 / (1 + x[i] * floatingLegTenor)
        # other rates are always forward rates
        if (i > 0):
            # solve current spot discount factor using previous spot discount factor
            # and current forward discount factor: DF(0,T) = DF(0,t) * DF(t,T), where (0 < t < T)
            DF *= (1 / (1 + x[i] * floatingLegTenor))
            if ((currentTimePoint - nextFixedLegCouponDate) == 0):
                Q += DF
                nextFixedLegCouponDate += 1
                
    return (-args[0] * Q + (1.0 - DF)) * args[3]

# array of given initial rates (the first given rate has to be spot cash rate)
# here, we give only one: 3M spot Libor fixing rate ('cash-to-first-future')
initialRates = np.array([0.006276])
scalingFactor = 1000000.0

# vanilla swap pricing functions as constraints
# args: swapRate, yearsToMaturity, floatingLegPaymentFrequency, notional, array of given initial rates
swaps = ({'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.0078775, 1, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.009295, 2, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.0103975, 3, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.01136, 4, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.012268, 5, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.013183, 6, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.01404, 7, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.014829, 8, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.01554, 9, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.016197, 10, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.016815, 11, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.017367, 12, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.018645, 15, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.01996, 20, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.020615, 25, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.02097, 30, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.021155, 40, 4, scalingFactor, initialRates]]},
        {'type': 'eq', 'fun': VanillaSwapPV, 'args': [[0.021017, 50, 4, scalingFactor, initialRates]]})

# initial guesses for Libor forward rates up to 50 years
# number of rates: 50 * 4 - 1 (years to maturity * floating leg payments per year - first given cash spot rate)
initialGuesses = np.full(199, 0.006)
model = opt.minimize(ObjectiveFunction, initialGuesses, args = ([initialRates, scalingFactor]), method = 'SLSQP', options = {'maxiter': 500}, constraints = swaps)

# print selected model results
print('Success: ' + str(model.success))
print('Message: ' + str(model.message))
print('Number of iterations: ' + str(model.nit))
print('Objective function (sum of squared errors): ' + str(model.fun))
#print('Changing variables (Libor forward rates): ' + str(model.x * 100))
pl.plot(model.x)
pl.show()

Alglib : Ho-Lee calibration in C#

A couple of years ago I published post on Ho-Lee short interest rate model calibration using Microsoft Solver Foundation (MSF). Solver is still available (somewhere) but not supported or developed any further, since Microsoft terminated that part of their product development (if I have understood the story correctly). As a personal note I have to say, that MSF was actually not the easiest package to use, from programming point of view.

Life goes on and very fortunately there are a lot of other alternatives available. This time I wanted to present Alglib, which is offering quite impressive package of numerical stuff already in their non-commercial (free of charge) edition. In optimization category, there is Levenberg-Marquardt algorithm (LMA) implementation available, which can be used for multivariate optimization tasks. Within this post, I am re-publishing short interest rate model calibration scheme, but only using Alglib LMA for performing the actual optimization routine.

In the case that Alglib is completely new package, it is recommended to check out some library documentation first, since it contains a lot of simple "Copy-Paste-Run"-type of examples, which will help to understand how the flow of information is processed. As a personal note I have to say, that time and effort needed to "getting it up and running" is extremely low. Just find Alglib download page, download proper C# version, create a new C# project and add reference to DLL file (alglibnet2.dll) included in download folder.

The outcome

Market prices of zero-coupon bonds, volatility (assumed to be estimated constant) and short rate are used, in order to solve time-dependent theta coefficients for Ho-Lee short interest rate model. When these coefficients have been solved, the model can be used to price different types of interest rate products. The original data and other necessary details have been already presented in here.

Alglib LMA package offers a lot of flexibility, concerning how the actual optimization task will be solved. Basically, there are three different schemes : V (using function vector only), VJ (using function vector and first order partial derivatives, known as Jacobian matrix) and FGH (using function vector, gradient and second order partial derivatives, known as Hessian matrix). Solved time-dependent theta coefficients for the first two schemes (V, VJ) are presented in the screenshot below.

Due to the usage of analytical first order partial derivatives instead of numerical equivalents (finite difference), Vector-Jacobian scheme (VJ) is a bit more accurate and computationally much faster. For calculating required partial derivatives for Jacobian/Hessian matrix, one might find this online tool quite handy.

The program

In order to test this calibration program, just create a new C# project (AlgLibTester), copy-paste the following code into a new CS file and add reference to Alglib DLL file.

A few words on this program. First of all, Alglib LMA has been wrapped inside static LevenbergMarquardtSolver class, in order to provide kind of a generic way to use that solver. Basically, LMA requires callback method(s) for calculating values for vector of functions and/or Jacobian matrix. In this program, these callback methods are first defined as interfaces, which (the both in this program) are going to be implemented in HoLeeZeroCouponCalibration class. This class is also storing all the "source data" (maturities, zero-coupon bond prices, volatility and short rate) to be used, as LMA is processing its data (coefficients to be solved) through these callback methods. Note, that LevenbergMarquardtSolver class has absolutely no information about external world, except that it can receive callback methods defined in IFunctionVector and IJacobianMatrix interfaces. With this scheme, we can separate all financial algorithms and data from the actual optimization task.

using System;
using System.Linq;
//
// 'type definitions' for making life a bit more easier with long class names
using LM = AlgLibTester.LevenbergMarquardtSolver;
using HoLee = AlgLibTester.HoLeeZeroCouponCalibration;

namespace AlgLibTester
{
    class Program
    {
        static void Main(string[] args)
        {
            // short rate volatility, short rate, vector of maturities, vector of zero-coupon bond prices
            double sigma = 0.00039;
            double r0 = 0.00154;
            double[] t = new double[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 };
            double[] zero = new double[] { 0.9964, 0.9838, 0.9611, 0.9344, 0.9059, 0.8769, 0.8478, 0.8189, 0.7905, 0.7626 };
            //
            // create callback functions wrapped inside a class, which 
            // implements IFunctionVector and IJacobianMatrix interfaces
            HoLee callback = new HoLee(sigma, r0, t, zero);
            //
            // container for initial guesses for theta coefficients, which are going to be solved
            double[] theta = null;
            //
            // Example 1 :
            // use function vector only (using 'V' mode of the Levenberg-Marquardt optimizer)
            theta = new double[] { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
            LM.Solve(ref theta, 10, (IFunctionVector)callback);
            Console.WriteLine("Using function vectors only");
            Console.WriteLine("iterations : {0}", LM.report.iterationscount);
            theta.ToList<double>().ForEach(it => Console.WriteLine(it));
            Console.WriteLine("");
            //
            // Example 2 :
            // use function vector and jacobian matrix (using 'VJ' mode of the Levenberg-Marquardt optimizer)
            theta = new double[] { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }; 
            LM.Solve(ref theta, 10, (IFunctionVector)callback, (IJacobianMatrix)callback);
            Console.WriteLine("Using function vectors and Jacobian matrix");
            Console.WriteLine("iterations : {0}", LM.report.iterationscount);
            theta.ToList<double>().ForEach(it => Console.WriteLine(it));
        }
    }
    //
    //
    //
    // static class, which is wrapper for Alglib Levenberg-Marquardt solver
    // in order to make the usage of this particular solver a bit more 'generic'
    public static class LevenbergMarquardtSolver
    {
        // alglib report is exposed to public
        public static alglib.minlmreport report;
        private static alglib.minlmstate state;
        //
        // create minimization solver using only function vector ('V' vector method)
        public static void Solve(ref double[] x, int m, IFunctionVector functionVector, 
            // optional parameters which are set to default values
            double diffstep = 1E-08, double epsg = 1E-15, double epsf = 0.0, double epsx = 0.0, int maxits = 0)
        {
            // create LM model, set termination conditions, perform optimization and get results
            alglib.minlmcreatev(m, x, diffstep, out state);
            alglib.minlmsetcond(state, epsg, epsf, epsx, maxits);
            alglib.minlmoptimize(state, functionVector.callback, null, null);
            alglib.minlmresults(state, out x, out report);
        }
        // method overloading : create minimization solver using function vector and Jacobian matrix ('VJ' vector-jacobian method)
        public static void Solve(ref double[] x, int m, IFunctionVector functionVector, IJacobianMatrix jacobianMatrix,
            // optional parameters which are set to default values
            double epsg = 1E-15, double epsf = 0.0, double epsx = 0.0, int maxits = 0)
        {
            // create LM model, set termination conditions, perform optimization and get results
            alglib.minlmcreatevj(m, x, out state);
            alglib.minlmsetcond(state, epsg, epsf, epsx, maxits);
            alglib.minlmoptimize(state, functionVector.callback, jacobianMatrix.callback, null, null);
            alglib.minlmresults(state, out x, out report);
        }
    }
    //
    //
    //
    // interface : definition for function vector callback method, required by LevenbergMarquardtSolver class
    public interface IFunctionVector
    {
        void callback(double[] x, double[] fi, object obj);
    }
    //
    //
    //
    // interface : definition for jacobian matrix callback method, optionally required by LevenbergMarquardtSolver class
    public interface IJacobianMatrix
    {
        void callback(double[] x, double[] fi, double[,] jac, object obj);
    }
    //
    //
    //
    // Ho-Lee calibration class, which implements IFunctionVector and IJacobianMatrix interfaces
    public class HoLeeZeroCouponCalibration : IFunctionVector, IJacobianMatrix
    {
        // parameters, which are required in order to calculate
        // zero-coupon bond prices using Ho-Lee short rate model
        private double sigma;
        private double r0;
        private double[] t;
        private double[] zero;
        //
        public HoLeeZeroCouponCalibration(double sigma, double r0, double[] t, double[] zero)
        {
            this.sigma = sigma;
            this.r0 = r0;
            this.t = t;
            this.zero = zero;
        }
        // callback method for calculating vector of values for functions
        void IFunctionVector.callback(double[] x, double[] fi, object obj)
        {
            // calculate squared differences of Ho-Lee prices (using theta coefficients) 
            // and market prices and then assign these differences into function vector fi
            for (int i = 0; i < fi.Count(); i++)
            {
                double HOLEE_ZERO = Math.Exp(-Math.Pow(t[i], 2.0) * x[i] / 2.0
                    + Math.Pow(sigma, 2.0) * Math.Pow(t[i], 3.0) / 6.0 - r0 * t[i]);
                fi[i] = Math.Pow(zero[i] - HOLEE_ZERO, 2.0);
            }
        }
        // callback method for calculating partial derivatives for jacobian matrix
        void IJacobianMatrix.callback(double[] x, double[] fi, double[,] jac, object obj)
        {
            double HOLEE_ZERO = 0.0;
            // part 1 : calculate squared differences of Ho-Lee prices (using theta coefficients) 
            // and market prices, then assign these squared differences into function vector fi
            for (int i = 0; i < fi.Count(); i++)
            {
                HOLEE_ZERO = Math.Exp(-Math.Pow(t[i], 2.0) * x[i] / 2.0
                    + Math.Pow(sigma, 2.0) * Math.Pow(t[i], 3.0) / 6.0 - r0 * t[i]);
                fi[i] = Math.Pow(zero[i] - HOLEE_ZERO, 2.0);
            }           
            // part 2 : calculate all partial derivatives for Jacobian square matrix
            // loop through m functions
            for (int i = 0; i < fi.Count(); i++)
            {
                // loop through n theta coefficients
                for (int j = 0; j < x.Count(); j++)
                {
                    double derivative = 0.0;
                    // partial derivative is non-zero only for diagonal cases
                    if (i == j)
                    {
                        HOLEE_ZERO = Math.Exp(-Math.Pow(t[j], 2.0) * x[j] / 2.0
                            + Math.Pow(sigma, 2.0) * Math.Pow(t[j], 3.0) / 6.0 - r0 * t[j]);
                        derivative = Math.Pow(t[j], 2.0) * (zero[j] - HOLEE_ZERO) * HOLEE_ZERO;
                    }
                    jac[i, j] = derivative;
                }
            }
        }
    }
}

Finally, thanks for reading my blog.
-Mike

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