Sunday, December 13, 2015

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

Sunday, August 23, 2015

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

Saturday, August 15, 2015

Bootstrapping Libor Discount Factor and Forward Rate Curves using C# and Excel-DNA

Defining a set of correct discount factors and forward rates is the cornerstone of valuing any Libor-based products. As we know, calculation of present value for a Libor cash flow seems to be almost too easy : first define a cash flow happening in the future by using forward rates, then use discount factors to get present value for that cash flow. Then we also know, that the real deal is always getting those forward rates and discount factors from somewhere, in the first place.

For the task described above, we usually have our front office system or tools provided by some third-party vendor, which are then performing all those required complex calculations quietly behind the scenes. Figures are popping up in the screen, finding their way into some strange reports and everyone are happy. Needless to say, we should be able to perform those calculations also by hand, if so required. This is only my personal opinion on the matter, of course. During this small project, we are creating a program, which takes in market data as input parameter. From this market data, the program is then performing all those complex calculations (bootstrapping, interpolations) and returning discount factors and simply-compounded Libor forward rates.

CURVOLOGY

When constructing discount factors and forward rates from market data, one has to make a lot of different decisions. Just for an example
  • How to construct short-end of the curve? For what maturities are we going to use cash securities?
  • How to construct mid-part of the curve? Are we using futures contracts on Libor or Forward Rate Agreements? 
  • Are we going to make adjustments for convexity effect?
  • How to construct long-end of the curve? For what maturity we start to use swap rates?
  • What kind of interpolation methods are we going to use?
  • Do we interpolate discount factors or spot rates?
Due to this high degrees of freedom embedded in this task, it should not be a miracle, that two persons (given the same set of market data) most probably will end up having completely different set of discount factors and forward rates. As I have been studying this topic and getting a bit more familiar with the methods and procedures used by some well-known authors, I have personally come to the conclusion that there are some common guidelines and general conventions, but absolutely not any one universally correct way to do this task.

MARKET DATA

For this project, I have been using dataset and curve creation procedures as presented in Richard Flavell's excellent book Swaps and Other Derivatives. In this book, Flavell is giving complete treatment, how to construct Libor zero curve (discount factors and forward rates). One of the things why I like so much this book is the fact, that Flavell is not assuming anything. Instead, he is always deriving everything from the beginning. Great thing is also, that the book has all the example calculations dissected in the attached Excel workbooks. Personally, these Excels have been sort of a goldmine for me, when I have been studying different valuation issues.

PROJECT OUTCOME

The result of this project is XLL addin for Excel, which calculates discount factors and simply-compounded forward rates for USD Libor. We will use Excel-DNA for interfacing our C# program with Excel. The scheme what has been used for this project, has been fully explained in this blog post. So, carefully follow the instructions given in that post, when interfacing the program with Excel. Needless to say, one should be able to use the core classes of the program (the classes which are actually creating the curves) in any other C# project, without Excel interfacing.

PROGRAM DESIGN

Rough UML class diagram for this program is shown below. The purpose of this diagram is to get some conceptual overview, how the program is working and how the objects are related with each other. Since I am not UML expert, I am presenting my deepest apologies for its incompleteness and possible errors.

In the first stage, Client (CreateUSDLiborCurves, which will be presented in the section Excel Interfacing) is requesting ExcelBuilder (specialization of abstract MarketDataBuilder) to build market data and return it as a CurveData object. Then, Client is creating and using USDLiborZeroCurve object (implementation of ICurve interface). Inside USDLiborZeroCurve object, the curves (discount factors and forward rates) are going to be created in different types of sequential calculations (data checking, data selection, bootstrapping, interpolation). Finally, Client is requesting ExcelPrinter (specialization of abstract CurvePrinter) object to print the resulting dataset into Excel worksheet. During the calculation process, USDLiborZeroCurve is also using two static classes (DateConvention and Interpolation, not shown in UML diagram), which contains different types of date-related and interpolation-related methods.

























THE PROGRAM

