QuantLib : Simulating HW1F paths using PathGenerator

Monte Carlo is bread and butter for so many purposes. Calculating payoffs for complex path-dependent products or simulating future exposures for calculating CVA are two excellent examples. The big question is always how to do this efficiently. Designing, implementing and setting up any non-trivial in-house tool to do the job is everything but not a simple afternoon exercise with a cup of coffee and Excel. Fortunately, QuantLib is offering pretty impressive tools for simulating stochastic paths. This time, I wanted to share the results of my woodshedding with QL PathGenerator class. 

Parallel lives

In order to really appreciate the tools offered by QL, let us see the results first. Some simulated paths using Hull-White One-Factor model are shown in the picture below.

If one really want to start from the scratch, there are a lot of things to do in order to produce these paths on a flexible manner and handling all the complexities of the task at the same time. Thanks for QL, those days are finally over.

Legoland

 

Setting up desired Stochastic Process and Gaussian Sequence Generator are two main components needed in order to get this thing up and running.

Along with required process parameters (reversion speed and rate volatility), HullWhiteProcess needs Handle to YieldTermStructure object, such as PiecewiseYieldCurve, as an input.

 // create Hull-White one-factor stochastic process
 Real reversionSpeed = 0.75;
 Real rateVolatility = 0.015;
 boost::shared_ptr<StochasticProcess1D> HW1F(
  new HullWhiteProcess(curveHandle, reversionSpeed, rateVolatility));

For this example, I have used my own PiecewiseCurveBuilder template class in order to make curve assembling a bit more easier. It should be noted, that the menu of one-dimensional stochastic processes in QL is covering pretty much all standard processes one needs for different asset classes.

Gaussian Sequence Generator (GSG) is assembled by using the following three classes : uniform random generator (MersenneTwisterUniformRng),  distributional transformer (CLGaussianRng) and RandomSequenceGenerator.

 // type definition for complex declaration
 typedef RandomSequenceGenerator<CLGaussianRng<MersenneTwisterUniformRng>> GSG;
 //
 // create mersenne twister uniform random generator
 unsigned long seed = 28749;
 MersenneTwisterUniformRng generator(seed);
 //
 // create gaussian generator by using central limit transformation method
 CLGaussianRng<MersenneTwisterUniformRng> gaussianGenerator(generator);
 //
 // define maturity, number of steps per path and create gaussian sequence generator
 Time maturity = 5.0;
 Size nSteps = 1250;
 GSG gaussianSequenceGenerator(nSteps, gaussianGenerator);
 //
 // create path generator using Hull-White process and gaussian sequence generator
 PathGenerator<GSG> pathGenerator(HW1F, maturity, nSteps, gaussianSequenceGenerator, false);

Finally, PathGenerator object is created by feeding desired process and generator objects in constructor method, along with the other required parameters (maturity, number of steps). After this, PathGenerator object is ready for producing stochastic paths for its client.

The program

 

Example program will first create relinkable handle to PiecewiseYieldCurve object. Remember to include required files into your project from here. After this, the program creates HW1F process object and Gaussian Sequence Generator object, which are feeded into PathGenerator object. Finally, the program creates 20 stochastic paths, which are saved into Matrix object and ultimately being printed into text file for further analysis (Excel chart).

#include "PiecewiseCurveBuilder.cpp"
#include <fstream>
#include <string>
//
// type definition for complex declaration
typedef RandomSequenceGenerator<CLGaussianRng<MersenneTwisterUniformRng>> GSG;
//
// function prototypes
RelinkableHandle<YieldTermStructure> CreateCurveHandle(Date settlementDate);
void PrintMatrix(const Matrix& matrix, std::string filePathName);
//
int main()
{
 // request handle for piecewise USD Libor curve
 Date tradeDate(22, January, 2016);
 Settings::instance().evaluationDate() = tradeDate;
 Date settlementDate = UnitedKingdom().advance(tradeDate, 2, Days);
 RelinkableHandle<YieldTermStructure> curveHandle = CreateCurveHandle(settlementDate);
 //
 // create Hull-White one-factor stochastic process
 Real reversionSpeed = 0.75;
 Real rateVolatility = 0.015;
 boost::shared_ptr<StochasticProcess1D> HW1F(
  new HullWhiteProcess(curveHandle, reversionSpeed, rateVolatility));
 //
 // create mersenne twister uniform random generator
 unsigned long seed = 28749;
 MersenneTwisterUniformRng generator(seed);
 //
 // create gaussian generator by using central limit transformation method
 CLGaussianRng<MersenneTwisterUniformRng> gaussianGenerator(generator);
 //
 // define maturity, number of steps per path and create gaussian sequence generator
 Time maturity = 5.0;
 Size nSteps = 1250;
 GSG gaussianSequenceGenerator(nSteps, gaussianGenerator);
 //
 // create path generator using Hull-White process and gaussian sequence generator
 PathGenerator<GSG> pathGenerator(HW1F, maturity, nSteps, gaussianSequenceGenerator, false);
 //
 // create matrix container for 20 generated paths
 Size nColumns = 20;
 Matrix paths(nSteps + 1, nColumns);
 for(unsigned int i = 0; i != paths.columns(); i++)
 {
  // request a new stochastic path from path generator
  QuantLib::Sample<Path> path = pathGenerator.next();
  //
  // save generated path into container
  for(unsigned int j = 0; j != path.value.length(); j++)
  {
   paths[j][i] = path.value.at(j);
  }
 }
 // finally, print matrix content into text file
 PrintMatrix(paths, "C:\\temp\\HW1F.txt");
 return 0;
}
//
void PrintMatrix(const Matrix& matrix, std::string filePathName)
{
 // open text file for input, loop through matrix rows
 std::ofstream file(filePathName);
 for(unsigned int i = 0; i != matrix.rows(); i++)
 {
  // concatenate column values into string separated by semicolon
  std::string stream;
  for(unsigned int j = 0; j != matrix.columns(); j++)
  {
   stream += (std::to_string(matrix[i][j]) + ";");
  }
  // print string into text file
  file << stream << std::endl;
 }
 // close text file
 file.close();
}
//
RelinkableHandle<YieldTermStructure> CreateCurveHandle(Date settlementDate)
{
 // create curve builder for piecewise USD Libor swap curve
 PiecewiseCurveBuilder<USDLibor> USDCurveBuilder(settlementDate, 
  UnitedKingdom(), Annual, Thirty360());
 //
 // add quotes directly into curve builder
 USDCurveBuilder.AddDeposit(0.0038975, 1 * Weeks);
 USDCurveBuilder.AddDeposit(0.004295, 1 * Months);
 USDCurveBuilder.AddDeposit(0.005149, 2 * Months);
 USDCurveBuilder.AddDeposit(0.006127, 3 * Months);
 USDCurveBuilder.AddFRA(0.008253, 3 * Months, 3 * Months);
 USDCurveBuilder.AddFRA(0.009065, 6 * Months, 3 * Months);
 USDCurveBuilder.AddFRA(0.01059, 9 * Months, 3 * Months);
 USDCurveBuilder.AddSwap(0.011459, 2 * Years);
 USDCurveBuilder.AddSwap(0.013745, 3 * Years);
 USDCurveBuilder.AddSwap(0.015475, 4 * Years);
 USDCurveBuilder.AddSwap(0.016895, 5 * Years);
 USDCurveBuilder.AddSwap(0.01813, 6 * Years);
 USDCurveBuilder.AddSwap(0.019195, 7 * Years);
 USDCurveBuilder.AddSwap(0.020115, 8 * Years);
 USDCurveBuilder.AddSwap(0.020905, 9 * Years);
 USDCurveBuilder.AddSwap(0.021595, 10 * Years);
 USDCurveBuilder.AddSwap(0.0222, 11 * Years);
 USDCurveBuilder.AddSwap(0.022766, 12 * Years);
 USDCurveBuilder.AddSwap(0.0239675, 15 * Years);
 USDCurveBuilder.AddSwap(0.025105, 20 * Years);
 USDCurveBuilder.AddSwap(0.025675, 25 * Years);
 USDCurveBuilder.AddSwap(0.026015, 30 * Years);
 USDCurveBuilder.AddSwap(0.026205, 40 * Years);
 USDCurveBuilder.AddSwap(0.026045, 50 * Years);
 //
 // return relinkable curve handle
 return USDCurveBuilder.GetCurveHandle();
}

Thanks for reading my blog.

-Mike

QuantLib : Builder Class for PiecewiseYieldCurve

Selection of high-quality benchmark securities and bootstrapping of valuation curve is the bread and butter in valuing financial transactions. In one of my blog, I opened up one possible framework for this procedure. Needless to say, that program was still far away from being really easy-to-use and state-of-the-art implementation. For this reason, I wanted to take a look at some external analytics libraries I have been getting introduced last year.

