Showing posts with label Boost library. Show all posts
Showing posts with label Boost library. Show all posts

Sunday, April 22, 2018

C++/CLI Interoperability : Using QuantLib in C#

Assume the following scenario : there is QuantLib program available, written in native C++ and you would like to use it from your C# program. What would you do? Now, even there are several extensions available, in this post we will implement specific scheme of using native C++ program in C# via C++/CLI wrapper program. The scheme presented here is widely known as Adapter design pattern.

Our native C++ program is using QuantLib libraries for creating custom implementations for instrument and pricing engine (zero-coupon bond). C++/CLI is then wrapping this native C++ program and exposing only selected methods for its client (C# program). For those completely new to C++/CLI language, there is a relatively comprehensive tutorial available here, written by Adam Sawicki. Needless to say, even this post is presenting relatively simple example of using QuantLib via wrapper, one is able to apply presented scheme to much more complex valuation scenarios. Sounds like a cliche, but "only the sky is the limit" applies well in here.

As far as I see, a programmer is facing complexities here on two fronts : on a pure program level (C++ language, QuantLib library) and on technical level (handling projects and their configurations specifically on Visual Studio). I assume, that the reader is already familiar with C++ language (including its modern features) and has experience on using QuantLib library. Now, in order to decrease the learning curve on that technical level, I decided to include all possible configurations-related steps in here. This means, that when all the steps described in this post have been implemented on a priestly manner, one should end up with succesfully compiled projects. In order to give some concrete back-up for this claim, I have specifically re-created this project succesfully from the scratch now two times.


C++/CLI wrapper project


Create a new C++ CRL Class library project. At this point, pre-installed and pre-built QuantLib and Boost libraries should be available. Create references to required Boost and QuantLib header files and libraries as follows.









In the case one may not have these libraries available, all required procedures for getting this part done correctly is well presented in here and here.

Next, add the following two new header files to this project.

Native C++ header file.

// MJZeroCouponBond.h
#pragma once
#include <ql/quantlib.hpp>
using namespace QuantLib;

namespace QuantLibCppNative {

 // implementation for zero-coupon bond instrument
 class MJZeroCouponBond : public Instrument {
 public:
  // forward class declarations
  class arguments;
  class engine;
  //
  // ctor and implementations for required base class methods
  MJZeroCouponBond(Real faceAmount, Date maturityDate);
  bool isExpired() const;
 private:
  void setupArguments(PricingEngine::arguments* args) const;
  // term sheet related information
  Real faceAmount_;
  Date maturityDate_;
 };

 // inner arguments class
 class MJZeroCouponBond::arguments : public PricingEngine::arguments{
 public:
  void validate() const;
  Real faceAmount;
  Date maturityDate;
 };

 // inner engine class
 class MJZeroCouponBond::engine
  : public GenericEngine < MJZeroCouponBond::arguments, MJZeroCouponBond::results > {
  // base class for all further engine implementations
 };

 // implementation for base class engine
 class MJZeroCouponBondEngine : public MJZeroCouponBond::engine {
 public:
  MJZeroCouponBondEngine(const Handle<YieldTermStructure>& curve);
  void calculate() const;
 private:
  Handle<YieldTermStructure> curve_;
 };

}

C++/CLI header file.

// Wrapper.h
#pragma once
#include "MJZeroCouponBond.h"
using namespace System;

namespace QuantLibCppWrapper {

 // C++/CLI wrapper class for native C++ QuantLib
 public ref class MJZeroCouponBondWrapper {
 public:
  MJZeroCouponBondWrapper(double riskFreeRate_, double faceAmount_, DateTime transactionDate_,
   DateTime maturityDate_, String^ calendar_, String^ daycounter_, int settlementDays_);
  !MJZeroCouponBondWrapper();
  ~MJZeroCouponBondWrapper();
  double PV();
 private:
  // note : class members as native pointers
  QuantLibCppNative::MJZeroCouponBond* bond;
  QuantLibCppNative::MJZeroCouponBondEngine* pricer;
 };
}

Next, add the following two implementation files to this project.

Native C++ implementation file.

// MJZeroCouponBond.cpp
#include "MJZeroCouponBond.h"
using namespace QuantLib;

namespace QuantLibCppNative {

 // implementations for MJZeroCouponBond instrument
 MJZeroCouponBond::MJZeroCouponBond(Real faceAmount, Date maturityDate)
  : faceAmount_(faceAmount), maturityDate_(maturityDate) { // ctor
 }

 bool MJZeroCouponBond::isExpired() const {
  return Settings::instance().evaluationDate() >= maturityDate_;
 }

 void MJZeroCouponBond::setupArguments(PricingEngine::arguments* args) const {
  MJZeroCouponBond::arguments* args_ = dynamic_cast<MJZeroCouponBond::arguments*>(args);
  QL_REQUIRE(args_ != nullptr, "arguments casting error");
  args_->faceAmount = faceAmount_;
  args_->maturityDate = maturityDate_;
 }

 void MJZeroCouponBond::arguments::validate() const {
  QL_REQUIRE(faceAmount > 0.0, "transaction face amount cannot be zero");
 }

 // implementations for MJZeroCouponBondEngine
 MJZeroCouponBondEngine::MJZeroCouponBondEngine(const Handle<YieldTermStructure>& curve)
  : curve_(curve) {
  // register observer (MJZeroCouponBondEngine) with observable (curve)
  registerWith(curve_);
 }

 void MJZeroCouponBondEngine::calculate() const {
  // extract required parameters from arguments or curve object
  // implement the actual pricing algorithm and store result
  Date maturityDate = arguments_.maturityDate;
  Real P = arguments_.faceAmount;
  DiscountFactor df = curve_->discount(maturityDate);
  Real PV = P * df;
  results_.value = PV;
 }

}

C++/CLI implementation file.

// Wrapper.cpp
#pragma once
#include "Wrapper.h"

namespace QuantLibCppWrapper {

 // ctor
 MJZeroCouponBondWrapper::MJZeroCouponBondWrapper(double riskFreeRate_, double faceAmount_, DateTime transactionDate_,
  DateTime maturityDate_, String^ calendar_, String^ daycounter_, int settlementDays_) {

  // conversion for calendar (string to QuantLib::Calendar)
  Calendar calendar;
  if (calendar_->ToUpper() == "TARGET") calendar = TARGET();
   else throw gcnew System::Exception("undefined calendar");
  
  // conversion for daycounter (string to QuantLib::DayCounter)
  DayCounter dayCounter;
  if (daycounter_->ToUpper() == "ACT360") dayCounter = Actual360();
   else throw gcnew System::Exception("undefined daycounter");

  // conversion : transaction and maturity dates (Datetime to QuantLib::Date)
  Date transactionDate(transactionDate_.ToOADate());
  Date maturityDate(maturityDate_.ToOADate());
  Date settlementDate = calendar.advance(transactionDate, Period(settlementDays_, Days));
  Settings::instance().evaluationDate() = settlementDate;

  // create flat discount curve
  auto riskFreeRate = boost::make_shared<SimpleQuote>(riskFreeRate_);
  Handle<Quote> riskFreeRateHandle(riskFreeRate);
  auto riskFreeRateTermStructure = boost::make_shared<FlatForward>(settlementDate, riskFreeRateHandle, dayCounter);
  Handle<YieldTermStructure> riskFreeRateTermStructureHandle(riskFreeRateTermStructure);

  // create instrument and pricer as native pointers
  bond = new QuantLibCppNative::MJZeroCouponBond(faceAmount_, maturityDate);
  pricer = new QuantLibCppNative::MJZeroCouponBondEngine(riskFreeRateTermStructureHandle);

  // pair instrument and pricer
  // note : cast pricer from native pointer to boost shared pointer
  bond->setPricingEngine(static_cast<boost::shared_ptr<QuantLibCppNative::MJZeroCouponBondEngine>>(pricer));
 }

 // finalizer
 MJZeroCouponBondWrapper::!MJZeroCouponBondWrapper() {
  delete bond;
  delete pricer;
 }

 // destructor
 MJZeroCouponBondWrapper::~MJZeroCouponBondWrapper() {
  this->!MJZeroCouponBondWrapper();
 }

 // return pv
 double MJZeroCouponBondWrapper::PV() {
  return bond->NPV();
 }

}

Next, some C++/CLI project settings needs to be modified.

Disable the use of pre-compiled headers.


Update properties to suppress some specific warnings.








Optionally, update properties to suppress (almost) all the other warnings.









After these steps, I have completed a succesfull built for my C++/CLI project.








C# client project


First, create a new C# console project into the existing solution.

Then, add reference to previously created C++/CLI project.








Next, implement the following program to C# project.

using System;

namespace Client {

    static class ZeroCouponBond {
        static void Main() {
            try {

                // create parameters
                double riskFreeRate = 0.01;
                double faceAmount = 1000000.0;
                DateTime transactionDate = new DateTime(2018, 4, 16);
                DateTime maturityDate = new DateTime(2020, 4, 16);
                string calendar = "TARGET";
                string daycounter = "ACT360";
                int settlementDays = 2;
                
                // use C++/CLI wrapper : create bond, request pv
                QuantLibCppWrapper.MJZeroCouponBondWrapper zero = new QuantLibCppWrapper.MJZeroCouponBondWrapper(
                    riskFreeRate, faceAmount, transactionDate, maturityDate, calendar, daycounter, settlementDays);
                double PV = zero.PV();
                Console.WriteLine(PV.ToString());
                GC.SuppressFinalize(zero);
            }
            catch (Exception e) {
                Console.WriteLine(e.Message);
            }
        }
    }
}

Finally, set C# project as start-up project (on Visual Studio, right-click selected C# project and select "Set as StartUp Project").

After completing all previous steps, I have completed a succesfull built for the both projects and got the result of 979953.65 as present value for this 2-year zero-coupon bond, evaluated by using our native C++ QuantLib program via C++/CLI wrapper class.








At this point, we are done.

Postlude : no sweeping under the carpet allowed


Clever eyes might be catching, that there is a SuppressFinalize method call made for Garbage Collector at the end of our C# program. If I will remove that call, I will get the following exception.






Now, the both destructor and finalizer in our C++/CLI project have been correctly implemented, as C++/CLI standard is suggesting. As I understand, the issue is coming from C# program side, as Garbage Collector will do its final memory release sweep before exit. The issue has been pretty completely chewed in here, here and here.

Interesting point is, that in the last given reference above, the author is actually suggesting to delete dynamically allocated member variables (bond and pricer) in destructor. Now, if I will modify that C++/CLI program accordingly as suggested (delete member variables in destructor, not in finalizer, remove trigger call from destructor to finalizer, remove method call for Garbage Collector in C# program), this program will work again as expected. By taking a look at author's past experience, at least I would come to the conclusion, that this suggested approach is also a way to go.

Finally, thanks a lot again for using your precious time for reading this blog. I hope you got what you were looking for.
-Mike

Wednesday, September 6, 2017

QuantLib : Hull-White one-factor model calibration

Sometimes during the last year I published one post on simulating Hull-White interest rate paths using Quantlib. My conclusion was, that with all the tools provided by this wonderful library, this task should be (relatively) easy thing to do. However, as we know the devil is always in the details - in this case, in the actual process parameters (mean reversion, sigma). Before processing any simulation, we have to get those parameters from somewhere. In this post, we use Quantlib tools for calibrating our model parameters to swaption prices existing in the market.

There are so many variations of this model out there (depending on time-dependencies of different parameters), but the following is the model we are going to use in this example.




Parameters alpha and sigma are constants and theta is time-dependent variable (usually calibrated to current yield curve). Swaption surface is presented in the picture below. Within the table below, time to swaption maturity has been defined in vertical axis, while tenor for underlying swap contract has been defined in horizontal axis. Co-terminal swaptions (used in calibration process) have been specifically marked with yellow colour.























Calibration results


Test cases for three different calibration schemes are included in the example program. More specifically, we :
  1. calibrate the both parameters
  2. calibrate sigma parameter and freeze reversion parameter to 0.05
  3. calibrate reversion parameter and freeze sigma parameter to 0.01
With the given data, we get the following results for these calibration schemes.









The program


As preparatory task, builder class for constructing yield curve should be implemented into a new project from here. After this task, the following header file (ModelCalibrator.h) and tester file (Tester.cpp) should be added to this new project.

First, piecewise yield curve (USD swap curve) and swaption volatilities (co-terminal swaptions) are created by using two free functions in tester. All the data has been hard-coded inside these functions. Needless to say, in any production-level program, this data feed should come from somewhere else. After required market data and ModelCalibrator object has been created, calibration helpers for diagonal swaptions are going to be created and added to ModelCalibrator. Finally, test cases for different calibration schemes are processed.

// ModelCalibrator.h
#pragma once
#include <ql\quantlib.hpp>
#include <algorithm>
//
namespace MJModelCalibratorNamespace
{
 using namespace QuantLib;
 //
 template <typename MODEL, typename OPTIMIZER = LevenbergMarquardt>
 class ModelCalibrator
 {
 public:
  // default implementation (given values) for end criteria class
  ModelCalibrator(const EndCriteria& endCriteria = EndCriteria(1000, 500, 0.0000001, 0.0000001, 0.0000001))
   : endCriteria(endCriteria)
  { }
  void AddCalibrationHelper(boost::shared_ptr<CalibrationHelper>& helper)
  {
   // add any type of calibration helper
   helpers.push_back(helper);
  }
  void Calibrate(boost::shared_ptr<MODEL>& model,
   const boost::shared_ptr<PricingEngine>& pricingEngine,
   const Handle<YieldTermStructure>& curve,
   const std::vector<bool> fixedParameters = std::vector<bool>())
  {
   // assign pricing engine to all calibration helpers
   std::for_each(helpers.begin(), helpers.end(),
    [&pricingEngine](boost::shared_ptr<CalibrationHelper>& helper)
   { helper->setPricingEngine(pricingEngine); });
   //
   // create optimization model for calibrating requested model
   OPTIMIZER solver;
   //
   if (fixedParameters.empty())
   {
    // calibrate all involved parameters
    model->calibrate(helpers, solver, this->endCriteria);
   }
   else
   {
    // calibrate all involved non-fixed parameters
    // hard-coded : vector for weights and constraint type
    model->calibrate(helpers, solver, this->endCriteria, NoConstraint(), std::vector<Real>(), fixedParameters);
   }
  }
 private:
  EndCriteria endCriteria;
  std::vector<boost::shared_ptr<CalibrationHelper>> helpers;
 };
}
//
//
//
//
// Tester.cpp
#include "PiecewiseCurveBuilder.cpp"
#include "ModelCalibrator.h"
#include <iostream>
//
using namespace MJPiecewiseCurveBuilderNamespace;
namespace MJ_Calibrator = MJModelCalibratorNamespace;
//
// declarations for two free functions used to construct required market data
void CreateCoTerminalSwaptions(std::vector<Volatility>& diagonal);
void CreateYieldCurve(RelinkableHandle<YieldTermStructure>& curve, 
 const Date& settlementDate, const Calendar& calendar);
//
//
//
int main()
{
 try
 {
  // dates
  Date tradeDate(4, September, 2017);
  Settings::instance().evaluationDate() = tradeDate;
  Calendar calendar = TARGET();
  Date settlementDate = calendar.advance(tradeDate, Period(2, Days), ModifiedFollowing);
  //
  // market data : create piecewise yield curve
  RelinkableHandle<YieldTermStructure> curve;
  CreateYieldCurve(curve, settlementDate, calendar);
  //
  // market data : create co-terminal swaption volatilities
  std::vector<Volatility> diagonal;
  CreateCoTerminalSwaptions(diagonal);
  //
  // create model calibrator
  MJ_Calibrator::ModelCalibrator<HullWhite> modelCalibrator;
  //
  // create and add calibration helpers to model calibrator
  boost::shared_ptr<IborIndex> floatingIndex(new USDLibor(Period(3, Months), curve));
  for (unsigned int i = 0; i != diagonal.size(); ++i)
  {
   int timeToMaturity = i + 1;
   int underlyingTenor = diagonal.size() - i;
   //
   // using 1st constructor for swaption helper class
   modelCalibrator.AddCalibrationHelper(boost::shared_ptr<CalibrationHelper>(new SwaptionHelper(
     Period(timeToMaturity, Years), // time to swaption maturity
     Period(underlyingTenor, Years), // tenor of the underlying swap
     Handle<Quote>(boost::shared_ptr<Quote>(new SimpleQuote(diagonal[i]))), // swaption volatility
     floatingIndex, // underlying floating index
     Period(1, Years), // tenor for underlying fixed leg
     Actual360(), // day counter for underlying fixed leg
     floatingIndex->dayCounter(), // day counter for underlying floating leg
     curve))); // term structure
  }
  //
  // create model and pricing engine, calibrate model and print calibrated parameters
  // case 1 : calibrate all involved parameters (HW1F : reversion, sigma)
  boost::shared_ptr<HullWhite> model(new HullWhite(curve));
  boost::shared_ptr<PricingEngine> jamshidian(new JamshidianSwaptionEngine(model));
  modelCalibrator.Calibrate(model, jamshidian, curve);
  std::cout << "calibrated reversion = " << model->params()[0] << std::endl;
  std::cout << "calibrated sigma = " << model->params()[1] << std::endl;
  std::cout << std::endl;
  //
  // case 2 : calibrate sigma and fix reversion to famous 0.05
  model = boost::shared_ptr<HullWhite>(new HullWhite(curve, 0.05, 0.0001));
  jamshidian = boost::shared_ptr<PricingEngine>(new JamshidianSwaptionEngine(model));
  std::vector<bool> fixedReversion = { true, false };
  modelCalibrator.Calibrate(model, jamshidian, curve, fixedReversion);
  std::cout << "fixed reversion = " << model->params()[0] << std::endl;
  std::cout << "calibrated sigma = " << model->params()[1] << std::endl;
  std::cout << std::endl;
  //
  // case 3 : calibrate reversion and fix sigma to 0.01
  model = boost::shared_ptr<HullWhite>(new HullWhite(curve, 0.05, 0.01));
  jamshidian = boost::shared_ptr<PricingEngine>(new JamshidianSwaptionEngine(model));
  std::vector<bool> fixedSigma = { false, true };
  modelCalibrator.Calibrate(model, jamshidian, curve, fixedSigma);
  std::cout << "calibrated reversion = " << model->params()[0] << std::endl;
  std::cout << "fixed sigma = " << model->params()[1] << std::endl;
 }
 catch (std::exception& e)
 {
  std::cout << e.what() << std::endl;
 }
 return 0;
}
//
void CreateCoTerminalSwaptions(std::vector<Volatility>& diagonal)
{
 // hard-coded data
 // create co-terminal swaptions 
 diagonal.push_back(0.3133); // 1x10
 diagonal.push_back(0.3209); // 2x9
 diagonal.push_back(0.3326); // 3x8
 diagonal.push_back(0.331); // 4x7
 diagonal.push_back(0.3281); // 5x6
 diagonal.push_back(0.318); // 6x5
 diagonal.push_back(0.3168); // 7x4
 diagonal.push_back(0.3053); // 8x3
 diagonal.push_back(0.2992); // 9x2
 diagonal.push_back(0.3073); // 10x1
}
//
void CreateYieldCurve(RelinkableHandle<YieldTermStructure>& curve,
 const Date& settlementDate, const Calendar& calendar) 
{
 // hard-coded data
 // create piecewise yield curve by using builder class
 DayCounter curveDaycounter = Actual360();
 PiecewiseCurveBuilder<ZeroYield, Linear> builder;
 //
 // cash rates
 pQuote q_1W(new SimpleQuote(0.012832));
 pIndex i_1W(new USDLibor(Period(1, Weeks)));
 builder.AddDeposit(q_1W, i_1W);
 //
 pQuote q_1M(new SimpleQuote(0.012907));
 pIndex i_1M(new USDLibor(Period(1, Months)));
 builder.AddDeposit(q_1M, i_1M);
 //
 pQuote q_3M(new SimpleQuote(0.0131611));
 pIndex i_3M(new USDLibor(Period(3, Months))); 
 builder.AddDeposit(q_3M, i_3M);
 //
 // futures
 Date IMMDate;
 pQuote q_DEC17(new SimpleQuote(98.5825));
 IMMDate = IMM::nextDate(settlementDate + Period(3, Months));
 builder.AddFuture(q_DEC17, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 pQuote q_MAR18(new SimpleQuote(98.5425));
 IMMDate = IMM::nextDate(settlementDate + Period(6, Months));
 builder.AddFuture(q_MAR18, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 pQuote q_JUN18(new SimpleQuote(98.4975));
 IMMDate = IMM::nextDate(settlementDate + Period(9, Months));
 builder.AddFuture(q_JUN18, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 pQuote q_SEP18(new SimpleQuote(98.4475));
 IMMDate = IMM::nextDate(settlementDate + Period(12, Months));
 builder.AddFuture(q_SEP18, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 pQuote q_DEC18(new SimpleQuote(98.375));
 IMMDate = IMM::nextDate(settlementDate + Period(15, Months));
 builder.AddFuture(q_DEC18, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 pQuote q_MAR19(new SimpleQuote(98.3425));
 IMMDate = IMM::nextDate(settlementDate + Period(18, Months));
 builder.AddFuture(q_MAR19, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 pQuote q_JUN19(new SimpleQuote(98.3025));
 IMMDate = IMM::nextDate(settlementDate + Period(21, Months));
 builder.AddFuture(q_JUN19, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 pQuote q_SEP19(new SimpleQuote(98.2675));
 IMMDate = IMM::nextDate(settlementDate + Period(24, Months));
 builder.AddFuture(q_SEP19, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 pQuote q_DEC19(new SimpleQuote(98.2125));
 IMMDate = IMM::nextDate(settlementDate + Period(27, Months));
 builder.AddFuture(q_DEC19, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 pQuote q_MAR20(new SimpleQuote(98.1775));
 IMMDate = IMM::nextDate(settlementDate + Period(30, Months));
 builder.AddFuture(q_MAR20, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 pQuote q_JUN20(new SimpleQuote(98.1425));
 IMMDate = IMM::nextDate(settlementDate + Period(33, Months));
 builder.AddFuture(q_JUN20, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
 //
 // swaps
 pIndex swapFloatIndex(new USDLibor(Period(3, Months))); 
 pQuote q_4Y(new SimpleQuote(0.01706));
 builder.AddSwap(q_4Y, Period(4, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_5Y(new SimpleQuote(0.0176325));
 builder.AddSwap(q_5Y, Period(5, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_6Y(new SimpleQuote(0.01874));
 builder.AddSwap(q_6Y, Period(6, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_7Y(new SimpleQuote(0.0190935));
 builder.AddSwap(q_7Y, Period(7, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_8Y(new SimpleQuote(0.02011));
 builder.AddSwap(q_8Y, Period(8, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_9Y(new SimpleQuote(0.02066));
 builder.AddSwap(q_9Y, Period(9, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_10Y(new SimpleQuote(0.020831));
 builder.AddSwap(q_10Y, Period(10, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_11Y(new SimpleQuote(0.02162));
 builder.AddSwap(q_11Y, Period(11, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_12Y(new SimpleQuote(0.0217435));
 builder.AddSwap(q_12Y, Period(12, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_15Y(new SimpleQuote(0.022659));
 builder.AddSwap(q_15Y, Period(15, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_20Y(new SimpleQuote(0.0238125));
 builder.AddSwap(q_20Y, Period(20, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_25Y(new SimpleQuote(0.0239385));
 builder.AddSwap(q_25Y, Period(25, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 pQuote q_30Y(new SimpleQuote(0.02435));
 builder.AddSwap(q_30Y, Period(30, Years), calendar, Annual, ModifiedFollowing, Actual360(), swapFloatIndex);
 //
 curve = builder.GetCurveHandle(settlementDate, curveDaycounter);
}


I have noticed, that there are some rule-of-thumbs, but the actual calibration for any of these models is not so straightforward as one may think. There are some subtle issues on market data quality and products under pricing, which are leaving us with relatively high degree of freedom. Further step towards the calibration red pill can be taken by checking out this excellent research paper, written by the guys from Mizuho Securities.

Finally, as always, thanks a lot again for reading this blog.
-Mike

Sunday, September 3, 2017

QuantLib : another implementation for piecewise yield curve builder class


Last year I published one possible implementation using Quantlib library for constructing piecewise yield curves. Within this second implementation, I have done a couple of changes in order to increase configurability. For allowing more flexible curve construction algorithm, Traits and Interpolations are now defined as template parameters. Previously, these were hard-coded inside the class. Secondly, all types of quotes for class methods are now wrapped inside shared pointers. Previously, there was an option to give quotes as rates. After some reconsideration, I have decided to give up this option completely. Finally, I have added a class method for allowing futures prices to be used in curve construction process.


Source data and results


The following source data from the book by Richard Flavell has been used for validation.



















Resulting zero-yield curve and discount factors are presented in the graph below. For this data used, the absolute maximum difference between Flavell-constructed and Quantlib-constructed zero yield curves is around one basis point.

















The program


// PiecewiseCurveBuilder.h
#pragma once
#include <ql/quantlib.hpp>
//
namespace MJPiecewiseCurveBuilderNamespace
{
 using namespace QuantLib;
 //
 // type alias definitions
 using pQuote = boost::shared_ptr<Quote>;
 using pIndex = boost::shared_ptr<IborIndex>;
 using pHelper = boost::shared_ptr<RateHelper>;
 //
 // T = traits, I = interpolation
 template<typename T, typename I>
 class PiecewiseCurveBuilder
 {
 public:
  PiecewiseCurveBuilder();
  void AddDeposit(const pQuote& quote, const pIndex& index);
  void AddFRA(const pQuote& quote, const Period& periodLengthToStart, const pIndex& index);
  void AddSwap(const pQuote& quote, const Period& periodLength, const Calendar& fixedCalendar, Frequency fixedFrequency,
   BusinessDayConvention fixedConvention, const DayCounter& fixedDayCount, const pIndex& floatIndex);
  //
  void AddFuture(const pQuote& quote, const Date& IMMDate, int lengthInMonths, const Calendar& calendar,
   BusinessDayConvention convention, const DayCounter& dayCounter, bool endOfMonth = true);
  //
  RelinkableHandle<YieldTermStructure> GetCurveHandle(const Date& settlementDate, const DayCounter& dayCounter);
 private:
  std::vector<pHelper> rateHelpers;

 };
}
//
//
//
//
// PiecewiseCurveBuilder.cpp
#pragma once
#include "PiecewiseCurveBuilder.h"
//
namespace MJPiecewiseCurveBuilderNamespace
{
 template<typename T, typename I>
 PiecewiseCurveBuilder<T, I>::PiecewiseCurveBuilder() { }
 //
 template<typename T, typename I>
 void PiecewiseCurveBuilder<T, I>::AddDeposit(const pQuote& quote, const pIndex& index)
 {
  pHelper rateHelper(new DepositRateHelper(Handle<Quote>(quote), index));
  rateHelpers.push_back(rateHelper);
 }
 //
 template<typename T, typename I>
 void PiecewiseCurveBuilder<T, I>::AddFRA(const pQuote& quote, const Period& periodLengthToStart, const pIndex& index)
 {
  pHelper rateHelper(new FraRateHelper(Handle<Quote>(quote),
   periodLengthToStart, pIndex(index)));
  rateHelpers.push_back(rateHelper);
 }
 //
 template<typename T, typename I>
 void PiecewiseCurveBuilder<T, I>::AddSwap(const pQuote& quote, const Period& periodLength, const Calendar& fixedCalendar,
  Frequency fixedFrequency, BusinessDayConvention fixedConvention, const DayCounter& fixedDayCount,
  const pIndex& floatIndex)
 {
  pHelper rateHelper(new SwapRateHelper(Handle<Quote>(quote), periodLength, fixedCalendar, fixedFrequency,
   fixedConvention, fixedDayCount, floatIndex));
  rateHelpers.push_back(rateHelper);
 }
 //
 template<typename T, typename I>
 void PiecewiseCurveBuilder<T, I>::AddFuture(const pQuote& quote, const Date& IMMDate, int lengthInMonths, const Calendar& calendar,
  BusinessDayConvention convention, const DayCounter& dayCounter, bool endOfMonth)
 {
  pHelper rateHelper(new FuturesRateHelper(Handle<Quote>(quote), IMMDate, lengthInMonths, calendar, convention, endOfMonth, dayCounter));
  rateHelpers.push_back(rateHelper);
 }
 //
 template<typename T, typename I>
 RelinkableHandle<YieldTermStructure> PiecewiseCurveBuilder<T, I>::GetCurveHandle(const Date& settlementDate, const DayCounter& dayCounter)
 {
  // T = traits, I = interpolation
  boost::shared_ptr<YieldTermStructure> yieldTermStructure(new PiecewiseYieldCurve<T, I>(settlementDate, rateHelpers, dayCounter));
  return RelinkableHandle<YieldTermStructure>(yieldTermStructure);
 }
}
//
//
//
//
// Tester.cpp
#include "PiecewiseCurveBuilder.cpp"
#include <iostream>
using namespace MJPiecewiseCurveBuilderNamespace;
//
int main()
{
 try
 {
  Date tradeDate(4, February, 2008);
  Settings::instance().evaluationDate() = tradeDate;
  Calendar calendar = TARGET();
  Date settlementDate = calendar.advance(tradeDate, Period(2, Days), ModifiedFollowing);
  DayCounter curveDaycounter = Actual360();
  PiecewiseCurveBuilder<ZeroYield, Linear> builder;
  //
  //
  // cash part of the curve
  pQuote q_1W(new SimpleQuote(0.032175));
  pIndex i_1W(new USDLibor(Period(1, Weeks)));
  builder.AddDeposit(q_1W, i_1W);
  //
  pQuote q_1M(new SimpleQuote(0.0318125));
  pIndex i_1M(new USDLibor(Period(1, Months)));
  builder.AddDeposit(q_1M, i_1M);
  //
  pQuote q_3M(new SimpleQuote(0.03145));
  pIndex i_3M(new USDLibor(Period(3, Months)));
  builder.AddDeposit(q_3M, i_3M);
  //
  //
  // futures part of the curve
  Date IMMDate;
  pQuote q_JUN08(new SimpleQuote(97.41));
  IMMDate = IMM::nextDate(settlementDate + Period(4, Months));
  builder.AddFuture(q_JUN08, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
  //
  pQuote q_SEP08(new SimpleQuote(97.52));
  IMMDate = IMM::nextDate(settlementDate + Period(7, Months));
  builder.AddFuture(q_SEP08, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
  //
  pQuote q_DEC08(new SimpleQuote(97.495));
  IMMDate = IMM::nextDate(settlementDate + Period(10, Months));
  builder.AddFuture(q_DEC08, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
  //
  pQuote q_MAR09(new SimpleQuote(97.395));
  IMMDate = IMM::nextDate(settlementDate + Period(13, Months));
  builder.AddFuture(q_MAR09, IMMDate, 3, calendar, ModifiedFollowing, Actual360());
  //
  //
  // swap part of the curve
  pIndex swapFloatIndex(new USDLibor(Period(3, Months)));
  pQuote q_2Y(new SimpleQuote(0.02795));
  builder.AddSwap(q_2Y, Period(2, Years), calendar, Annual,
   ModifiedFollowing, Actual360(), swapFloatIndex);
  //
  pQuote q_3Y(new SimpleQuote(0.03035));
  builder.AddSwap(q_3Y, Period(3, Years), calendar, Annual,
   ModifiedFollowing, Actual360(), swapFloatIndex);
  //
  pQuote q_4Y(new SimpleQuote(0.03275));
  builder.AddSwap(q_4Y, Period(4, Years), calendar, Annual,
   ModifiedFollowing, Actual360(), swapFloatIndex);
  //
  pQuote q_5Y(new SimpleQuote(0.03505));
  builder.AddSwap(q_5Y, Period(5, Years), calendar, Annual,
   ModifiedFollowing, Actual360(), swapFloatIndex);
  //
  pQuote q_6Y(new SimpleQuote(0.03715));
  builder.AddSwap(q_6Y, Period(6, Years), calendar, Annual,
   ModifiedFollowing, Actual360(), swapFloatIndex);
  //
  pQuote q_7Y(new SimpleQuote(0.03885));
  builder.AddSwap(q_7Y, Period(7, Years), calendar, Annual,
   ModifiedFollowing, Actual360(), swapFloatIndex);
  //
  pQuote q_8Y(new SimpleQuote(0.04025));
  builder.AddSwap(q_8Y, Period(8, Years), calendar, Annual,
   ModifiedFollowing, Actual360(), swapFloatIndex);
  //
  pQuote q_9Y(new SimpleQuote(0.04155));
  builder.AddSwap(q_9Y, Period(9, Years), calendar, Annual,
   ModifiedFollowing, Actual360(), swapFloatIndex);
  //
  pQuote q_10Y(new SimpleQuote(0.04265));
  builder.AddSwap(q_10Y, Period(10, Years), calendar, Annual,
   ModifiedFollowing, Actual360(), swapFloatIndex);
  //
  pQuote q_12Y(new SimpleQuote(0.04435));
  builder.AddSwap(q_12Y, Period(12, Years), calendar, Annual,
   ModifiedFollowing, Actual360(), swapFloatIndex);
  //
  //
  // get curve handle and print out discount factors
  RelinkableHandle<YieldTermStructure> curve = builder.GetCurveHandle(settlementDate, curveDaycounter);
  std::cout << curve->discount(Date(11, February, 2008)) << std::endl;
  std::cout << curve->discount(Date(4, March, 2008)) << std::endl;
  std::cout << curve->discount(Date(4, May, 2008)) << std::endl;
  std::cout << curve->discount(Date(6, August, 2008)) << std::endl;
  std::cout << curve->discount(Date(6, November, 2008)) << std::endl;
  std::cout << curve->discount(Date(6, February, 2009)) << std::endl;
  std::cout << curve->discount(Date(8, February, 2010)) << std::endl;
  std::cout << curve->discount(Date(7, February, 2011)) << std::endl;
  std::cout << curve->discount(Date(6, February, 2012)) << std::endl;
  std::cout << curve->discount(Date(6, February, 2013)) << std::endl;
  std::cout << curve->discount(Date(6, February, 2014)) << std::endl;
  std::cout << curve->discount(Date(6, February, 2015)) << std::endl;
  std::cout << curve->discount(Date(8, February, 2016)) << std::endl;
  std::cout << curve->discount(Date(6, February, 2017)) << std::endl;
  std::cout << curve->discount(Date(6, February, 2018)) << std::endl;
 }
 catch (std::exception& e)
 {
  std::cout << e.what() << std::endl;
 }
 return 0;
}


Data updating


Since all rates have been wrapped inside quotes (which have been wrapped inside handles), those can be accessed only by using dynamic pointer casting. Below is an example for updating 3M cash quote.

boost::dynamic_pointer_cast<SimpleQuote>(q_3M)->setValue(0.03395);

After this quote updating shown above, a new requested value (zero rate, forward rate or discount factor) from constructed curve object will be re-calculated by using updated quote.

For those readers who are completely unfamiliar with this stuff presented, there are three extremely well-written slides available in Quantlib documentation page, written by Dimitri Reiswich. These slides are offering excellent way to get very practical hands-on overview on QuantLib library. Also, Luigi Ballabio has finally been finishing his book on QuantLib implementation, which can be purchased from Leanpub. This book is offering deep diving experience into the abyss of Quantlib architechture.

Finally, as usual, thanks for spending your precious time here and reading this blog.
-Mike

Wednesday, April 27, 2016

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