The part of the program, which is creating the outcome curves (discount factors and forward rates), is shown below. Create a new Class project and just copyPaste everything into a new cs-file.


using System;
using System.Collections;
using System.Collections.Generic;
using System.Linq;
using ExcelDna.Integration;

namespace CurveAddin
{
    //
    // public enumerators and delegate methods
    public enum ENUM_MARKETDATATYPE { CASH = 1, FRA = 2, FUTURE = 3, SWAP = 4, DF = 5, SPOTRATE = 6, FORWARDRATE = 7 }
    public enum ENUM_PERIODTYPE { MONTHS = 1, YEARS = 2 }
    public delegate double DayCountFactor(DateTime start, DateTime end);
    public delegate DateTime AddPeriod(DateTime start, ENUM_PERIODTYPE periodType, int period);
    public delegate double Interpolator(DayCountFactor dayCountFactor, Dictionary<DateTime, double> data, DateTime key);
    public delegate double ConvexityAdjustment(double rateVolatility, double start, double end);
    //
    //
    // class hierarchy for curve printer
    public abstract class CurvePrinter
    {
        public abstract void Print();
    }
    public class ExcelPrinter : CurvePrinter
    {
        private static dynamic Excel;
        private dynamic[,] data;
        public ExcelPrinter(dynamic[,] data)
        {
            this.data = data;
        }
        public override void Print()
        {
            // Create Excel application
            Excel = ExcelDnaUtil.Application;
            //
            // clear old data from output range, resize output range 
            // and finally print data to Excel worksheet
            Excel.Range["_USDLiborZeroCurve"].CurrentRegion = "";
            Excel.Range["_USDLiborZeroCurve"].Resize[data.GetLength(0), data.GetLength(1)] = data;
        }
    }
    //
    //
    // class hierarchy for market data builder
    public abstract class MarketDataBuilder
    {
        public abstract CurveData Build();
    }
    public class ExcelBuilder : MarketDataBuilder
    {
        private static dynamic Excel;
        private DateTime settlementDate;
        public DateTime SettlementDate { get { return this.settlementDate; } }
        //
        public override CurveData Build()
        {
            // Create Excel application
            Excel = ExcelDnaUtil.Application;
            //
            // read settlement date from Excel worksheet
            settlementDate = DateTime.FromOADate((double)Excel.Range("_settlementDate").Value2);
            //
            // read source security data from Excel worksheet
            object[,] ExcelSourceData = (object[,])Excel.Range["_marketData"].CurrentRegion.Value2;
            //
            // create curve data object from source security data
            CurveData marketData = new CurveData(Interpolation.LinearInterpolation);
            int rows = ExcelSourceData.GetUpperBound(0);
            for (int i = 1; i <= rows; i++)
            {
                DateTime maturity = DateTime.FromOADate((double)ExcelSourceData[i, 1]);
                double rate = (double)ExcelSourceData[i, 2];
                string instrumentType = ((string)ExcelSourceData[i, 3]).ToUpper();
                ENUM_MARKETDATATYPE marketDataType = (ENUM_MARKETDATATYPE)Enum.Parse(typeof(ENUM_MARKETDATATYPE), instrumentType);
                marketData.AddCurveDataElement(maturity, rate, marketDataType);
            }
            return marketData;
        }
    }
    //
    //
    // interface for all curve objects
    public interface ICurve
    {
        void Create();
        double GetDF(DateTime maturity);
        double GetFWD(DateTime start);
        Dictionary<DateTime, double> GetDF(DateTime start, int nYears);
        Dictionary<DateTime, double> GetFWD(DateTime start, int nYears);
        CurveData DiscountCurve { get; }
        CurveData ForwardCurve { get; }
    }
    //
    // implementation for USD Libor curve
    public class USDLiborZeroCurve : ICurve
    {
        public readonly DayCountFactor dayCountFactor;
        public readonly AddPeriod addPeriod;
        public readonly int basis;
        public readonly Interpolator interpolator;
        public readonly DateTime settlementDate;
        public CurveData DiscountCurve { get { return this.discountCurve; } }
        public CurveData ForwardCurve { get { return this.forwardCurve; } }
        //
        private CurveData marketData;
        private CurveData curveDataSelection;
        private CurveData bootstrapCurve;
        private CurveData spotCurve;
        private CurveData discountCurve;
        private CurveData forwardCurve;
        private int nCash;
        private int nFuturesOrFRAs;
        private bool adjustmentForConvexity;
        private ConvexityAdjustment convexityAdjustment;
        private double rateVolatility;
        //
        public USDLiborZeroCurve(CurveData marketData, Interpolator interpolator, AddPeriod addPeriod,
            DayCountFactor dayCountFactor, DateTime settlementDate, int nCash, int nFuturesOrFRAs,
            bool adjustmentForConvexity = false, ConvexityAdjustment convexityAdjustment = null, double rateVolatility = 0.0)
        {
            this.marketData = marketData;
            this.interpolator = interpolator;
            this.addPeriod = addPeriod;
            this.dayCountFactor = dayCountFactor;
            this.settlementDate = settlementDate;
            this.nCash = nCash;
            this.nFuturesOrFRAs = nFuturesOrFRAs;
            this.basis = 3; // HARD-CODED !! for USD Libor curve
            this.adjustmentForConvexity = adjustmentForConvexity; // optional parameter
            this.convexityAdjustment = convexityAdjustment; // optional parameter
            this.rateVolatility = rateVolatility; // optional parameter
        }
        public void Create()
        {
            // sequence of private methods for creating spot discount curve and
            // simply-compounded forward rate curve for a given set of market data
            checkMarketData();
            selectCurveData();
            bootstrapDiscountFactors();
            createSpotCurve();
            createDiscountCurve();
            createForwardCurve();
        }
        // get discount factor for a given maturity date
        public double GetDF(DateTime maturity)
        {
            return discountCurve.GetMarketRate(ENUM_MARKETDATATYPE.DF, maturity, dayCountFactor);
        }
        // get dictionary consisting of date and discount factor for a date schedule
        public Dictionary<DateTime, double> GetDF(DateTime start, int nYears)
        {
            List<DateTime> schedule = DateConvention.CreateDateSchedule(start, nYears, basis, ENUM_PERIODTYPE.MONTHS, addPeriod);
            Dictionary<DateTime, double> curve = new Dictionary<DateTime, double>();
            schedule.ForEach(it => curve.Add(it, GetDF(it)));
            return curve;
        }
        // get simply-compounded forward rate for a given start date
        public double GetFWD(DateTime start)
        {
            return forwardCurve.GetMarketRate(ENUM_MARKETDATATYPE.FORWARDRATE, start, dayCountFactor);
        }
        // get dictionary consisting of date and simply-compounded forward rate for a date schedule
        public Dictionary<DateTime, double> GetFWD(DateTime start, int nYears)
        {
            List<DateTime> schedule = DateConvention.CreateDateSchedule(start, nYears, basis, ENUM_PERIODTYPE.MONTHS, addPeriod);
            Dictionary<DateTime, double> curve = new Dictionary<DateTime, double>();
            schedule.ForEach(it => curve.Add(it, GetFWD(it)));
            return curve;
        }
        // use interpolated spot discount factor curve for calculating
        // simply-compounded forward rates for all required maturities (basis)
        // note : maturity element of the forward curve stores the information
        // on when the 3-month period starts for a given forward rate element
        private void createForwardCurve()
        {
            forwardCurve = new CurveData(interpolator);
            int n = discountCurve.Count();
            DateTime maturity;
            double dt = 0.0;
            double fdf = 0.0;
            double f = 0.0;
            //
            for (int i = 0; i < n; i++)
            {
                if (i == 0)
                {
                    // first forward rate is the first spot rate
                    maturity = discountCurve[i].MaturityDate;
                    fdf = discountCurve[i].Rate;
                    dt = dayCountFactor(settlementDate, maturity);
                    f = ((1 / fdf) - 1) / dt;
                    forwardCurve.AddCurveDataElement(settlementDate, f, ENUM_MARKETDATATYPE.FORWARDRATE);
                }
                else
                {
                    // other forward rates are calculated recursively
                    // from previous spot discount factors
                    maturity = discountCurve[i].MaturityDate;
                    DateTime previousMaturity = discountCurve[i - 1].MaturityDate;
                    fdf = discountCurve[i].Rate / discountCurve[i - 1].Rate;
                    dt = dayCountFactor(previousMaturity, maturity);
                    f = ((1 / fdf) - 1) / dt;
                    forwardCurve.AddCurveDataElement(previousMaturity, f, ENUM_MARKETDATATYPE.FORWARDRATE);
                }
            }
        }
        // use continuously compounded spot rate curve for interpolating 
        // continuously compounded spot rates for all required maturities 
        // and convert these spot rates back to spot discount factors
        private void createDiscountCurve()
        {
            discountCurve = new CurveData(interpolator);
            DateTime finalCurveDate = spotCurve.ElementAt(spotCurve.Count() - 1).MaturityDate;
            DateTime t;
            int counter = 0;
            double dt = 0.0;
            double r = 0.0;
            double df = 0.0;
            //
            do
            {
                counter++;
                t = addPeriod(settlementDate, ENUM_PERIODTYPE.MONTHS, basis * counter);
                dt = dayCountFactor(settlementDate, t);
                r = spotCurve.GetMarketRate(ENUM_MARKETDATATYPE.SPOTRATE, t, dayCountFactor);
                df = Math.Exp(-r * dt);
                discountCurve.AddCurveDataElement(t, df, ENUM_MARKETDATATYPE.DF);
            } while (t < finalCurveDate);
        }
        // create continuously compounded spot rate curve 
        // from bootstrapped discount factors
        private void createSpotCurve()
        {
            spotCurve = new CurveData(interpolator);
            double t = 0.0;
            double r = 0.0;
            int n = bootstrapCurve.Count();
            for (int i = 0; i < n; i++)
            {
                t = dayCountFactor(settlementDate, bootstrapCurve.ElementAt(i).MaturityDate);
                r = -Math.Log(bootstrapCurve.ElementAt(i).Rate) / t;
                spotCurve.AddCurveDataElement(bootstrapCurve.ElementAt(i).MaturityDate, r, ENUM_MARKETDATATYPE.SPOTRATE);
            }
        }
        // use bootstrap algorithm to create spot discount factors
        // from all selected curve data elements
        private void bootstrapDiscountFactors()
        {
            bootstrapCurve = new CurveData(interpolator);
            double dt = 0.0;
            double r = 0.0;
            double df = 0.0;
            double Q = 0.0;
            int n = curveDataSelection.Count();
            //
            for (int i = 0; i < n; i++)
            {
                if (curveDataSelection[i].InstrumentType == ENUM_MARKETDATATYPE.CASH)
                {
                    dt = dayCountFactor(settlementDate, curveDataSelection[i].MaturityDate);
                    r = curveDataSelection[i].Rate;
                    df = 1 / (1 + r * dt);
                    bootstrapCurve.AddCurveDataElement(curveDataSelection[i].MaturityDate, df, ENUM_MARKETDATATYPE.DF);
                }
                if ((curveDataSelection[i].InstrumentType == ENUM_MARKETDATATYPE.FRA) |
                    (curveDataSelection[i].InstrumentType == ENUM_MARKETDATATYPE.FUTURE))
                {
                    dt = dayCountFactor(curveDataSelection[i - 1].MaturityDate, curveDataSelection[i].MaturityDate);
                    r = curveDataSelection[i].Rate;
                    df = bootstrapCurve.ElementAt(i - 1).Rate / (1 + r * dt);
                    bootstrapCurve.AddCurveDataElement(curveDataSelection[i].MaturityDate, df, ENUM_MARKETDATATYPE.DF);
                    //
                    if ((curveDataSelection[i + 1].InstrumentType == ENUM_MARKETDATATYPE.SWAP))
                        Q += bootstrapCurve.ElementAt(i).Rate * dayCountFactor(settlementDate, curveDataSelection[i].MaturityDate);
                }
                //
                if (curveDataSelection[i].InstrumentType == ENUM_MARKETDATATYPE.SWAP)
                {
                    r = curveDataSelection[i].Rate;
                    dt = dayCountFactor(bootstrapCurve.ElementAt(i - 1).MaturityDate, curveDataSelection[i].MaturityDate);
                    df = (1 - r * Q) / (r * dt + 1);
                    bootstrapCurve.AddCurveDataElement(curveDataSelection[i].MaturityDate, df, ENUM_MARKETDATATYPE.DF);
                    Q += (df * dt);
                }
            }
        }
        // select rate instruments to be used from a given set of curve data elements
        private void selectCurveData()
        {
            curveDataSelection = new CurveData(interpolator);
            int counter = 0;
            double rate = 0.0;
            DateTime maturityDate;
            //
            // select cash securities
            for (int i = 1; i <= nCash; i++)
            {
                counter++;
                maturityDate = addPeriod(settlementDate, ENUM_PERIODTYPE.MONTHS, basis * counter);
                // check if cash rate for required maturity exists
                if (!marketData.ElementLookup(ENUM_MARKETDATATYPE.CASH, maturityDate))
                    throw new Exception("USDLiborZeroCurve error : required cash securities are missing");
                rate = marketData.GetMarketRate(ENUM_MARKETDATATYPE.CASH, maturityDate, dayCountFactor);
                curveDataSelection.AddCurveDataElement(maturityDate, rate, ENUM_MARKETDATATYPE.CASH);
            }
            // select fra or futures contracts
            if (marketData.ElementLookup(ENUM_MARKETDATATYPE.FRA))
            {
                for (int i = 1; i <= nFuturesOrFRAs; i++)
                {
                    if (i > 1) counter++;
                    maturityDate = addPeriod(settlementDate, ENUM_PERIODTYPE.MONTHS, basis * counter);
                    // check if fra rate for required maturity exists
                    if (!marketData.ElementLookup(ENUM_MARKETDATATYPE.FRA, maturityDate))
                        throw new Exception("USDLiborZeroCurve error : required FRA contracts are missing");
                    rate = marketData.GetMarketRate(ENUM_MARKETDATATYPE.FRA, maturityDate, dayCountFactor);
                    curveDataSelection.AddCurveDataElement(addPeriod(maturityDate, ENUM_PERIODTYPE.MONTHS, basis), rate, ENUM_MARKETDATATYPE.FRA);
                }
            }
            else
            {
                for (int i = 1; i <= nFuturesOrFRAs; i++)
                {
                    if (i > 1) counter++;
                    maturityDate = addPeriod(settlementDate, ENUM_PERIODTYPE.MONTHS, basis * counter);
                    // check if implied futures rate for required maturity exists
                    if (!marketData.ElementLookup(ENUM_MARKETDATATYPE.FUTURE, maturityDate))
                        throw new Exception("USDLiborZeroCurve error : required futures contracts are missing");
                    rate = marketData.GetMarketRate(ENUM_MARKETDATATYPE.FUTURE, maturityDate, dayCountFactor);
                    //
                    // forward rate = futures rate - convexity adjustment
                    if (adjustmentForConvexity)
                    {
                        double t1 = dayCountFactor(settlementDate, maturityDate);
                        double t2 = t1 + (basis / 12.0);
                        rate -= convexityAdjustment(rateVolatility, t1, t2);
                    }
                    curveDataSelection.AddCurveDataElement(addPeriod(maturityDate, ENUM_PERIODTYPE.MONTHS, basis), rate, ENUM_MARKETDATATYPE.FUTURE);
                }
            }
            // select swap contracts
            DateTime lastSwapYear = marketData[marketData.Count() - 1].MaturityDate;
            DateTime lastFRAOrFutureYear = curveDataSelection[curveDataSelection.Count() - 1].MaturityDate;
            int nSwaps = (lastSwapYear.Year - lastFRAOrFutureYear.Year);
            for (int i = 1; i <= nSwaps; i++)
            {
                counter++;
                maturityDate = addPeriod(settlementDate, ENUM_PERIODTYPE.YEARS, i + 1);
                // check if swap rate for required maturity exists
                if (!marketData.ElementLookup(ENUM_MARKETDATATYPE.SWAP, maturityDate))
                    throw new Exception("USDLiborZeroCurve error : required swap contracts are missing");
                rate = marketData.GetMarketRate(ENUM_MARKETDATATYPE.SWAP, maturityDate, dayCountFactor);
                curveDataSelection.AddCurveDataElement(maturityDate, rate, ENUM_MARKETDATATYPE.SWAP);
            }
        }
        // rough diagnostics : check for completely non-existing market data
        // requirement : all three rate categories (cash, FRA/futures, swaps) 
        // must be provided by the client in order to create the curves
        private void checkMarketData()
        {
            // cash securities
            if (!marketData.ElementLookup(ENUM_MARKETDATATYPE.CASH))
                throw new Exception("LiborZeroCurve error : cash securities are required to build the curve");
            //
            // fra/futures contracts
            if ((!marketData.ElementLookup(ENUM_MARKETDATATYPE.FUTURE)) && (!marketData.ElementLookup(ENUM_MARKETDATATYPE.FRA)))
                throw new Exception("LiborZeroCurve error : FRA or futures contracts are required to build the curve");
            //
            // swap contracts
            if (!marketData.ElementLookup(ENUM_MARKETDATATYPE.SWAP))
                throw new Exception("LiborZeroCurve error : swap contracts are required to build the curve");
        }
    }
    //
    //
    // container class for holding multiple curve data elements
    public class CurveData : IEnumerable<CurveDataElement>
    {
        private List<CurveDataElement> curveDataElements;
        private Interpolator interpolator;
        //
        public CurveData(Interpolator interpolator)
        {
            this.interpolator = interpolator;
            curveDataElements = new List<CurveDataElement>();
        }
        public void AddCurveDataElement(DateTime maturity, double rate, ENUM_MARKETDATATYPE instrumentType)
        {
            curveDataElements.Add(new CurveDataElement(maturity, rate, instrumentType));
        }
        public void AddCurveDataElement(CurveDataElement curveDataElement)
        {
            curveDataElements.Add(curveDataElement);
        }
        // implementation for generic IEnumerable
        public IEnumerator<CurveDataElement> GetEnumerator()
        {
            foreach (CurveDataElement curveDataElement in curveDataElements)
            {
                yield return curveDataElement;
            }
        }
        // implementation for non-generic IEnumerable
        IEnumerator IEnumerable.GetEnumerator()
        {
            return GetEnumerator();
        }
        // read-only indexer
        public CurveDataElement this[int index]
        {
            get
            {
                return curveDataElements[index];
            }
        }
        public double GetMarketRate(ENUM_MARKETDATATYPE instrumentType, DateTime maturity, DayCountFactor dayCountFactor)
        {
            // filter required market data elements by instrument type
            List<CurveDataElement> group = curveDataElements.Where(it => it.InstrumentType == instrumentType).ToList<CurveDataElement>();
            //
            // extract maturity and rate into dictionary object
            Dictionary<DateTime, double> data = new Dictionary<DateTime, double>();
            group.ForEach(it => data.Add(it.MaturityDate, it.Rate));
            //
            // get market rate for a given date by using given interpolation delegate method
            return interpolator(dayCountFactor, data, maturity);
        }
        // check if elements with specific instrument type and maturity exists
        public bool ElementLookup(ENUM_MARKETDATATYPE instrumentType, DateTime maturity)
        {
            // first, filter required market data elements
            List<CurveDataElement> group = curveDataElements.Where(it => it.InstrumentType == instrumentType).ToList<CurveDataElement>();
            //
            // then, check if maturity lies between min and max maturity of filtered group
            bool hasElement = ((maturity >= group.Min(it => it.MaturityDate)) && (maturity <= group.Max(it => it.MaturityDate)));
            return hasElement;
        }
        // check if elements with only specific instrument type exists
        public bool ElementLookup(ENUM_MARKETDATATYPE instrumentType)
        {
            int elements = curveDataElements.Count(it => it.InstrumentType == instrumentType);
            bool hasElements = false;
            if (elements > 0) hasElements = true;
            return hasElements;
        }
    }
    //
    //
    // class holding information on one curve data element
    public class CurveDataElement
    {
        private DateTime maturityDate;
        private double rate;
        private ENUM_MARKETDATATYPE rateType;
        //
        public DateTime MaturityDate { get { return this.maturityDate; } }
        public double Rate { get { return this.rate; } }
        public ENUM_MARKETDATATYPE InstrumentType { get { return this.rateType; } }
        //
        public CurveDataElement(DateTime maturity, double rate, ENUM_MARKETDATATYPE rateType)
        {
            this.maturityDate = maturity;
            this.rate = rate;
            this.rateType = rateType;
        }
    }
    //
    //
    // static library class for handling date-related convention calculations
    public static class DateConvention
    {
        // calculate time difference between two dates by using ACT/360 convention
        public static double ACT360(DateTime start, DateTime end)
        {
            return (end - start).TotalDays / 360;
        }
        // create a list of scheduled dates for a given basis and date convention
        public static List<DateTime> CreateDateSchedule(DateTime start, int nYears, int basis,
            ENUM_PERIODTYPE periodType, AddPeriod addPeriod)
        {
            List<DateTime> schedule = new List<DateTime>();
            int nPeriods = nYears * (12 / basis);
            for (int i = 1; i <= nPeriods; i++)
            {
                schedule.Add(addPeriod(start, periodType, (basis * i)));
            }
            return schedule;
        }
        // add period into a given date by using modified following convention
        public static DateTime AddPeriod_ModifiedFollowing(DateTime start, ENUM_PERIODTYPE periodType, int period)
        {
            DateTime dt = new DateTime();
            //
            switch (periodType)
            {
                case ENUM_PERIODTYPE.MONTHS:
                    dt = start.AddMonths(period);
                    break;
                case ENUM_PERIODTYPE.YEARS:
                    dt = start.AddYears(period);
                    break;
            }
            //
            switch (dt.DayOfWeek)
            {
                case DayOfWeek.Saturday:
                    dt = dt.AddDays(2.0);
                    break;
                case DayOfWeek.Sunday:
                    dt = dt.AddDays(1.0);
                    break;
            }
            return dt;
        }
        // calculate value for convexity adjustment for a given time period
        public static double SimpleConvexityApproximation(double rateVolatility, double start, double end)
        {
            return 0.5 * (rateVolatility * rateVolatility * start * end);
        }
    }
    //
    //
    // static library class for storing interpolation methods to be used by delegates
    public static class Interpolation
    {
        public static double LinearInterpolation(DayCountFactor dayCountFactor,
            Dictionary<DateTime, double> data, DateTime 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) throw new Exception("Interpolation error : lower bound violation");
                if (key > data.ElementAt(data.Count - 1).Key) throw new Exception("Interpolation error : upper bound violation");
            }
            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 = dayCountFactor(data.ElementAt(i).Key, data.ElementAt(i + 1).Key);
                        double w = dayCountFactor(data.ElementAt(i).Key, key) / t;
                        value = data.ElementAt(i).Value * (1 - w) + data.ElementAt(i + 1).Value * w;
                        break;
                    }
                }
            }
            return value;
        }
    }
}