This time, I wanted to share some fruits of woodshedding QuantLib way to build up piecewise term structure. As I was initially tackling through the example files found in QL downloadable package and the other relevant examples found with Google, there was one obvious issue : with all those rates, quotes and handles involved, a lot of administrative code writing needs to be done first, in order to get any piecewise curve up and running. Since we need valuation curves anyway for everything we do with QL, there should be a manageable way to handle this issue.

One template to rule them all

For the reason mentioned above, I came up with an idea of auxiliary piecewise curve builder class. The purpose of this class would be simple : it could give a client a chance to assemble piecewise curve by adding arbitrary amount of different types of quotes (deposit, FRA or swap) and finally a client could request handle for assembled curve. The resulting curve handle could then be directly used in other QL programs. Effectively, this class would be wrapping some of the required administrative code work away from a client.

Walkthrough

As an example of the usage of this PiecewiseCurveBuilder template class, let us go through some parts of the example program. In order to keep this example as simple as possible, we will operate only with one rate for each category (deposit, FRA, swap). First, we create stand-alone quotes from market data. Memory for a SimpleQuote (inherits from Quote base class) will be reserved by using Boost shared pointer :

  boost::shared_ptr<Quote> usd_quote_deposit_3M(new SimpleQuote(0.006127));
  boost::shared_ptr<Quote> usd_quote_fra_3M6M(new SimpleQuote(0.008253));
  boost::shared_ptr<Quote> usd_quote_swap_2Y(new SimpleQuote(0.011459));

Next, we create curve builder class instance for assembling USD Libor swap curve :

  Date settlementDate = UnitedKingdom().advance(tradeDate, 2, Days);
  PiecewiseCurveBuilder<USDLibor> USDCurveBuilder(settlementDate, UnitedKingdom(), Annual, Thirty360());

After this, we add quotes into USD curve builder :

  USDCurveBuilder.AddDeposit(usd_quote_deposit_3M, 3 * Months);
  USDCurveBuilder.AddFRA(usd_quote_fra_3M6M, 3 * Months, 3 * Months);
  USDCurveBuilder.AddSwap(usd_quote_swap_2Y, 2 * Years);

Finally, we request relinkable curve handle from USD curve builder :

  RelinkableHandle<YieldTermStructure> curveHandle = USDCurveBuilder.GetCurveHandle();
  DoSomethingWithCurveHandle(curveHandle);

Requested RelinkableHandle (inherits from Handle base class) can then be directly used in other QL programs. One should be aware, that all the changes we made in stand-alone quotes (SimpleQuotes) will have an effect on the curve. For an example, we could modify the existing quotes by shocking their rate up by 100 bps :

  boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_deposit_3M)->setValue(usd_quote_deposit_3M->value() + 0.01);
  boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_fra_3M6M)->setValue(usd_quote_fra_3M6M->value() + 0.01);
  boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_swap_2Y)->setValue(usd_quote_swap_2Y->value() + 0.01);
  DoSomethingWithCurveHandle(curveHandle);

Additional add methods for different quote types are enabling even easier way to assemble a curve. Instead of giving a rate wrapped inside Quote object, it can be given directly into curve builder. First, we create curve builder for CHF Libor swap curve :

  settlementDate = UnitedKingdom().advance(tradeDate, 2, Days);
  PiecewiseCurveBuilder<CHFLibor> CHFCurveBuilder(settlementDate, UnitedKingdom(), Annual, Thirty360());

Next, we add market rates directly into CHF curve builder :

  CHFCurveBuilder.AddDeposit(-0.006896, 6 * Months);
  CHFCurveBuilder.AddFRA(-0.007103, 6 * Months, 6 * Months);
  CHFCurveBuilder.AddSwap(-0.0068355, 2 * Years);

Finally, we request relinkable curve handle from CHF curve builder :

  curveHandle = CHFCurveBuilder.GetCurveHandle();
  DoSomethingWithCurveHandle(curveHandle);

This last option would be suitable in situations, where a client has no need for any auxiliary rate
updating.

The complete example program has been presented below. The program will first create updateable USD Libor swap curve, print all the rates, modify quotes and re-prints the rates. After this, the program will create updateable CHF Libor swap curve and prints a set of discount factors. Finally, the program will create non-updateable CHF Libor swap curve by using another set of add methods and prints a set of discount factors.

Have a great start for the year 2016. Thanks for reading my blog.
-Mike

The program

// PiecewiseCurveBuilder.h
#pragma once
#include <ql/quantlib.hpp>
using namespace QuantLib;
//
template <class T>
class PiecewiseCurveBuilder
{
private:
 Date settlementDate;
 Calendar calendar;
 Frequency fixedLegFrequency;
 DayCounter dayCounter;
 DayCounter fixedLegDayCounter;
 BusinessDayConvention businessDayConvention;
 std::vector<boost::shared_ptr<RateHelper>> rateHelpers;
 boost::shared_ptr<YieldTermStructure> yieldTermStructure;
public:
 PiecewiseCurveBuilder(Date settlementDate, Calendar calendar, 
  Frequency fixedLegFrequency, DayCounter fixedLegDayCounter);
 void AddDeposit(boost::shared_ptr<Quote>& quote, Period periodLength);
 void AddDeposit(Rate quote, Period periodLength);
 void AddFRA(boost::shared_ptr<Quote>& quote, Period periodLengthToStart, Period periodLength);
 void AddFRA(Rate quote, Period periodLengthToStart, Period periodLength);
 void AddSwap(boost::shared_ptr<Quote>& quote, Period periodLength);
 void AddSwap(Rate quote, Period periodLength);
 RelinkableHandle<YieldTermStructure> GetCurveHandle();
};
//
//
//
// PiecewiseCurveBuilder.cpp
#pragma once
#include "PiecewiseCurveBuilder.h"
//
template <class T>
PiecewiseCurveBuilder<T>::PiecewiseCurveBuilder(Date settlementDate, 
 Calendar calendar, Frequency fixedLegFrequency, DayCounter fixedLegDayCounter)
{
 this->settlementDate = settlementDate;
 this->calendar = calendar;
 this->fixedLegFrequency = fixedLegFrequency;
 this->fixedLegDayCounter = fixedLegDayCounter;
 boost::shared_ptr<IborIndex> index(new T(3 * Months));
 dayCounter = index->dayCounter();
 businessDayConvention = index->businessDayConvention();
}
template <class T>
void PiecewiseCurveBuilder<T>::AddDeposit(boost::shared_ptr<Quote>& quote, Period periodLength)
{
 // using third DepositRateHelper constructor
 boost::shared_ptr<RateHelper> rateHelper(new DepositRateHelper(
  Handle<Quote>(quote), boost::shared_ptr<IborIndex>(new T(periodLength))));
 rateHelpers.push_back(rateHelper);
}
template <class T>
void PiecewiseCurveBuilder<T>::AddDeposit(Rate quote, Period periodLength)
{
 // using second DepositRateHelper constructor
 boost::shared_ptr<RateHelper> rateHelper(new DepositRateHelper(
  quote, boost::shared_ptr<IborIndex>(new T(periodLength))));
 rateHelpers.push_back(rateHelper);
}
template <class T>
void PiecewiseCurveBuilder<T>::AddFRA(boost::shared_ptr<Quote>& quote, 
 Period periodLengthToStart, Period periodLength)
{
 // using seventh FraRateHelper constructor
 boost::shared_ptr<RateHelper> rateHelper(new FraRateHelper(
  Handle<Quote>(quote), periodLengthToStart, 
  boost::shared_ptr<IborIndex>(new T(periodLength))));
 rateHelpers.push_back(rateHelper); 
}
template <class T>
void PiecewiseCurveBuilder<T>::AddFRA(Rate quote, 
 Period periodLengthToStart, Period periodLength)
{
 // using third FraRateHelper constructor
 boost::shared_ptr<RateHelper> rateHelper(new FraRateHelper(
  quote, periodLengthToStart, 
  boost::shared_ptr<IborIndex>(new T(periodLength))));
 rateHelpers.push_back(rateHelper); 
}
template <class T>
void PiecewiseCurveBuilder<T>::AddSwap(boost::shared_ptr<Quote>& quote, 
 Period periodLength)
{
 // using fifth SwapRateHelper constructor
 boost::shared_ptr<IborIndex> index(new T(periodLength));
 boost::shared_ptr<RateHelper> rateHelper(new SwapRateHelper(
  Handle<Quote>(quote), periodLength, calendar, fixedLegFrequency, 
  businessDayConvention, fixedLegDayCounter, index));
 rateHelpers.push_back(rateHelper);
}
template <class T>
void PiecewiseCurveBuilder<T>::AddSwap(Rate quote, 
 Period periodLength)
{
 // using fourth SwapRateHelper constructor
 boost::shared_ptr<IborIndex> index(new T(periodLength));
 boost::shared_ptr<RateHelper> rateHelper(new SwapRateHelper(
  quote, periodLength, calendar, fixedLegFrequency, 
  businessDayConvention, fixedLegDayCounter, index));
 rateHelpers.push_back(rateHelper);
}
template <class T>
RelinkableHandle<YieldTermStructure> PiecewiseCurveBuilder<T>::GetCurveHandle()
{
 if(yieldTermStructure == NULL)
 {
  yieldTermStructure = boost::shared_ptr<YieldTermStructure>(
   new PiecewiseYieldCurve<Discount, LogLinear>(
   settlementDate, rateHelpers, dayCounter));
 }
 return RelinkableHandle<YieldTermStructure>(yieldTermStructure);
}
//
//
//
// Tester.cpp
#include "PiecewiseCurveBuilder.cpp"
//
// prototype : template method for calculating a fair swap rate
template <typename T> Rate GetSwapRate(Period swapMaturity, Date settlementDate, Period floatingLegTenor, 
 Handle<YieldTermStructure>& curveHandle, Calendar calendar, DayCounter fixedLegDayCount, Period fixedLegTenor);
