QuantLib : implementing Equity-linked note using Monte Carlo framework

Within the last post, an implementation for a simple custom instrument and analytical pricing engine was presented. This post is presenting a bit more complex implementation for an equity-linked note, using custom pricing engine implementation built on the top of QuantLib Monte Carlo framework.

Term sheet

The most essential transaction details have been presented within the following screenshot below.

In a nutshell, 3-year holding period has been divided into three annual periods. For each period, running cumulative (but capped) coupon will be calculated based on simulated index fixing values (Geometric Brownian Motion). At the end of each period, period payoff (floored) will be calculated. Pricing-wise, this period payoff will then be discounted to present valuation date. Structure total PV is the sum of all discounted period payoffs.

In order to get more familiar with QuantLib Monte Carlo framework, one may take a look at Implementing QuantLib blog or get the actual book on QuantLib implementation from Leanpub. Also, MoneyScience is organizing Introduction to QuantLib Development course regularly. The most essential parts of the library are thoroughly covered during this 3-day training course, hosted by QuantLib lead developer Luigi Ballabio. Thanks for reading this blog. Merry Christmas and Happy New Year for everybody.

-Mike

// EquityLinkedNote.h
#pragma once
#include <ql/quantlib.hpp>
using namespace QuantLib;
//
// instrument implementation for equity-linked note
class EquityLinkedNote : public Instrument {
public:
 // forward class declarations
 class arguments;
 class engine;
 //
 // ctor and implementations for required base class methods
 EquityLinkedNote(Real notional, Real initialFixing, const std::vector<Date>& fixingDates, 
  const std::vector<Date>& paymentDates, Real cap, Real floor);
 bool isExpired() const;
private:
 void setupArguments(PricingEngine::arguments* args) const;
 // term sheet information
 Real notional_;
 Real initialFixing_;
 std::vector<Date> fixingDates_;
 std::vector<Date> paymentDates_;
 Real cap_;
 Real floor_;
};
// inner arguments class
class EquityLinkedNote::arguments : public PricingEngine::arguments{
public:
 void validate() const;
 Real notional;
 Real initialFixing;
 std::vector<Date> fixingDates;
 std::vector<Date> paymentDates;
 Real cap;
 Real floor;
};
// inner engine class
class EquityLinkedNote::engine
 : public GenericEngine<EquityLinkedNote::arguments, EquityLinkedNote::results> {
 // base class for all further engine implementations
};
//
//
// path pricer implementation for equity-linked note
class EquityLinkedNotePathPricer : public PathPricer<Path> {
public:
 EquityLinkedNotePathPricer(Real notional, Real initialFixing, const std::vector<Date>& fixingDates,
  const std::vector<Date>& paymentDates, Real cap, Real floor, const Handle<YieldTermStructure>& curve);
 Real operator()(const Path& path) const;
private:
 Real notional_;
 Real initialFixing_;
 std::vector<Date> fixingDates_;
 std::vector<Date> paymentDates_;
 Real cap_;
 Real floor_;
 Handle<YieldTermStructure> curve_;
};
//
//
// monte carlo framework engine implementation for base engine class
template <typename RNG = PseudoRandom, typename S = Statistics>
class EquityLinkedNoteMonteCarloPricer : public EquityLinkedNote::engine, private McSimulation<SingleVariate, RNG, S> {
public:
 // ctor
 EquityLinkedNoteMonteCarloPricer(const boost::shared_ptr<StochasticProcess>& process,
  const Handle<YieldTermStructure>& curve, bool antitheticVariate, Size requiredSamples,
  Real requiredTolerance, Size maxSamples, BigNatural seed)
  : McSimulation<SingleVariate, RNG, S>(antitheticVariate, false), process_(process), curve_(curve),
  requiredSamples_(requiredSamples), requiredTolerance_(requiredTolerance), maxSamples_(maxSamples), seed_(seed) {
   // register observer (pricer) with observables (stochastic process, curve)
   registerWith(process_);
   registerWith(curve_);
  }
 // implementation for required base engine class method
 void calculate() const {
  // the actual simulation process will be performed within the following method
  McSimulation<SingleVariate, RNG, S>::calculate(requiredTolerance_, requiredSamples_, maxSamples_);
  this->results_.value = this->mcModel_->sampleAccumulator().mean();
  //
  if (RNG::allowsErrorEstimate) {
   this->results_.errorEstimate = this->mcModel_->sampleAccumulator().errorEstimate();
  }
  else {
   this->results_.errorEstimate = Null<Real>();
  }
 }
private:
 // type definitions
 typedef McSimulation<SingleVariate, RNG, S> simulation;
 typedef typename simulation::path_pricer_type path_pricer_type;
 typedef typename simulation::path_generator_type path_generator_type;
 //
 // implementation for McSimulation class virtual method
 TimeGrid timeGrid() const {
  // create time grid based on a set of given index fixing dates
  Date referenceDate = curve_->referenceDate();
  DayCounter dayCounter = curve_->dayCounter();
  std::vector<Time> fixingTimes(arguments_.fixingDates.size());
  for (Size i = 0; i != fixingTimes.size(); ++i) {
   fixingTimes[i] = dayCounter.yearFraction(referenceDate, arguments_.fixingDates[i]);
  }
  return TimeGrid(fixingTimes.begin(), fixingTimes.end());
 }
 // implementation for McSimulation class virtual method
 boost::shared_ptr<path_generator_type> pathGenerator() const {
  // create time grid and get information concerning number of simulation steps for a path
  TimeGrid grid = timeGrid();
  Size steps = (grid.size() - 1);
  // create random sequence generator and return path generator
  typename RNG::rsg_type generator = RNG::make_sequence_generator(steps, seed_);
  return boost::make_shared<path_generator_type>(process_, grid, generator, false);
 }
 // implementation for McSimulation class virtual method
 boost::shared_ptr<path_pricer_type> pathPricer() const {
  // create path pricer implementation for equity-linked note
  return boost::make_shared<EquityLinkedNotePathPricer>(arguments_.notional, arguments_.initialFixing, 
   arguments_.fixingDates, arguments_.paymentDates, arguments_.cap, arguments_.floor, this->curve_);
 }
 // private data members
 boost::shared_ptr<StochasticProcess> process_;
 Handle<YieldTermStructure> curve_;
 Size requiredSamples_;
 Real requiredTolerance_;
 Size maxSamples_;
 BigNatural seed_;
};
//
//
//
//
//
// EquityLinkedNote.cpp
#include "EquityLinkedNote.h"
#include <algorithm>
//
// implementations for equity-linked note methods
EquityLinkedNote::EquityLinkedNote(Real notional, Real initialFixing, const std::vector<Date>& fixingDates,
 const std::vector<Date>& paymentDates, Real cap, Real floor)
 : notional_(notional), initialFixing_(initialFixing), fixingDates_(fixingDates), 
 paymentDates_(paymentDates), cap_(cap), floor_(floor) {
 // ctor
}
bool EquityLinkedNote::isExpired() const {
 Date valuationDate = Settings::instance().evaluationDate();
 // note is expired, if valuation date is greater than the last fixing date
 if (valuationDate > fixingDates_.back())
  return true;
 return false;
}
void EquityLinkedNote::setupArguments(PricingEngine::arguments* args) const {
 EquityLinkedNote::arguments* args_ = dynamic_cast<EquityLinkedNote::arguments*>(args);
 QL_REQUIRE(args_ != nullptr, "arguments casting error");
 args_->notional = notional_;
 args_->initialFixing = initialFixing_;
 args_->fixingDates = fixingDates_;
 args_->paymentDates = paymentDates_;
 args_->cap = cap_;
 args_->floor = floor_;
}
void EquityLinkedNote::arguments::validate() const {
 // checks for some general argument properties
 QL_REQUIRE(notional > 0.0, "notional must be greater than zero");
 QL_REQUIRE(initialFixing > 0.0, "initial fixing must be greater than zero");
 QL_REQUIRE(cap > floor, "cap must be greater than floor");
 // check for date consistency : all payment dates have to be included 
 // within a set of given fixing dates
 for (int i = 0; i != paymentDates.size(); ++i) {
  if (std::find(fixingDates.begin(), fixingDates.end(), paymentDates[i]) == fixingDates.end()) {
   QL_REQUIRE(false, "payment date has to be included within given fixing dates");
  }
 }
}
//
// implementations for equity-linked path pricer methods
EquityLinkedNotePathPricer::EquityLinkedNotePathPricer(Real notional, Real initialFixing, const std::vector<Date>& fixingDates,
 const std::vector<Date>& paymentDates, Real cap, Real floor, const Handle<YieldTermStructure>& curve)
 : notional_(notional), initialFixing_(initialFixing), fixingDates_(fixingDates),
 paymentDates_(paymentDates), cap_(cap), floor_(floor), curve_(curve) {
 // ctor
}
// the actual pricing algorithm for a simulated path is implemented in this method
Real EquityLinkedNotePathPricer::operator()(const Path& path) const {
 Real coupon = 0.0;
 Real cumulativeCoupon = 0.0;
 Real aggregatePathPayoff = 0.0;
 int paymentDateCounter = 0;
 //
 // loop through fixing dates
 for (int i = 1; i != fixingDates_.size(); ++i) {
  // calculate floating coupon, based on simulated index fixing values
  coupon = std::min(path.at(i) / path.at(i - 1) - 1, cap_);
  // add floating coupon to cumulative coupon
  cumulativeCoupon += coupon;
  //
  // calculate period payoff for each payment date
  if (fixingDates_[i] == paymentDates_[paymentDateCounter]) {
   // calculate discounted payoff for current period, add value to aggregate path payoff
   aggregatePathPayoff += std::max(cumulativeCoupon, floor_) * notional_ * curve_->discount(fixingDates_[i]);
   // re-initialize cumulative coupon to zero, look for the next payment date
   cumulativeCoupon = 0.0;
   paymentDateCounter++;
  }
 }
 return aggregatePathPayoff;
}
//
//
//
//
//
// main.cpp
#include "EquityLinkedNote.h"
//
int main() {
 try {
  // common data : calendar, daycounter, dates
  Calendar calendar = TARGET();
  DayCounter dayCounter = Actual360();
  Date transactionDate(30, October, 2017);
  Natural settlementDays = 2;
  Date settlementDate = calendar.advance(transactionDate, Period(settlementDays, Days));
  Settings::instance().evaluationDate() = settlementDate;
  //
  // term sheet parameters
  Real notional = 1000000.0;
  Real initialFixing = 3662.18;
  Real cap = 0.015;
  Real floor = 0.0;
  std::vector<Date> fixingDates {
   Date(30, November, 2017), Date(30, December, 2017), Date(30, January, 2018), 
   Date(28, February, 2018), Date(30, March, 2018), Date(30, April, 2018),
   Date(30, May, 2018), Date(30, June, 2018), Date(30, July, 2018), 
   Date(30, August, 2018), Date(30, September, 2018), Date(30, October, 2018),
   Date(30, November, 2018), Date(30, December, 2018), Date(30, January, 2019), 
   Date(28, February, 2019), Date(30, March, 2019), Date(30, April, 2019),
   Date(30, May, 2019), Date(30, June, 2019), Date(30, July, 2019), 
   Date(30, August, 2019), Date(30, September, 2019), Date(30, October, 2019),
   Date(30, November, 2019), Date(30, December, 2019), Date(30, January, 2020), 
   Date(29, February, 2020), Date(30, March, 2020), Date(30, April, 2020),
   Date(30, May, 2020), Date(30, June, 2020), Date(30, July, 2020), 
   Date(30, August, 2020), Date(30, September, 2020), Date(30, October, 2020)
  };
  std::vector<Date> paymentDates {
   Date(30, October, 2018), Date(30, October, 2019), Date(30, October, 2020)
  };
  //
  // create structured equity-linked note
  auto note = boost::make_shared<EquityLinkedNote>(notional, initialFixing, fixingDates, paymentDates, cap, floor);
  //
  // market data
  // create discount curve
  Real riskFreeRate = 0.01;
  auto riskFreeRateQuote = boost::make_shared<SimpleQuote>(riskFreeRate);
  Handle<Quote> riskFreeRateHandle(riskFreeRateQuote);
  auto riskFreeRateTermStructure = boost::make_shared<FlatForward>(settlementDate, riskFreeRateHandle, dayCounter);
  Handle<YieldTermStructure> riskFreeRateTermStructureHandle(riskFreeRateTermStructure);
  //
  // create stochastic process
  Handle<Quote> initialFixingHandle(boost::shared_ptr<Quote>(new SimpleQuote(initialFixing)));
  Real volatility = 0.16;
  auto volatilityQuote = boost::make_shared<SimpleQuote>(volatility);
  Handle<Quote> volatilityHandle(volatilityQuote);
  Handle<BlackVolTermStructure> volatilityTermStructureHandle(boost::shared_ptr<BlackVolTermStructure>
   (new BlackConstantVol(settlementDays, calendar, volatilityHandle, dayCounter)));
  auto process = boost::make_shared<BlackScholesProcess>(initialFixingHandle, riskFreeRateTermStructureHandle, volatilityTermStructureHandle);
  //
  // create simulation-related attributes
  bool useAntitheticVariates = false;
  Size requiredSamples = 1000;
  Real requiredTolerance = Null<Real>();
  Size maxSamples = 1000;
  BigNatural seed = 0;
  //
  auto engine = boost::make_shared<EquityLinkedNoteMonteCarloPricer<PseudoRandom, Statistics>>
   (process, riskFreeRateTermStructureHandle, useAntitheticVariates, requiredSamples, requiredTolerance, maxSamples, seed);
  note->setPricingEngine(engine);
  std::cout << note->NPV() << std::endl;
 }
 catch (std::exception& e) {
  std::cout << e.what() << std::endl;
 }
 return 0;
}