EXCEL INTERFACING

For interfacing the previous C# program with Excel, carefully follow all the instructions given in this blog post. Add another Class module (a new cs-file) into this project and copyPaste the following program into this file. This is the program (CreateUSDLiborCurves), which will be called from our Excel worksheet.


using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Windows.Forms;
using ExcelDna.Integration;

namespace CurveAddin
{
    public static class CurveAddin
    {
        public static void CreateUSDLiborCurves()
        {
            try
            {
                // build market data from Excel worksheet
                MarketDataBuilder builder = new ExcelBuilder();
                CurveData marketData = builder.Build();
                DateTime settlementDate = ((ExcelBuilder)builder).SettlementDate;
                //
                // construct USD Libor curve object
                // HARD-CODED parameters : 
                // interpolation method, date conventions, number of contracts for cash (n=1) and futures (n=3)
                ICurve curve = new USDLiborZeroCurve(marketData, Interpolation.LinearInterpolation, 
                    DateConvention.AddPeriod_ModifiedFollowing, DateConvention.ACT360, settlementDate, 1, 3);
                curve.Create();
                //
                // read discount factor and forward rate data into 2d-array
                int rows = curve.DiscountCurve.Count();
                int cols = 3;
                dynamic[,] data = new dynamic[rows, cols];
                for (int i = 0; i < rows; i++)
                {
                    data[i, 0] = curve.DiscountCurve[i].MaturityDate.ToOADate();
                    data[i, 1] = curve.DiscountCurve[i].Rate;
                    data[i, 2] = curve.ForwardCurve[i].Rate;
                }
                //
                // print curve data into Excel worksheet
                (new ExcelPrinter(data)).Print();
            }
            catch (Exception e)
            {
                MessageBox.Show(e.Message);
            }
        }
    }
}