// prototype : hard-coded printer for USD rates
void PrintUSDRates(Date settlementDate, Handle<YieldTermStructure>& curveHandle);
// prototype : hard-coded printer for CHF discount factors
void PrintCHFDiscountFactors(Date settlementDate, Handle<YieldTermStructure>& curveHandle);
//
int main()
{
 // create trade date
 Date tradeDate(5, January, 2016);
 Settings::instance().evaluationDate() = tradeDate;
 //
 // create relinkable handle for curve
 RelinkableHandle<YieldTermStructure> curveHandle;
 //
 // create stand-alone quotes from USD market data
 boost::shared_ptr<Quote> usd_quote_deposit_1W(new SimpleQuote(0.0038975));
 boost::shared_ptr<Quote> usd_quote_deposit_1M(new SimpleQuote(0.004295));
 boost::shared_ptr<Quote> usd_quote_deposit_2M(new SimpleQuote(0.005149));
 boost::shared_ptr<Quote> usd_quote_deposit_3M(new SimpleQuote(0.006127));
 boost::shared_ptr<Quote> usd_quote_fra_3M6M(new SimpleQuote(0.008253));
 boost::shared_ptr<Quote> usd_quote_fra_6M9M(new SimpleQuote(0.009065));
 boost::shared_ptr<Quote> usd_quote_fra_9M12M(new SimpleQuote(0.01059));
 boost::shared_ptr<Quote> usd_quote_swap_2Y(new SimpleQuote(0.011459));
 boost::shared_ptr<Quote> usd_quote_swap_3Y(new SimpleQuote(0.013745));
 boost::shared_ptr<Quote> usd_quote_swap_4Y(new SimpleQuote(0.015475));
 boost::shared_ptr<Quote> usd_quote_swap_5Y(new SimpleQuote(0.016895));
 boost::shared_ptr<Quote> usd_quote_swap_6Y(new SimpleQuote(0.01813));
 boost::shared_ptr<Quote> usd_quote_swap_7Y(new SimpleQuote(0.019195));
 boost::shared_ptr<Quote> usd_quote_swap_8Y(new SimpleQuote(0.020115));
 boost::shared_ptr<Quote> usd_quote_swap_9Y(new SimpleQuote(0.020905));
 boost::shared_ptr<Quote> usd_quote_swap_10Y(new SimpleQuote(0.021595));
 boost::shared_ptr<Quote> usd_quote_swap_11Y(new SimpleQuote(0.0222));
 boost::shared_ptr<Quote> usd_quote_swap_12Y(new SimpleQuote(0.022766));
 boost::shared_ptr<Quote> usd_quote_swap_15Y(new SimpleQuote(0.0239675));
 boost::shared_ptr<Quote> usd_quote_swap_20Y(new SimpleQuote(0.025105));
 boost::shared_ptr<Quote> usd_quote_swap_25Y(new SimpleQuote(0.025675));
 boost::shared_ptr<Quote> usd_quote_swap_30Y(new SimpleQuote(0.026015));
 boost::shared_ptr<Quote> usd_quote_swap_40Y(new SimpleQuote(0.026205));
 boost::shared_ptr<Quote> usd_quote_swap_50Y(new SimpleQuote(0.026045));
 //
 // create curve builder for USD Libor swap curve
 Date settlementDate = UnitedKingdom().advance(tradeDate, 2, Days);
 PiecewiseCurveBuilder<USDLibor> USDCurveBuilder(settlementDate, 
  UnitedKingdom(), Annual, Thirty360());
 //
 // add quotes (reference to shared pointer) into USD curve builder
 USDCurveBuilder.AddDeposit(usd_quote_deposit_1W, 1 * Weeks);
 USDCurveBuilder.AddDeposit(usd_quote_deposit_1M, 1 * Months);
 USDCurveBuilder.AddDeposit(usd_quote_deposit_2M, 2 * Months);
 USDCurveBuilder.AddDeposit(usd_quote_deposit_3M, 3 * Months);
 USDCurveBuilder.AddFRA(usd_quote_fra_3M6M, 3 * Months, 3 * Months);
 USDCurveBuilder.AddFRA(usd_quote_fra_6M9M, 6 * Months, 3 * Months);
 USDCurveBuilder.AddFRA(usd_quote_fra_9M12M, 9 * Months, 3 * Months);
 USDCurveBuilder.AddSwap(usd_quote_swap_2Y, 2 * Years);
 USDCurveBuilder.AddSwap(usd_quote_swap_3Y, 3 * Years);
 USDCurveBuilder.AddSwap(usd_quote_swap_4Y, 4 * Years);
 USDCurveBuilder.AddSwap(usd_quote_swap_5Y, 5 * Years);
 USDCurveBuilder.AddSwap(usd_quote_swap_6Y, 6 * Years);
 USDCurveBuilder.AddSwap(usd_quote_swap_7Y, 7 * Years);
 USDCurveBuilder.AddSwap(usd_quote_swap_8Y, 8 * Years);
 USDCurveBuilder.AddSwap(usd_quote_swap_9Y, 9 * Years);
 USDCurveBuilder.AddSwap(usd_quote_swap_10Y, 10 * Years);
 USDCurveBuilder.AddSwap(usd_quote_swap_11Y, 11 * Years);
 USDCurveBuilder.AddSwap(usd_quote_swap_12Y, 12 * Years);
 USDCurveBuilder.AddSwap(usd_quote_swap_15Y, 15 * Years);
 USDCurveBuilder.AddSwap(usd_quote_swap_20Y, 20 * Years);
 USDCurveBuilder.AddSwap(usd_quote_swap_25Y, 25 * Years);
 USDCurveBuilder.AddSwap(usd_quote_swap_30Y, 30 * Years);
 USDCurveBuilder.AddSwap(usd_quote_swap_40Y, 40 * Years);
 USDCurveBuilder.AddSwap(usd_quote_swap_50Y, 50 * Years);
 //
 // get relinkable curve handle from USD curve builder and print rates
 curveHandle = USDCurveBuilder.GetCurveHandle();
 PrintUSDRates(settlementDate, curveHandle);
 //
 // modify existing USD quotes by shocking rates up by 100 bps
 Real bump = 0.01;
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_deposit_1W)->setValue(usd_quote_deposit_1W->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_deposit_1M)->setValue(usd_quote_deposit_1M->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_deposit_2M)->setValue(usd_quote_deposit_2M->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_deposit_3M)->setValue(usd_quote_deposit_3M->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_fra_3M6M)->setValue(usd_quote_fra_3M6M->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_fra_6M9M)->setValue(usd_quote_fra_6M9M->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_fra_9M12M)->setValue(usd_quote_fra_9M12M->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_swap_2Y)->setValue(usd_quote_swap_2Y->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_swap_3Y)->setValue(usd_quote_swap_3Y->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_swap_4Y)->setValue(usd_quote_swap_4Y->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_swap_5Y)->setValue(usd_quote_swap_5Y->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_swap_6Y)->setValue(usd_quote_swap_6Y->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_swap_7Y)->setValue(usd_quote_swap_7Y->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_swap_8Y)->setValue(usd_quote_swap_8Y->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_swap_9Y)->setValue(usd_quote_swap_9Y->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_swap_10Y)->setValue(usd_quote_swap_10Y->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_swap_11Y)->setValue(usd_quote_swap_11Y->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_swap_12Y)->setValue(usd_quote_swap_12Y->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_swap_15Y)->setValue(usd_quote_swap_15Y->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_swap_20Y)->setValue(usd_quote_swap_20Y->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_swap_25Y)->setValue(usd_quote_swap_25Y->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_swap_30Y)->setValue(usd_quote_swap_30Y->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_swap_40Y)->setValue(usd_quote_swap_40Y->value() + bump);
 boost::dynamic_pointer_cast<SimpleQuote>(usd_quote_swap_50Y)->setValue(usd_quote_swap_50Y->value() + bump);
 //
 // re-print USD rates
 PrintUSDRates(settlementDate, curveHandle);
 //
 //
 //
 // create stand-alone quotes from CHF market data
 boost::shared_ptr<Quote> chf_quote_deposit_1W(new SimpleQuote(-0.00827));
 boost::shared_ptr<Quote> chf_quote_deposit_1M(new SimpleQuote(-0.00798));
 boost::shared_ptr<Quote> chf_quote_deposit_2M(new SimpleQuote(-0.00785));
 boost::shared_ptr<Quote> chf_quote_deposit_3M(new SimpleQuote(-0.00756));
 boost::shared_ptr<Quote> chf_quote_deposit_6M(new SimpleQuote(-0.006896));
 boost::shared_ptr<Quote> chf_quote_fra_6M12M(new SimpleQuote(-0.007103));
 boost::shared_ptr<Quote> chf_quote_swap_2Y(new SimpleQuote(-0.0068355));
 boost::shared_ptr<Quote> chf_quote_swap_3Y(new SimpleQuote(-0.006125));
 boost::shared_ptr<Quote> chf_quote_swap_4Y(new SimpleQuote(-0.0050195));
 boost::shared_ptr<Quote> chf_quote_swap_5Y(new SimpleQuote(-0.003725));
 boost::shared_ptr<Quote> chf_quote_swap_6Y(new SimpleQuote(-0.002425));
 boost::shared_ptr<Quote> chf_quote_swap_7Y(new SimpleQuote(-0.0011885));
 boost::shared_ptr<Quote> chf_quote_swap_8Y(new SimpleQuote(-0.000094));
 boost::shared_ptr<Quote> chf_quote_swap_9Y(new SimpleQuote(0.0008755));
 boost::shared_ptr<Quote> chf_quote_swap_10Y(new SimpleQuote(0.0016365));
 boost::shared_ptr<Quote> chf_quote_swap_12Y(new SimpleQuote(0.003));
 boost::shared_ptr<Quote> chf_quote_swap_15Y(new SimpleQuote(0.004525));
 boost::shared_ptr<Quote> chf_quote_swap_20Y(new SimpleQuote(0.0063));
 boost::shared_ptr<Quote> chf_quote_swap_25Y(new SimpleQuote(0.00735));
 boost::shared_ptr<Quote> chf_quote_swap_30Y(new SimpleQuote(0.007825));
 //
 // create curve builder for CHF Libor swap curve
 settlementDate = UnitedKingdom().advance(tradeDate, 2, Days);
 PiecewiseCurveBuilder<CHFLibor> CHFCurveBuilder(settlementDate, 
  UnitedKingdom(), Annual, Thirty360());
 //
 // add quotes (reference to shared pointer) into CHF curve builder
 CHFCurveBuilder.AddDeposit(chf_quote_deposit_1W, 1 * Weeks);
 CHFCurveBuilder.AddDeposit(chf_quote_deposit_1M, 1 * Months);
 CHFCurveBuilder.AddDeposit(chf_quote_deposit_2M, 2 * Months);
 CHFCurveBuilder.AddDeposit(chf_quote_deposit_3M, 3 * Months);
 CHFCurveBuilder.AddDeposit(chf_quote_deposit_6M, 6 * Months);
 CHFCurveBuilder.AddFRA(chf_quote_fra_6M12M, 6 * Months, 6 * Months);
 CHFCurveBuilder.AddSwap(chf_quote_swap_2Y, 2 * Years);
 CHFCurveBuilder.AddSwap(chf_quote_swap_3Y, 3 * Years);
 CHFCurveBuilder.AddSwap(chf_quote_swap_4Y, 4 * Years);
 CHFCurveBuilder.AddSwap(chf_quote_swap_5Y, 5 * Years);
 CHFCurveBuilder.AddSwap(chf_quote_swap_6Y, 6 * Years);
 CHFCurveBuilder.AddSwap(chf_quote_swap_7Y, 7 * Years);
 CHFCurveBuilder.AddSwap(chf_quote_swap_8Y, 8 * Years);
 CHFCurveBuilder.AddSwap(chf_quote_swap_9Y, 9 * Years);
 CHFCurveBuilder.AddSwap(chf_quote_swap_10Y, 10 * Years);
 CHFCurveBuilder.AddSwap(chf_quote_swap_12Y, 12 * Years);
 CHFCurveBuilder.AddSwap(chf_quote_swap_15Y, 15 * Years);
 CHFCurveBuilder.AddSwap(chf_quote_swap_20Y, 20 * Years);
 CHFCurveBuilder.AddSwap(chf_quote_swap_25Y, 25 * Years);
 CHFCurveBuilder.AddSwap(chf_quote_swap_30Y, 30 * Years);
 //
 // get relinkable curve handle from CHF curve builder and print discount factors
 curveHandle = CHFCurveBuilder.GetCurveHandle();
 PrintCHFDiscountFactors(settlementDate, curveHandle);
 //
 //
 //
 // create another curve builder for CHF Libor swap curve
 settlementDate = UnitedKingdom().advance(tradeDate, 2, Days);
 PiecewiseCurveBuilder<CHFLibor> CHFCurveBuilder2(settlementDate, 
  UnitedKingdom(), Annual, Thirty360());
 //
 // add market rates directly into CHF curve builder
 CHFCurveBuilder2.AddDeposit(-0.00827, 1 * Weeks);
 CHFCurveBuilder2.AddDeposit(-0.00798, 1 * Months);
 CHFCurveBuilder2.AddDeposit(-0.00785, 2 * Months);
 CHFCurveBuilder2.AddDeposit(-0.00756, 3 * Months);
 CHFCurveBuilder2.AddDeposit(-0.006896, 6 * Months);
 CHFCurveBuilder2.AddFRA(-0.007103, 6 * Months, 6 * Months);
 CHFCurveBuilder2.AddSwap(-0.0068355, 2 * Years);
 CHFCurveBuilder2.AddSwap(-0.006125, 3 * Years);
 CHFCurveBuilder2.AddSwap(-0.0050195, 4 * Years);
 CHFCurveBuilder2.AddSwap(-0.003725, 5 * Years);
 CHFCurveBuilder2.AddSwap(-0.002425, 6 * Years);
 CHFCurveBuilder2.AddSwap(-0.0011885, 7 * Years);
 CHFCurveBuilder2.AddSwap(-0.000094, 8 * Years);
 CHFCurveBuilder2.AddSwap(0.0008755, 9 * Years);
 CHFCurveBuilder2.AddSwap(0.0016365, 10 * Years);
 CHFCurveBuilder2.AddSwap(0.003, 12 * Years);
 CHFCurveBuilder2.AddSwap(0.004525, 15 * Years);
 CHFCurveBuilder2.AddSwap(0.0063, 20 * Years);
 CHFCurveBuilder2.AddSwap(0.00735, 25 * Years);
 CHFCurveBuilder2.AddSwap(0.007825, 30 * Years);
 //
 // get relinkable curve handle from CHF curve builder and re-print discount factors
 curveHandle = CHFCurveBuilder2.GetCurveHandle();
 PrintCHFDiscountFactors(settlementDate, curveHandle);
 return 0;
}
//
template <typename T> Rate GetSwapRate(Period swapMaturity, Date settlementDate, Period floatingLegTenor, 
 Handle<YieldTermStructure>& curveHandle, Calendar calendar, DayCounter fixedLegDayCount, Period fixedLegTenor)
{
 // using quantlib factory method for building vanilla swap
 boost::shared_ptr<IborIndex> index(new T(floatingLegTenor, curveHandle));
 VanillaSwap swap = MakeVanillaSwap(swapMaturity, index)
  .withEffectiveDate(settlementDate)
  .withFixedLegCalendar(calendar)
  .withFixedLegConvention(index->businessDayConvention())
  .withFixedLegDayCount(fixedLegDayCount)
  .withFixedLegTenor(fixedLegTenor)
  .withFloatingLegCalendar(calendar);
 return swap.fairRate();
}
//
void PrintUSDRates(Date settlementDate, Handle<YieldTermStructure>& curveHandle)
{
 Calendar calendar = UnitedKingdom();
 std::cout << std::endl;
 // print USD deposit rates
 std::cout << curveHandle->zeroRate(calendar.adjust(settlementDate + 1 * Weeks), 
  Actual360(), Simple).rate() << std::endl;
 std::cout << curveHandle->zeroRate(calendar.adjust(settlementDate + 1 * Months), 
  Actual360(), Simple).rate() << std::endl;
 std::cout << curveHandle->zeroRate(calendar.adjust(settlementDate + 2 * Months), 
  Actual360(), Simple).rate() << std::endl;
 std::cout << curveHandle->zeroRate(calendar.adjust(settlementDate + 3 * Months), 
  Actual360(), Simple).rate() << std::endl;
 // print USD forward rates
 std::cout << curveHandle->forwardRate(calendar.adjust(settlementDate + 3 * Months), 
  calendar.adjust(settlementDate + 6 * Months), Actual360(), Simple).rate() << std::endl;
 std::cout << curveHandle->forwardRate(calendar.adjust(settlementDate + 6 * Months), 
  calendar.adjust(settlementDate + 9 * Months), Actual360(), Simple).rate() << std::endl;
 std::cout << curveHandle->forwardRate(calendar.adjust(settlementDate + 9 * Months), 
  calendar.adjust(settlementDate + 12 * Months), Actual360(), Simple).rate() << std::endl;
 // print USD swap rates
 std::cout << GetSwapRate<USDLibor>(2 * Years, settlementDate, 3 * Months, 
  curveHandle, calendar, Thirty360(), 1 * Years) << std::endl;
 std::cout << GetSwapRate<USDLibor>(3 * Years, settlementDate, 3 * Months, 
  curveHandle, calendar, Thirty360(), 1 * Years) << std::endl;
 std::cout << GetSwapRate<USDLibor>(4 * Years, settlementDate, 3 * Months, 
  curveHandle, calendar, Thirty360(), 1 * Years) << std::endl;
 std::cout << GetSwapRate<USDLibor>(5 * Years, settlementDate, 3 * Months, 
  curveHandle, calendar, Thirty360(), 1 * Years) << std::endl;
 std::cout << GetSwapRate<USDLibor>(6 * Years, settlementDate, 3 * Months, 
  curveHandle, calendar, Thirty360(), 1 * Years) << std::endl;
 std::cout << GetSwapRate<USDLibor>(7 * Years, settlementDate, 3 * Months, 
  curveHandle, calendar, Thirty360(), 1 * Years) << std::endl;
 std::cout << GetSwapRate<USDLibor>(8 * Years, settlementDate, 3 * Months, 
  curveHandle, calendar, Thirty360(), 1 * Years) << std::endl;
 std::cout << GetSwapRate<USDLibor>(9 * Years, settlementDate, 3 * Months, 
  curveHandle, calendar, Thirty360(), 1 * Years) << std::endl;
 std::cout << GetSwapRate<USDLibor>(10 * Years, settlementDate, 3 * Months, 
  curveHandle, calendar, Thirty360(), 1 * Years) << std::endl;
 std::cout << GetSwapRate<USDLibor>(11 * Years, settlementDate, 3 * Months, 
  curveHandle, calendar, Thirty360(), 1 * Years) << std::endl;
 std::cout << GetSwapRate<USDLibor>(12 * Years, settlementDate, 3 * Months, 
  curveHandle, calendar, Thirty360(), 1 * Years) << std::endl;
 std::cout << GetSwapRate<USDLibor>(15 * Years, settlementDate, 3 * Months, 
  curveHandle, calendar, Thirty360(), 1 * Years) << std::endl;
 std::cout << GetSwapRate<USDLibor>(20 * Years, settlementDate, 3 * Months, 
  curveHandle, calendar, Thirty360(), 1 * Years) << std::endl;
 std::cout << GetSwapRate<USDLibor>(25 * Years, settlementDate, 3 * Months, 
  curveHandle, calendar, Thirty360(), 1 * Years) << std::endl;
 std::cout << GetSwapRate<USDLibor>(30 * Years, settlementDate, 3 * Months, 
  curveHandle, calendar, Thirty360(), 1 * Years) << std::endl;
 std::cout << GetSwapRate<USDLibor>(40 * Years, settlementDate, 3 * Months, 
  curveHandle, calendar, Thirty360(), 1 * Years) << std::endl;
 std::cout << GetSwapRate<USDLibor>(50 * Years, settlementDate, 3 * Months, 
  curveHandle, calendar, Thirty360(), 1 * Years) << std::endl;
}
//
void PrintCHFDiscountFactors(Date settlementDate, Handle<YieldTermStructure>& curveHandle)
{
 // print CHF discount factors for every 6 months up to 30 years
 std::cout << std::endl;
 for(unsigned int i = 1; i != 61; i++)
 {
  std::cout << curveHandle->discount(UnitedKingdom()
   .adjust(settlementDate + (6 * i) * Months)) << std::endl;
 }
}

