---

CurveFitting.cpp

// C++ includes

#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>

using std::cout;
using std::cerr;
using std::endl;
using std::setw;

// personal headers
#include <ql/date.hpp>
#include <ql/array.hpp>
#include <ql/Math/matrix.hpp>
#include <ql/errors.hpp>

#include <portability.h>
#include <leastsquare.h>
#include <timer.h>

using QuantLib::Date;
using QuantLib::Array;
using QuantLib::Math::Matrix;
using QuantLib::Math::transpose;
using QuantLib::Error;



inline
std::ostream& operator<< (std::ostream& os, const QuantLib::Array& v)
{
  size_t i, sz = v.size();
 
  os << "[";
  for (i=0; i<sz-1; ++i)
    os << v[i] << " ";
  os << v[sz-1] 
     << "]" 
     << std::endl;
  
  return os;
}

class BraceVolatilityFunction {
  double a_, b_, c_, d_;
public :
  inline BraceVolatilityFunction(double a, double b, double c, double d)
    : a_(a), b_(b), c_(c), d_(d) {}
  inline ~BraceVolatilityFunction() {}

  inline double operator() (double t) 
  { return (a_+b_*t)*exp(-c_*t) + d_; }

  inline double da(double t) 
  { return exp(-c_*t); }

  inline double db(double t) 
  { return t*exp(-c_*t); }

  inline double dc(double t) 
  { return -t*(a_+b_*t)*exp(-c_*t); }

  inline double dd(double t) 
  { return 1.; }
};


class CurveFittingProblem 
  : public LeastSquareProblem<Array,Matrix > {
  const Array &ttm_;
  const Array &target_;
public :
  CurveFittingProblem(const Array& ttm, const Array& target)
    : ttm_(ttm), target_(target) {}
  virtual ~CurveFittingProblem() {}

  virtual int size() { return ttm_.size(); }

  virtual void targetAndValue(const Array& abcd, 
                  Array& target, 
                  Array& fct2fit)
  {
    BraceVolatilityFunction bvf(abcd[0],abcd[1],abcd[2],abcd[3]);
    target = target_;// target values
    for (int i=0; i<ttm_.size(); ++i)
      fct2fit[i] = bvf(ttm_[i]);
  }

  virtual void targetValueAndfirstDerivative(const Array& abcd, 
                         Matrix& grad_fct2fit, 
                         Array& target, 
                         Array& fct2fit)
  {
    BraceVolatilityFunction bvf(abcd[0],abcd[1],abcd[2],abcd[3]);
    target = target_;// target values
    for (int i=0; i<ttm_.size(); ++i) {
      // function value at current point abcd
      fct2fit[i] = bvf(ttm_[i]);
      /*
    matrix of first derivatives :
    the derivatives with respect to the parameter a,b,c,d
    are stored by row.
      */
    grad_fct2fit[i][0] = bvf.da(ttm_[i]);
    grad_fct2fit[i][1] = bvf.db(ttm_[i]);
    grad_fct2fit[i][2] = bvf.dc(ttm_[i]);
    grad_fct2fit[i][3] = bvf.dd(ttm_[i]);
    }

  }
  };

/*
  We define here an inverse problem to show how to fit
  parametric function to data. 
*/
int main()
{
  /*
    Parameter values that produce the volatility hump.
    Consider it as optimal values of the curve fitting
    problem.
  */
  double abcd_[4] = {0.147014, 0.057302, 0.249964, 0.148556 };
  Array abcd(4);
  abcd[0] = abcd_[0];
  abcd[1] = abcd_[1];
  abcd[2] = abcd_[2];
  abcd[3] = abcd_[3];

  cout << "Optimal values  : " << abcd << endl;

  // Define the target volatility function
  BraceVolatilityFunction bvf(abcd[0],abcd[1],abcd[2],abcd[3]);

  Timer t;
  t.start();
  try {
 
    double t, 
      startDate = 0.0, // start date of volatilty
      endDate = 20., // end date of volatility
      period = 0.25; // period length between values (in year fraction : quarterly)

    int i, periodNumber = (int)(endDate / period); // number of period

    Array targetValue(periodNumber), timeToMaturity(periodNumber);
    // Fill target and time to maturity arrays
    for (i=0; i<periodNumber; ++i) {
      t = startDate + i*period;
      timeToMaturity[i] = t;
      targetValue[i] = bvf(t);
    }

    // Accuracy of the optimization method
    double accuracy = 1e-10;// It is the square of the accuracy
    // Maximum number of iterations
    int maxiter = 10000;

    Array initialValue(4, 0.1);
        
    // Least square optimizer
    NonLinearLeastSquare<Array > lsqnonlin(accuracy,maxiter);

    // Define the least square problem
    CurveFittingProblem cfp(timeToMaturity, targetValue);
     
    // Set initial values
    lsqnonlin.setInitialValue(initialValue);
    // perform fitting
    Array solution = lsqnonlin.Perform(cfp);
        
    cout << endl 
     << "Optimal values  : " << abcd
     << "Solution values : " << solution
     << "Solution Error  : " << (solution - abcd ) << endl;

    // save in given format
    BraceVolatilityFunction bvf2(solution[0],solution[1],solution[2],solution[3]);

    // Gnuplot format
    /*
      std::ofstream file("curve.dta");
      for (i=0; i<periodNumber; ++i) {
      t = startDate + i*period;
      file << setw(8) << t << setw(20) << bvf(t) << setw(20) << bvf2(t) << endl;
      }
    */
        
    std::ofstream file("curve.csv");
    const char separator[] = ", ";
    // CSV format
    for (i=0; i<periodNumber; ++i) {
      t = startDate + i*period;
      file << t << separator << bvf(t) << separator << bvf2(t) << endl;
    }
    file.close();

  }
 catch(Error & e)
   {
     cerr << e.what() << endl;;
   }
 t.stop();

 cout << "Ellapsed time : " << (double) t << " seconds." << endl;


 return 1;
}

---