Excel

The following screenshots are presenting parameters, result, path simulations and coupon calculations in Excel for this structured note. PV is an average of all discounted path payoffs.

QuantLib : custom instrument and pricing engine implementation

Using in-built QuantLib instruments and pricing engines is almost trivial, since the library is actually performing all required groundwork on behalf of a lucky user. Hair may start to get gray, as soon as one wants to implement a new instrument and related pricing engine. The example presented below is implementing a simple zero-coupon bond instrument and corresponding (analytical) pricing engine.

// MJZeroCouponBond.h
#pragma once
#include <ql/quantlib.hpp>
using namespace QuantLib;
//
// 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_;
};
//
//
//
//
//
// MJZeroCouponBond.cpp
#include "MJZeroCouponBond.h"
using namespace QuantLib;
//
// 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;
}
//
//
//
//
//
// Tester.cpp
#pragma once
#include "MJZeroCouponBond.h"
using namespace QuantLib;
//
int main() {
 //
 try {
  // common data : calendar, daycounter, dates
  Calendar calendar = TARGET();
  DayCounter dayCounter = Actual360();
  Date transactionDate(18, December, 2017);
  Date maturityDate(18, December, 2030);
  Natural settlementDays = 2;
  Date settlementDate = calendar.advance(transactionDate, Period(settlementDays, Days));
  Settings::instance().evaluationDate() = settlementDate;
  //
  // create flat discount curve
  auto riskFreeRate = boost::make_shared<SimpleQuote>(0.01);
  Handle<Quote> riskFreeRateHandle(riskFreeRate);
  auto riskFreeRateTermStructure = boost::make_shared<FlatForward>(settlementDate, riskFreeRateHandle, dayCounter);
  Handle<YieldTermStructure> riskFreeRateTermStructureHandle(riskFreeRateTermStructure);
  //
  // create zero coupon bonds (QL in-built and custom implementation)
  Real faceAmount = 250000.0;
  auto QLZero = boost::make_shared<ZeroCouponBond>(settlementDays, calendar, faceAmount, maturityDate);
  auto MJZero = boost::make_shared<MJZeroCouponBond>(faceAmount, maturityDate);
  //
  // create pricing engines  (QL in-built and custom implementation)
  auto QLEngine = boost::make_shared<DiscountingBondEngine>(riskFreeRateTermStructureHandle);
  auto MJEngine = boost::make_shared<MJZeroCouponBondEngine>(riskFreeRateTermStructureHandle);
  //
  // set engines to instruments and price the both transaction
  QLZero->setPricingEngine(QLEngine);
  MJZero->setPricingEngine(MJEngine);
  //
  std::cout << "pre-shift QL zero " << QLZero->NPV() << std::endl;
  std::cout << "pre-shift MJ zero " << MJZero->NPV() << std::endl;
  //
  // shift valuation curve and re-price the both transactions
  riskFreeRate->setValue(0.0101);
  std::cout << "post-shift QL zero " << QLZero->NPV() << std::endl;
  std::cout << "post-shift MJ zero " << MJZero->NPV() << std::endl;
 }
 catch (std::exception& e) {
  std::cout << e.what() << std::endl;
 }
 return 0;
}