Implementing statistics gatherer design by using Boost library in C++

As a result of Monte Carlo simulation process, we get a lot of simulated values. After this, we usually want to calculate a set of desired statistics from those simulated values. A great statistics gatherer design has been thoroughly presented in The book by Mark Joshi. However, just for the curious and for the sake of trying to create something original (while respecting traditions), I wanted to implement my own statistics gatherer design, in which several different statistical measures could be calculated independently by using the same set of processed data. Moreover, the number of such calculation objects and their implementations for algorithmic details would be fully configurable. Finally, I wanted to avoid any direct coupling between data (object, which holds the source data) and algorithms (objects, which are performing the actual statistical calculations).

Why bother?

In the following example program below, a sample of discounted payoffs from Monte Carlo process has been hard-coded and two statistical measures are then calculated and printed out.

void MonteCarloProcess()
{
 // process has been producing the following discounted payoffs
 // for this sample : average = 10.4881, standard error = 1.58502
 double discountedPayoffs[] = { 18.5705, 3.31508, 0.0, 3.64361, 0.0, 0.0, 
  47.2563, 10.6534, 85.5559, 0.0, 30.2363, 0.0, 17.8391, 2.15396, 0.0, 
  66.587, 0.0, 9.19303, 0.0, 0.0, 24.2946, 29.6556, 0.0, 0.0, 0.0, 65.926, 
  0.0, 14.0329, 1.43328, 0.0, 0.0, 0.0, 1.37088, 0.0, 2.49095, 21.4755, 
  36.5432, 0.0, 16.8795, 0.0, 0.0, 0.0, 19.8927, 11.3132, 37.3946, 10.2666, 
  26.1932, 0.0, 0.551356, 29.7159, 0.0, 31.5357, 0.0, 0.0, 4.64357, 4.45376, 
  21.6076, 12.693, 16.0065, 0.0, 0.0, 0.0, 0.0, 25.9665, 18.7169, 2.55222, 
  25.6431, 8.5027, 0.0, 0.0, 29.8704, 0.0, 22.7266, 22.8463, 0.0, 0.0, 0.0, 
  0.0, 4.90832, 13.2787, 0.0, 0.0, 9.77076, 24.5855, 12.6094, 0.0, 0.0, 1.92343, 
  5.66301, 0.0, 0.0, 13.6968, 0.0, 0.0, 35.2159, 0.0, 8.44648, 7.21964, 0.0, 19.2949 };
 //
 // dump array data into vector
 std::vector<double> data(discountedPayoffs, std::end(discountedPayoffs));
 //
 // calculate and print mean and standard error
 double mean = AlgorithmLibrary::Mean(data);
 double standardError = AlgorithmLibrary::StandardError(data);
 std::cout << mean << std::endl;
 std::cout << standardError << std::endl;
}