EXCEL WORKSHEET SETTINGS

Prepare the following market data (Flavell's book and Excel workbook on chapter three) and named ranges (marked with yellow color) into Excel worksheet. The program will take settlement date and market data range (date, rate, rate type) as input parameters. Note, that some of the parameters needed to create the curves have been hard-coded in the program (interpolation method, date conventions, number of contracts for cash and futures). However, it should be fairly straightforward to include all these parameters to be read directly from Excel worksheet. Finally, insert a form control button into worksheet and assign a macro for it (CreateUSDLiborCurves).
















Finally, a couple of notes concerning the market data setting. When creating Curve object, Client has to provide number of cash securities and number futures/FRA contracts as input parameters in constructor. Example program has been hard-coded to use one cash security and three futures contracts and it will hereby use Libor swap contracts starting on the year two. Now, for setting market data for this program, there are some general rules. First, the latest maturity for cash securities has to be later than the first maturity for futures/FRA contracts. Also, the last future/FRA maturity has to be later than the first swap contract maturity date minus one year.

TEST RUN

After I press button in my Excel, I will get the following results (columns H to J) from the program. The range consists of dates (quarterly up to 8.2.2038), discount factors and simply-compounded Libor forward rates for USD. The first result row should be interpreted as follows : discount factor (0.99220) gives factor for discounting three-month cash flow. Simply-compounded Libor forward rate (3.1450 %) gives three-month forward rate for a period, which is starting on 6.2.2008 and ending 6.5.2008. Similarly, Libor forward rate (2.7837 %) gives three-month forward rate for a period, which is starting on 6.5.2008 and ending 6.8.2008.



















After this, using generated curves (discount factors and forward rates) is straightforward. As an example, I have calculated PV for 2-year cap on 3-month USD Libor. After creating date schedule for this security, getting forward rates and discount factors can be requested by using familiar Excel worksheet functions.


















AFTERTHOUGHTS

The procedure of creating discount factors and Libor forward rates programmatically in C#, has been fully opened in this blog post. Source market data and creation procedures are following examples taken from the book written by Richard Flavell. With the tools presented in this blog, one should also be able to interface this program with Excel, if so desired.

I would like to thank Govert Van Drimmelen again for his amazing Excel-DNA, what I am always using for interfacing my C# program with Excel. For learning more things about Excel-DNA, check out its homepage. Getting more information and examples with your problems, the main source is Excel-DNA google group. Remember also, that Excel-DNA is an open-source project, and we (the happy users) can invest its future development by making a donation.

Finally, Thanks for spending your precious time in here and reading my blog.

-Mike Juniperhill