In order to get more familiar with QuantLib Monte Carlo framework, one may take a look at Implementing QuantLib blog or get the actual book on QuantLib implementation from Leanpub. Also, MoneyScience is organizing Introduction to QuantLib Development course regularly. The most essential parts of the library are thoroughly covered during this 3-day training course, hosted by QuantLib lead developer Luigi Ballabio. Thanks for reading this blog. Merry Christmas and Happy New Year for everybody.
-Mike

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

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

C++11 : modelling one-factor processes using functional programming paradigm

There was an interesting technical article on July 2017 Wilmott magazine written by Daniel Duffy and Avi Palley. This multi-page article was giving an overview on some "game changing" Boost libraries, which have been accepted as a part of C++11/14 standard, such as smart pointers, function wrappers, lambda expressions and tuples. The second part of this article series will be published in the next Wilmott magazine and it will present (according to editor) design and implementation for Monte Carlo option pricing framework. It is also an important point of mentioning, that (according to Amazon.com) Daniel Duffy is about to publish long-awaited second edition of his book on pricing derivatives using C++. The book was initially published somewhere back in 2004 and the landscape has changed quite dramatically since these days.

Within the last chapter of this article, functional programming paradigm was nicely applied for modelling one-factor stochastic differential equations, generally used in Finance. By applying more functional programming paradigm, the usually observed code bloating can be substantially reduced. As an example of such code bloat, I reviewed my own implementation for path generator, which models one-factor processes for Monte Carlo purposes. There is an abstract base class (OneFactorProcess) and implementations for GBM and Vasicek processes. Even there is nothing fundamentally wrong with this approach (class hierarchy), one may ask, whether there would be a bit more flexible ways to implement this kind of a scheme.