There is nothing fundamentally wrong in this example, but a great deal of it could be done in a bit different manner. The scheme presented in this example offers a chance to explore some modern tools for implementing flexible and configurable C++ programs.

Chinese walls

Personally, I would prefer to separate data (discounted payoffs) and algorithms (mean, standard error) completely. Ultimately, I would like to have a design, in which a process (monte carlo) is generating results and adding those results into a separate container. Whenever this process is finished, it will send a reference of this container for several entities, which will calculate all required statistics independently.

There are many ways to implement this kind of a scheme, but I have been heavily influenced by delegates stuff I have learned from C#. Needless to say, the equivalent mechanism for C# delegate in C++ is a function pointer. However, instead of raw function pointer I will use Boost.Function and Boost.Bind libraries.

My new design proposal will have the following components :

• AlgorithmLibrary for calculating different statistical measures. The header file for this contains collection of methods for calculating different statistical measures.
• ResultContainer class for storing processed results and function pointers (boost function) which are sending processed results for further calculations.
• StatisticsElement class for storing a value of statistical measure and function pointer (boost function) for an algorithm, which can be used for calculating required statistical measure.

 

Strangers in the night

In the first stage, StatisticsElement object will be created. This object will host a single statistical measure and a pointer to an algorithm (boost function) for calculating this specific measure. By giving required algorithm as a pointer (boost function), the object is not tied with any hard-coded algorithm. In the case there would be a need for another type for algorithm implementation for calculating required statistical measure, the scheme is flexible enough to be adjusted. In the constructor of this class, we are giving this pointer to an algorithm (boost function). Moreover, this object will be indirectly connected with ResultContainer object with a function pointer (boost function). Function pointer (boost function) will be created and binded (boost bind) with the specific method (process) of StatisticsElement object.

In the second stage, ResultContainer object will be created and all previously created function pointers (boost function, boost bind) for processing calculation results will be added into ResultContainer. This object is ultimately being shared with a process. Process will generate its results (double) and these results will be added into container object. When a process is finished, container object method (sendResults) will be called. The call for this method will trigger a loop, which will iterate through a vector of function pointers (boost function) and sending a reference for a result vector to all connected StatisticsElement objects.

Finally, the client program (in this example : main) will request the calculated statistical measures directly from all StatisticsElement objects. It should be stressed, that these two objects described above, do not have any knowledge about each other at any point. ResultContainer is just storing updates from a process and finally sending results "to somewhere" when the processing is over. StatisticsElement objects are processing their own calculation procedures as soon as they will receive results "from somewhere". It should also be noted, that this design actually implements observer pattern, where ResultContainer object is Observable and StatisticalElement objects are Observers.

Proposal

// ResultContainer.h
#pragma once
#include <vector>
#include <boost\function.hpp>
//
// class for storing processed results and function pointers 
// which are sending results for further processing
class ResultContainer
{
private:
 // container for storing processed results
 std::vector<double> results;
 // container for storing function pointers
 std::vector<boost::function<void
  (const std::vector<double>&)>> resultSenders;
public:
 // method for adding one processed result into container
 void addResult(double result);
 // method for adding one function pointer into container
 void addResultSender(boost::function<void
  (const std::vector<double>&)> resultSender);
 // method for sending all processed results
 void sendResults() const;
};
//
//
//
// StatisticsElement.h
#pragma once
#include <vector>
#include <boost\function.hpp>
//
// class for storing a value of statistical measure and 
// function pointer for an algorithm which can be used 
// for calculating required statistical measure
class StatisticsElement
{
private:
 // statistical measure
 double statisticsValue;
 // function pointer to an algorithm which calculates 
 // required statistical measure
 boost::function<double(const std::vector<double>&)> algorithm;
public:
 // parameter constructor
 StatisticsElement(boost::function<double
  (const std::vector<double>&)> algorithm);
 // method for processing data in order to 
 // calculate required statistical measure
 void process(const std::vector<double>& data);
 // method (overloaded operator) for accessing 
 // calculated statistical measure
 double operator()() const;
};
//
//
//
// AlgorithmLibrary.h
#pragma once
#include <vector>
#include <numeric>
//
// algorithms library for calculating statistical measures
namespace AlgorithmLibrary
{
 // calculate arithmetic average
 double Mean(const std::vector<double>& data) 
 { 
  return std::accumulate(data.begin(), data.end(), 0.0) 
   / data.size(); 
 }
 // calculate standard error estimate
 double StandardError(const std::vector<double>& data) 
 { 
  double mean = AlgorithmLibrary::Mean(data);
  double squaredSum = std::inner_product(data.begin(), 
   data.end(), data.begin(), 0.0);
  return std::sqrt(squaredSum / data.size() - mean * mean) 
   / std::sqrt(data.size());
 }
}
//
//
//
// ResultContainer.cpp
#include "ResultContainer.h"
//
void ResultContainer::addResult(double result)
{
 results.push_back(result);
}
void ResultContainer::addResultSender(boost::function<void
 (const std::vector<double>&)> resultSender)
{
 resultSenders.push_back(resultSender);
}
void ResultContainer::sendResults() const
{
 std::vector<boost::function<void
  (const std::vector<double>&)>>::const_iterator it;
 for(it = resultSenders.begin(); it != resultSenders.end(); it++)
 {
  (*it)(results);
 }
}
//
//
//
// StatisticsElement.cpp
#include "StatisticsElement.h"
//
StatisticsElement::StatisticsElement(boost::function<double
 (const std::vector<double>&)> algorithm) 
 : statisticsValue(0.0), algorithm(algorithm) 
{
 //
}
void StatisticsElement::process(const std::vector<double>& data)
{
 if(algorithm != NULL) statisticsValue = algorithm(data);
}
double StatisticsElement::operator()() const
{
 return statisticsValue;
}
//
//
//
// MainProgram.cpp
#include <boost\bind.hpp>
#include <iostream>
#include "AlgorithmLibrary.h"
#include "ResultContainer.h"
#include "StatisticsElement.h"
//
void MonteCarloProcess(ResultContainer& resultContainer)
{
 // process has been producing the following discounted payoffs
 // for this sample : average = 10.4881, standard error = 1.58502
 double discountedPayoffs[] = { 18.5705, 3.31508, 0.0, 3.64361, 0.0, 0.0, 
  47.2563, 10.6534, 85.5559, 0.0, 30.2363, 0.0, 17.8391, 2.15396, 0.0, 
  66.587, 0.0, 9.19303, 0.0, 0.0, 24.2946, 29.6556, 0.0, 0.0, 0.0, 65.926, 
  0.0, 14.0329, 1.43328, 0.0, 0.0, 0.0, 1.37088, 0.0, 2.49095, 21.4755, 
  36.5432, 0.0, 16.8795, 0.0, 0.0, 0.0, 19.8927, 11.3132, 37.3946, 10.2666, 
  26.1932, 0.0, 0.551356, 29.7159, 0.0, 31.5357, 0.0, 0.0, 4.64357, 4.45376, 
  21.6076, 12.693, 16.0065, 0.0, 0.0, 0.0, 0.0, 25.9665, 18.7169, 2.55222, 
  25.6431, 8.5027, 0.0, 0.0, 29.8704, 0.0, 22.7266, 22.8463, 0.0, 0.0, 0.0, 
  0.0, 4.90832, 13.2787, 0.0, 0.0, 9.77076, 24.5855, 12.6094, 0.0, 0.0, 1.92343, 
  5.66301, 0.0, 0.0, 13.6968, 0.0, 0.0, 35.2159, 0.0, 8.44648, 7.21964, 0.0, 19.2949 };
 //
 // dump array data into vector
 std::vector<double> data(discountedPayoffs, std::end(discountedPayoffs));
 // create vector iterator, loop through data and add items into result container object
 std::vector<double>::const_iterator it;
 for(it = data.begin(); it != data.end(); it++)
 {
  resultContainer.addResult(*it);
 }
 // trigger result processing for all 'connected' statistical element objects
 resultContainer.sendResults();
}
//
int main()
{
 // create : function pointer to mean algorithm, statistics element 
 // and function pointer to process method of this statistics element
 boost::function<double(const std::vector<double>&)> meanAlgorithm = 
  AlgorithmLibrary::Mean;
 StatisticsElement mean(meanAlgorithm);
 boost::function<void(const std::vector<double>&)> resultSenderForMean = 
  boost::bind(&StatisticsElement::process, &mean, _1);
 //
 // create : function pointer to standard error algorithm, statistics element and 
 // function pointer to process method of this statistics element
 boost::function<double(const std::vector<double>&)> standardErrorAlgorithm = 
  AlgorithmLibrary::StandardError;
 StatisticsElement standardError(standardErrorAlgorithm);
 boost::function<void(const std::vector<double>&)> resultSenderForStandardError = 
  boost::bind(&StatisticsElement::process, &standardError, _1);
 //
 // create : result container and add previously created function 
 // pointers (senders) into container
 ResultContainer resultContainer;
 resultContainer.addResultSender(resultSenderForMean);
 resultContainer.addResultSender(resultSenderForStandardError);
 //
 // run (hard-coded) monte carlo process
 MonteCarloProcess(resultContainer);
 //
 // print results from the both statistics elements
 std::cout << mean() << std::endl;
 std::cout << standardError() << std::endl;
 return 0;
}

Help

Concerning the actual installation and configuration of Boost libraries with compiler, there is a great tutorial by eefelix available in youtube. For using Boost libraries, there is a document available, written by Dimitri Reiswich. Personally, I would like to present my appreciations for these persons for their great contribution.

Thanks for reading my blog. Have a pleasant wait for the Christmas.
-Mike

Simulating Discount Factor and Forward Rate Curves using Monte Carlo in C#

The process of creating discount factor and forward rate curves with traditional bootstrapping algorithm was presented in the last post. In this post we are going to do the same thing, but following a bit different approach.

There are two ways to use the market data when creating these curves. The first approach is to fit the curve exactly into the current market data (calibration, bootstrapping). Another approach is to select an interest rate model, estimate all required parameters from the current market data and finally build the curves by using Monte Carlo method, for example. The advantage of the first approach is, that we are able to reprice (or replicate) the market with the resulting curves, while this (almost surely) will never happen with the second approach.

Monte Carlo method itself is widely used for a wide array of complex calculation tasks. Just for an example, calculating CVA for an interest rate swap requires us to calculate swap PV for a number of chosen timepoints happening in the future, before the expiration of a swap contract. For this complex calculation task, we can select an interest rate model and use Monte Carlo to calculate swap PV for any point in time in the future.