Within this post, I have been re-designing modelling part for one-factor processes by applying functional programming paradigm, as presented in that article. Reduction in code bloating is present, since there is currently only one single class for modelling different types of processes (before, there was a class hierarchy). Moreover, since the both functions for handling drift and diffusion terms will be constructed outside of this class, their construction process is now much more flexible than before.

The program

 

Implement the following program (two header files and one implementation file for tester) into a new project. For brevity reasons, I have re-designed only the part of the program, which models one-factor processes. Monte Carlo part has been implemented as free function in tester implementation file. First, one-factor process object (Process) will be created by using ProcessBuilder object (Builder Pattern). Within this example, I have implemented a builder for constructing Process object by using console. However, the flexibility in this program allows different types of builders to be implemented. As soon as Process object is created, "skeleton path" (std::vector) will be sent to Monte Carlo method (MonteCarloLite), along with all simulation-related attributes and Process object. As a result, this method will fill vector with the values from one simulation path for chosen parameters and applied stochastic process. Finally, a path will be printed back to console.

#pragma once
#include <functional>
// OneFactorSDE.h
namespace FunctionalOneFactorProcessExampleNamespace
{
 // blueprint for a function, which models drift or diffusion
 using Function = std::function<double(double s, double t)>;
 // wrapper for drift and diffusion function components
 using Functions = std::tuple<Function, Function>;
 //
 // class for modeling one-factor processes by employing drift and diffusion functions
 class Process
 {
 public:
  Process(Functions functions, double initialCondition)
   : driftFunction(std::get<0>(functions)), diffusionFunction(std::get<1>(functions)),
   initialCondition(initialCondition) { }
  double drift(double s, double t) { return this->driftFunction(s, t); }
  double diffusion(double s, double t) { return this->diffusionFunction(s, t); }
  double InitialCondition() const { return this->initialCondition; }
  //
 private:
  double initialCondition; // spot price for a process
  Function driftFunction; // function for modelling process drift
  Function diffusionFunction;  // function for modelling process diffusion
 };
}
//
//
//
// ProcessFactory.h
#pragma once
#include <iostream>
#include <string>
#include <memory>
#include "OneFactorSDE.h"
//
namespace FunctionalOneFactorProcessExampleNamespace
{
 // abstract base class for all process builders
 class ProcessBuilder
 {
 public:
  // return process object which is wrapped inside shared pointer
  // let implementation classes decide, how and from where to built process object
  virtual std::shared_ptr<Process> Build() = 0;
 };
 //
 // specific implementation for console process builder
 // process type and corresponding parameters will be requested from a client by using console
 class ConsoleProcessBuilder : public ProcessBuilder
 {
 public:
  std::shared_ptr<Process> Build() override
  {
   Functions functions;
   double initialCondition = 0.0;
   std::string consoleSelection;
   std::cout << "Select process [1 = GBM, 2 = Vasicek] > ";
   // if conversion cannot be performed, stoi will throw invalid argument exception
   std::getline(std::cin, consoleSelection);
   int processID = std::stoi(consoleSelection);
   //
   switch (processID)
   {
   // GBM process
   case 1:
   {
    // receive client inputs
    std::cout << "spot price > ";
    std::getline(std::cin, consoleSelection);
    initialCondition = std::stod(consoleSelection);
    //
    std::cout << "risk-free rate > ";
    std::getline(std::cin, consoleSelection);
    double r = std::stod(consoleSelection);
    //
    std::cout << "volatility > ";
    std::getline(std::cin, consoleSelection);
    double v = std::stod(consoleSelection);
    //
    // build drift and diffusion functions for GBM process object
    auto driftFunction = [r](double S, double T){ return r * S; };
    auto diffusionFunction = [v](double S, double T){ return v * S; };
    // wrap drift and diffusion functions into tuple
    functions = std::make_tuple(driftFunction, diffusionFunction);
    break;
   }
   case 2:
   {
    // receive client inputs
    std::cout << "spot price > ";
    std::getline(std::cin, consoleSelection);
    initialCondition = std::stod(consoleSelection);
    //
    std::cout << "reversion > ";
    std::getline(std::cin, consoleSelection);
    double reversion = std::stod(consoleSelection);
    //
    std::cout << "long-term rate > ";
    std::getline(std::cin, consoleSelection);
    double longTermRate = std::stod(consoleSelection);
    //
    std::cout << "rate volatility > ";
    std::getline(std::cin, consoleSelection);
    double v = std::stod(consoleSelection);
    //
    // build drift and diffusion functions for Vasicek process object
    auto driftFunction = [reversion, longTermRate](double S, double T)
     { return reversion * (longTermRate - S); };
    auto diffusionFunction = [v](double S, double T){ return v; };
    // wrap drift and diffusion functions into tuple
    functions = std::make_tuple(driftFunction, diffusionFunction);
    break;
   }
   default:
    // if selected process is not configured, program will throw invalid argument exception
    throw std::invalid_argument("invalid process ID");
    break;
   }
   // build and return constructed process object for a client
   // wrapped into shared pointer
   return std::shared_ptr<Process>(new Process(functions, initialCondition));
  }
 };
}
//
//
//
// Tester.cpp
#include <vector>
#include <random>
#include <algorithm>
#include <chrono>
#include "ProcessFactory.h"
namespace MJ = FunctionalOneFactorProcessExampleNamespace;
//
void MonteCarloLite(std::vector<double>& path, 
 std::shared_ptr<MJ::Process>& process, const double maturity);