PROJECT OUTCOME

In this project, we are following this second approach and use Vasicek one-factor interest rate model in our example. The resulting C# program uses Monte Carlo for simulating discount curve (a set of zero-coupon bonds). After this, client can request discount factor DF(0, T) for any given maturity or forward rate FWD(t, T) for any given period in the future. The resulting example program is a console project. However, with all the information available in this blog, one should be able to interface the program with Excel, if so desired. Needless to say, the part of the program which actually creates the curves can be used as a part of any other C# project.

PROGRAM DESIGN

For the purpose of getting some conceptual overview, about how the objects in the program are connected with each other, a rough UML class diagram is shown in the picture below.

YELLOW
First of all, Client is requesting MonteCarloEngine object to simulate the short rate paths. MonteCarloEngine is sending all simulated paths for MonteCarloCurveProcessor, using delegate method. For the simulation task, MonteCarloEngine needs (aggregation) three different objects : RandomNumber, SDE and Discretization. These objects are going to be created by HardCodedExampleBuilder (aggregation), which is a specialization for abstract MonteCarloBuilder class.

BLUE
All the needed information concerning the short rate model is embedded in EulerDiscretization (specialization for abstract Discretization) and Vasicek (specialization for abstract SDE) objects. These objects are aggregated in MonteCarloEngine.

GREEN
For performing Monte Carlo simulation task, random numbers are needed. NAG_BasicRandomNumber object (specialization for abstract RandomNumber) is aggregated in MonteCarloEngine. This object is using (aggregation) two other objects to perform the task. NAG_BasicGenerator (specialization for abstract UniformRandomGenerator) and StandardNormal (specialization for abstract InverseTransformation).

With the design used in the program, a client is able to change

• uniform random number generator
• the class, which is processing generated uniform randon numbers
• inverse transformation method for mapping generated uniform random numbers into a given probability distribution
• stochastic differential equation
• discretization scheme for a given stochastic differential equation

So, we are having a lot of flexibility with this presented design.

It should be noted, that the basic task of this program is to simulate paths for any given (one-factor) stochastic differential equation, using a given discretization scheme. MonteCarloEngine is kind of a Mediator object, which is performing this task. MonteCarloCurveProcessor object, which has been added into this design later, is only using the results generated by MonteCarloEngine for creating a set of discount factors and calculating forward rates, when requested by the client.

THE PROGRAM

Create a new C# console project and copyPaste the following program into a new cs-file. Remember to add reference to NAG .NET library. If you are not registrated NAG user, create a new implementation for UniformRandomGenerator and completely remove the current uniform generator implementation.

using System;
using System.Collections.Generic;
using System.Linq;
using NagLibrary;

namespace RandomDesign
{
    // public delegates
    public delegate double Interpolator(Dictionary<double, double> data, double key);
    public delegate void PathSender(double[] path);
    public delegate void ProcessStopper();
    public delegate int Seeder();
    //
    // Client
    class Program
    {
        static void Main(string[] args)
        {
            try
            {
                // create required objects, data and run monte carlo engine
                double tenor = 0.25;
                double maturity = 10.0;
                int paths = 250000;
                int steps = 1000;
                MonteCarloEngine engine = new MonteCarloEngine(new HardCodedExampleBuilder(), paths, steps);
                MonteCarloCurveProcessor curve = new MonteCarloCurveProcessor(maturity, tenor, Interpolation.LinearInterpolation);
                engine.sendPath += curve.Process;
                engine.stopProcess += curve.Calculate;
                engine.run();
                //
                // after engine run, discount factors and simply-compounded 
                // forward rates can be requested by the client
                Console.WriteLine("printing quarterly discount factors >");
                int n = Convert.ToInt32(maturity / tenor);
                for (int i = 1; i <= n; i++)
                {
                    double T = (tenor * i);
                    double df = curve.GetDF(T);
                    Console.WriteLine(Math.Round(df, 4));
                }
                Console.WriteLine();
                //
                Console.WriteLine("printing quarterly forward rates >");
                for (int i = 2; i <= n; i++)
                {
                    double t = (tenor * (i - 1));
                    double T = (tenor * i);
                    double f = curve.GetFWD(t, T);
                    Console.WriteLine(Math.Round(f, 4));
                }
                Console.WriteLine();
            }
            catch (Exception e)
            {
                Console.WriteLine(e.Message);
            }
        }
    }
    //
    //
    //
    // class for :
    // - processing stochastic paths sent by monte carlo engine
    // - calculating zero-coupon bond prices using collected information on stochastic paths
    // - retrieving discount factors or simply-compounded forward rates for any given period
    public class MonteCarloCurveProcessor
    {
        private Dictionary<double, double> discountCurve;
        private Interpolator interpolator;
        private double maturity;
        private double tenor;
        private int nBonds;
        private int nSteps;
        private double[] paths;
        private bool hasDimension = false;
        private int counter;
        //
        public MonteCarloCurveProcessor(double maturity, double tenor, Interpolator interpolator)
        {
            discountCurve = new Dictionary<double, double>();
            this.interpolator = interpolator;
            this.maturity = maturity;
            this.tenor = tenor;
            nBonds = Convert.ToInt32(this.maturity / this.tenor);
            counter = 0;
        }
        public void Process(double[] path)
        {
            if (!hasDimension)
            {
                nSteps = path.GetLength(0);
                paths = new double[nSteps];
                hasDimension = true;
            }
            // add path to all existing paths
            for (int i = 0; i < nSteps; i++)
            {
                paths[i] += path[i];
            }
            counter++;
        }
        public void Calculate()
        {
            // calculate the average path
            for (int i = 0; i < nSteps; i++)
            {
                paths[i] /= counter;
            }
            // integrate zero-coupon bond prices
            double dt = maturity / (nSteps - 1);
            int takeCount = Convert.ToInt32(tenor / dt);
            for (int i = 0; i < nBonds; i++)
            {
                double integral = paths.ToList<double>().Take(takeCount * (i + 1) + 1).Sum();
                discountCurve.Add(tenor * (i + 1), Math.Exp(-integral * dt));
            }
        }
        public double GetDF(double T)
        {
            // interpolate discount factor from discount curve
            return interpolator(discountCurve, T);
        }
        public double GetFWD(double t, double T)
        {
            // using discount factors, calculate forward discount 
            // factor and convert it into simply-compounded forward rate
            double df_T = interpolator(discountCurve, T);
            double df_t = interpolator(discountCurve, t);
            double fdf = (df_T / df_t);
            return ((1 / fdf) - 1) / (T - t);
        }
    }
    //
    //
    //
    // class hierarchy for simulation object builders
    public abstract class MonteCarloBuilder
    {
        public abstract Tuple<SDE, Discretization, RandomNumber> build();
    }
    // implementation for hard-coded example
    public class HardCodedExampleBuilder : MonteCarloBuilder
    {
        public override Tuple<SDE, Discretization, RandomNumber> build()
        {
            // create objects for generating standard normally distributed random numbers
            UniformRandomGenerator uniformRandom = new NAG_BasicGenerator(Seed.GetGUIDSeed);
            InverseTransformation transformer = new StandardNormal();
            RandomNumber nag = new NAG_BasicRandomNumber(uniformRandom, transformer); 
            //
            // create stochastic differential equation and discretization objects
            SDE vasicek = new Vasicek(0.1, 0.29, 0.0185);
            Discretization euler = new EulerDiscretization(vasicek, 0.02, 10.0);
            //
            return new Tuple<SDE, Discretization, RandomNumber>(vasicek, euler, nag);
        }
    }
    //
    //
    //
    // class hierarchy for stochastic differential equations
    public abstract class SDE
    {
        public abstract double drift(double s, double t);
        public abstract double diffusion(double s, double t);
    }
    // implementation for Vasicek one-factor interest rate model
    public class Vasicek : SDE
    {
        private double longTermRate;
        private double reversionSpeed;
        private double rateVolatility;
        //
        public Vasicek(double longTermRate, double reversionSpeed, double rateVolatility)
        {
            this.longTermRate = longTermRate;
            this.reversionSpeed = reversionSpeed;
            this.rateVolatility = rateVolatility;
        }
        public override double drift(double s, double t)
        {
            return reversionSpeed * (longTermRate - s);
        }
        public override double diffusion(double s, double t)
        {
            return rateVolatility;
        }
    }
    //
    //
    //
    // class hierarchy for discretization schemes
    public abstract class Discretization
    {
        protected SDE sde;
        protected double initialPrice;
        protected double expiration;
        //
        // read-only properties for initial price and expiration
        public double InitialPrice { get { return initialPrice; } }
        public double Expiration { get { return expiration; } }
        public Discretization(SDE sde, double initialPrice, double expiration)
        {
            this.sde = sde;
            this.initialPrice = initialPrice;
            this.expiration = expiration;
        }
        public abstract double next(double s, double t, double dt, double rnd);
    }
    // implementation for Euler discretization scheme
    public class EulerDiscretization : Discretization
    {
        public EulerDiscretization(SDE sde, double initialPrice, double expiration)
            : base(sde, initialPrice, expiration) 
        { 
            //
        }
        public override double next(double s, double t, double dt, double rnd)
        {
            return s + sde.drift(s, t) * dt + sde.diffusion(s, t) * Math.Sqrt(dt) * rnd;
        }
    }
    //
    //
    //
    // mediator class for creating stochastic paths
    public class MonteCarloEngine
    {
        private SDE sde;
        private Discretization discretization;
        private RandomNumber random;
        private long paths;
        private int steps;
        public event PathSender sendPath;
        public event ProcessStopper stopProcess;
        //
        public MonteCarloEngine(MonteCarloBuilder builder, long paths, int steps)
        {
            Tuple<SDE, Discretization, RandomNumber> components = builder.build();
            this.sde = components.Item1;
            this.discretization = components.Item2;
            this.random = components.Item3;
            this.paths = paths;
            this.steps = steps;
        }
        public void run()
        {
            // create skeleton for path values
            double[] path = new double[steps + 1]; 
            double dt = discretization.Expiration / steps;
            double vOld = 0.0; double vNew = 0.0;
            //
            for (int i = 0; i < paths; i++)
            {
                // request random numbers for a path
                double[] rand = random.GetPath(steps);
                path[0] = vOld = discretization.InitialPrice;
                //
                for (int j = 1; j <= steps; j++)
                {
                    // get next path value using a given discretization scheme
                    vNew = discretization.next(vOld, (dt * j), dt, rand[j - 1]);
                    path[j] = vNew; vOld = vNew;
                }
                sendPath(path); // send one path to client (MonteCarloCurveProcessor) for processing
            }
            stopProcess(); // notify client (MonteCarloCurveProcessor)
        }
    }
    //
    //
    //
    // class hierarchy for uniform random generators
    public abstract class UniformRandomGenerator
    {
        public abstract double[] Get(int n);
    }
    // implementation for NAG basic uniform pseudorandom number generator
    public class NAG_BasicGenerator : UniformRandomGenerator
    {
        private Seeder seeder;
        private int[] seed;
        private const int baseGenerator = 1; // hard-coded NAG basic generator
        private const int subGenerator = 1; // hard-coded sub-generator (not referenced)
        //
        public NAG_BasicGenerator(Seeder seeder)
        {
            this.seeder = seeder;
        }
        public override double[] Get(int n)
        {
            seed = new int[] { seeder() };
            double[] rnd = new double[n];
            int errorState = 0;
            G05.G05State g05State = new G05.G05State(baseGenerator, subGenerator, seed, out errorState);
            if (errorState != 0) throw new Exception("NAG : g05State error");
            G05.g05sq(n, 0.0, 1.0, g05State, rnd, out errorState);
            if (errorState != 0) throw new Exception("NAG : g05sq error");
            return rnd;
        }
    }
    //
    //
    //
    // class hierarchy for random number processors 
    public abstract class RandomNumber
    {
        public readonly UniformRandomGenerator generator;
        public readonly InverseTransformation transformer;
        //
        public RandomNumber(UniformRandomGenerator generator, 
            InverseTransformation transformer = null)
        {
            this.generator = generator;
            this.transformer = transformer;
        }
        public abstract double[] GetPath(int nSteps);
    }
    // implementation for class processing NAG basic random numbers
    public class NAG_BasicRandomNumber : RandomNumber
    {
        public NAG_BasicRandomNumber(UniformRandomGenerator generator,
            InverseTransformation transformer = null)
            : base(generator, transformer)
        {
            // initialize base class data members
        }
        public override double[] GetPath(int nSteps)
        {
            // create and return one random number path
            double[] path = base.generator.Get(nSteps);
            //
            // conditionally, map uniform random numbers to a given distribution
            if (base.transformer != null)
            {
                for (int i = 0; i < nSteps; i++)
                {
                    path[i] = base.transformer.Inverse(path[i]);
                }
            }
            return path;
        }
    }
    //
    //
    //
    // class hierarchy for all inverse transformation classes
    public abstract class InverseTransformation
    {
        public abstract double Inverse(double x);
    }
    // mapping uniform random variate to standard normal distribution
    public class StandardNormal : InverseTransformation
    {
        public override double Inverse(double x)
        {
            const double a1 = -39.6968302866538; const double a2 = 220.946098424521;
            const double a3 = -275.928510446969; const double a4 = 138.357751867269;
            const double a5 = -30.6647980661472; const double a6 = 2.50662827745924;
            //
            const double b1 = -54.4760987982241; const double b2 = 161.585836858041;
            const double b3 = -155.698979859887; const double b4 = 66.8013118877197;
            const double b5 = -13.2806815528857;
            //
            const double c1 = -7.78489400243029E-03; const double c2 = -0.322396458041136;
            const double c3 = -2.40075827716184; const double c4 = -2.54973253934373;
            const double c5 = 4.37466414146497; const double c6 = 2.93816398269878;
            //
            const double d1 = 7.78469570904146E-03; const double d2 = 0.32246712907004;
            const double d3 = 2.445134137143; const double d4 = 3.75440866190742;
            //
            const double p_low = 0.02425; const double p_high = 1.0 - p_low;
            double q = 0.0; double r = 0.0; double c = 0.0;
            //
            if ((x > 0.0) && (x < 1.0))
            {
                if (x < p_low)
                {
                    // rational approximation for lower region
                    q = Math.Sqrt(-2.0 * Math.Log(x));
                    c = (((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6)
                        / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1);
                }
                if ((x >= p_low) && (x <= p_high))
                {
                    // rational approximation for middle region
                    q = x - 0.5; r = q * q;
                    c = (((((a1 * r + a2) * r + a3) * r + a4) * r + a5) * r + a6) * q
                        / (((((b1 * r + b2) * r + b3) * r + b4) * r + b5) * r + 1);
                }
                if (x > p_high)
                {
                    // rational approximation for upper region
                    q = Math.Sqrt(-2 * Math.Log(1 - x));
                    c = -(((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6)
                        / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1);
                }
            }
            else
            {
                // throw if given x value is out of bounds
                // (less or equal to 0.0, greater or equal to 1.0)
                throw new Exception("StandardNormal : out of bounds error");
            }
            return c;
        }
    }
    //
    //
    //
    // static library class consisting of seeding methods for uniform random number generators
    public static class Seed
    {
        public static int GetGUIDSeed() { return Math.Abs(Guid.NewGuid().GetHashCode()); }
    }
    //
    //
    //
    // static library class for interpolation methods
    public static class Interpolation
    {
        public static double LinearInterpolation(Dictionary<double, double> data, double key)
        {
            double value = 0.0;
            int n = data.Count;
            //
            // boundary checkings
            if ((key < data.ElementAt(0).Key) || (key > data.ElementAt(data.Count - 1).Key))
            {
                if (key < data.ElementAt(0).Key) value = data.ElementAt(0).Value;
                if (key > data.ElementAt(data.Count - 1).Key) value = data.ElementAt(data.Count - 1).Value;
            }
            else
            {
                // iteration through all existing elements
                for (int i = 0; i < n; i++)
                {
                    if ((key >= data.ElementAt(i).Key) && (key <= data.ElementAt(i + 1).Key))
                    {
                        double t = data.ElementAt(i + 1).Key - data.ElementAt(i).Key;
                        double w = (key - data.ElementAt(i).Key) / t;
                        value = data.ElementAt(i).Value * (1 - w) + data.ElementAt(i + 1).Value * w;
                        break;
                    }
                }
            }
            return value;
        }
    }
}

VASICEK PARAMETERS ESTIMATION

Parameters (long-term mean rate, reversion speed and rate volatility) for Vasicek model have to be estimated from the market data. For this task, we can use OLS method or Maximum Likelihood method. In this case, the both methods should end up with the same set of result parameters. Concrete hands-on example on parameter estimation with the both methods for Ornstein-Uhlenbeck model (Vasicek), has been presented in this great website by Thijs van den Berg.

When using OLS method, the task can easily be performed in Excel using Data Analysis tools. In the screenshot presented below, I have replicated the parameter estimation using Excel regression tool and example dataset from the website given above.

The parameter estimation procedure with OLS method is straightforward in Excel and should be relatively easy to perform on any new chosen set of market data. Needless to say, the real deal in the parameter estimation approach is linked with the market data sample properties. The central question is the correct period to include into a sample for estimating such parameters.

Thanks a lot for using your precious time and reading my blog.

-Mike Juniperhill