//
int main()
{
 try
 {
  // create process object by using console builder
  MJ::ConsoleProcessBuilder builder;
  std::shared_ptr<MJ::Process> process = builder.Build();
  //
  // create simulation-related attributes
  int nSteps = 100;
  double timeToMaturity = 1.25;
  // create, process and print one path
  std::vector<double> path(nSteps);
  MonteCarloLite(path, process, timeToMaturity);
  std::for_each(path.begin(), path.end(), 
   [](double v) { std::cout << v << std::endl; });
 }
 catch (std::exception e)
 {
  std::cout << e.what() << std::endl;
 }
 return 0;
}
//
void MonteCarloLite(std::vector<double>& path, 
 std::shared_ptr<MJ::Process>& process, const double maturity)
{
 // lambda method for seeding uniform random generator
 std::function<unsigned long(void)> seeder =
  [](void) -> unsigned long { return static_cast<unsigned long>
  (std::chrono::steady_clock::now().time_since_epoch().count()); };
 //
 // create uniform generator, distribution and random generator function
 std::mt19937 uniformGenerator(seeder());
 std::normal_distribution<double> distribution;
 std::function<double(double)> randomGenerator;
 //
 // lambda method for processing standard normal random numbers
 randomGenerator = [uniformGenerator, distribution](double x) mutable -> double
 {
  x = distribution(uniformGenerator);
  return x;
 };
 //
 // create vector of standard normal random numbers
 // use lambda method created earlier
 std::transform(path.begin(), path.end(), path.begin(), randomGenerator);
 //
 double dt = maturity / (path.size() - 1);
 double dw = 0.0;
 double s = (*process).InitialCondition();
 double t = 0.0;
 path[0] = s; // 1st path element is always the current spot price
 //
 // transform random number vector into a path containing asset prices
 for (auto it = path.begin() + 1; it != path.end(); ++it)
 {
  t += dt;
  dw = (*it) * std::sqrt(dt);
  (*it) = s + (*process).drift(s, t) * dt + (*process).diffusion(s, t) * dw;
  s = (*it);
 }
}

As always, thanks a lot for reading this blog.
-Mike