Showing posts with label Discount factors. Show all posts
Showing posts with label Discount factors. Show all posts

Wednesday, October 30, 2019

QuantLib-Python: Note on ForwardCurve Construction

In this post, we will construct QuantLib ForwardCurve instance and investigate the resulting term structure of discount factors. Python program and Excel can be downloaded from my GitHub page. First, let us assume the following market data containing [A] one cash deposit and three interpolated rates from futures strip.


















After completing bootstrapping procedure, we should get [B] simple spot discount factors and [C] simple forward rates. Next, let us implement this scheme in QuantLib. Import QuantLib libraries and create relevant date objects. Then, set forward rates (from [C]), dates and create QuantLib forward curve instance.

import QuantLib as ql
today = ql.Date(30, 10, 2019)
ql.Settings.instance().evaluationDate = today  
settlementDate = ql.TARGET().advance(today, ql.Period(2, ql.Days))

rts = [0.03145, 0.03145, 0.0278373626373627, 0.0253076923076923, 0.0249373626373629]
dts = [settlementDate, ql.Date(1,2,2020), ql.Date(1,5,2020), ql.Date(1,8,2020), ql.Date(1,11,2020)]
c = ql.ForwardCurve(dts, rts, ql.Actual360(), ql.NullCalendar(), ql.BackwardFlat())
df = [c.discount(d) for d in dts]
print(df)

We get the following set of spot discount factors from the curve instance.
[1.0, 0.9919949898918935, 0.9851153255552918, 0.9787646299045335, 0.9725469122728971]
Note, that these are not what we expected to see. "The correct ones" are calculated in the column [B] in the example above. The issue here is, that when constructing ForwardCurve object, the given forward rates inputs are implicitly defined as continuously compounded rates. Assuming we would like to get the same discount factor term structure as in our example above in [B], we need to convert our forward rate inputs from simple to continuously compounded rates. Let us implement this in QuantLib. Convert simple rates to continuous rates. Then, reset forward rates and re-create QuantLib forward curve instance.

for i in range(len(rts) - 1):
    r_simple = ql.InterestRate(rts[i + 1], ql.Actual360(), ql.Simple, ql.Once)
    t = ql.Actual360().yearFraction(dts[i], dts[i + 1])
    r_continuous = r_simple.equivalentRate(ql.Continuous, ql.NoFrequency, t)
    rts[i + 1] = r_continuous.rate()
    
# set rate for the first node
rts[0] = rts[1]
c = ql.ForwardCurve(dts, rts, ql.Actual360(), ql.NullCalendar(), ql.BackwardFlat())
df = [c.discount(d) for d in dts]
print(df)

We get the following "correct" set of spot discount factors from the curve instance.
[1.0, 0.9920268596783518, 0.9851707210231435, 0.9788400520709886, 0.9726415228427708]
Inputs (dates and forward rates) for ForwardCurve object are defined in yellow box in column [E] below.
















ForwardCurve object interpolation scheme is backward flat. This means, that for the final period (from 1-Aug-20 to 1-Nov-20), forward rate of 2.48582% will be applied. Implicitly this means, that there is no practical use for the first given forward rate (at the node, where the date is 1-Nov-19). However, this rate has to be defined anyway.

Note, that QuantLib logic of constructing ForwardCurve object is completely correct. One should just be aware of the fact, that input forward rates given are handled implicitly as continuously compounded rates.

Finally, there are some other types of issues concerning forward rate interpolation, which might be worth of knowing when working with ForwardCurve object. These issues are fully explained in this video given by QuantLib lead developer Luigi Ballabio.

As always, thanks for reading this blog.
-Mike

Tuesday, December 4, 2018

QuantLib-Python: Term Structure Simulation Using HW1F Model

This post is presenting Python program, which uses QuantLib tools for simulating yield term structure for the chosen one-factor interest rate model. Further comparison results are also showing, that simulation method is able to replicate the initial yield curve, without any notable differences.

The idea is this: create yield curve object by using current market data (flat forward) and 1-D stochastic process for short rate dynamics (Hull-White). Then, use separate method (GeneratePaths) for generating desired amount of paths for a chosen stochastic process. Next, integrate short-rate for all simulated paths and calculate average zero-coupon bond prices. Finally, create a new yield curve object by using previously simulated zero-coupon bond prices and compare the resulting set of discount factors with the ones requested from the original yield curve. It should be noted, that there is a separate class (Grid), which is used for hosting all schedule-related information (such as schedule, dates and times) and their conversions (from schedule to times, from schedule to dates) in one compact place.

From the screenshot below we can conclude, that as
  • the data and the other parameters are within "sensible ranges", 
  • the number of paths is large enough and 
  • discretization error is minimized by selecting a small enough time step, 

simulated yield curve is able to replicate the initial yield curve without any notable differences. Some relevant issues around this particular topic has been chewed in here and here. The same stuff (and a lot more) has also been published in QuantLib Python Cookbook by the blog author Gouthaman Balaraman and QuantLib lead developer Luigi Ballabio.

Finally, outside of being a nice QuantLib exercise itself, there is not much point to simulate zero-coupon bond prices. Needless to say, the essence of Monte Carlo method (simulate a path, create term structure from it, price a product) can be used for much more interesting valuation problems.

Thanks for reading my blog.
-Mike



































%config IPCompleter.greedy = True
import math as Math
from QuantLib import *
import numpy as Numpy
import matplotlib.pyplot as Matplotlib

def main():

    # create grid object for 30Y, having time step of 1 day
    startDate = Date(3, December, 2018)
    endDate = Date(3, December, 2048)
    tenor = Period(1, Days)
    grid = Grid(startDate, endDate, tenor)

    # create yield curve and Hull-White one-factor interest rate model
    curve = YieldTermStructureHandle(FlatForward(startDate, 0.04875825, Actual365Fixed()))
    reversionSpeed = 0.05
    rateVolatility = 0.00586
    process = HullWhiteProcess(curve, reversionSpeed, rateVolatility)

    # request paths from generator method
    nPaths = 25000
    paths = GeneratePaths(process, grid.GetTimeGrid(), nPaths)

    # container for simulated zero-coupon bonds
    zeros = Numpy.zeros(shape = (grid.GetSize()))
    dt = grid.GetDt()
    gridSize = grid.GetSize()

    # process short-rate path integrations for all simulated paths
    for i in range(nPaths):
        integral = 0.0
        for j in range(gridSize):
            integral = integral + paths[i, j]
            if(j == 0):
                # zero-coupon bond price today is 1.0
                zeros[j] = 1.0 * nPaths
            else:
                zeros[j] = zeros[j] + Math.exp(-integral * dt)

    # calculate averages for all simulated zero-coupon bond prices
    zeros = zeros / nPaths

    # create yield term structure object from simulated bond prices
    times = grid.GetTimes()
    dates = grid.GetDates()
    simulatedCurve = DiscountCurve(dates, zeros, Actual365Fixed(), NullCalendar())

    # get discount factors for simulated and original yield curves
    dfs = Numpy.zeros(shape = (gridSize))
    simulatedDfs = Numpy.zeros(shape = (gridSize))
    for i in range(gridSize):
        simulatedDfs[i] = simulatedCurve.discount(times[i])
        dfs[i] = curve.discount(times[i])

    # plot simulated and original discount factors
    Matplotlib.title('discount factors')
    Matplotlib.plot(times, simulatedDfs, linestyle = 'dashed', label = 'simulated curve')
    Matplotlib.plot(times, dfs, linestyle = 'solid', label = 'original curve')
    Matplotlib.legend()
    Matplotlib.show()

    # plot difference between simulated and original discount factors in basis points
    Matplotlib.title('difference (bps)')
    Matplotlib.plot(times, (dfs - simulatedDfs) * 10000)
    Matplotlib.show()

    
# path generator method for uncorrelated and correlated 1-D stochastic processes
def GeneratePaths(process, timeGrid, n):

    # correlated processes, use GaussianMultiPathGenerator
    if(type(process) == StochasticProcessArray):
        times = []; [times.append(timeGrid[t]) for t in range(len(timeGrid))]        
        nGridSteps = (len(times) - 1) * process.size()
        sequenceGenerator = UniformRandomSequenceGenerator(nGridSteps, UniformRandomGenerator())
        gaussianSequenceGenerator = GaussianRandomSequenceGenerator(sequenceGenerator)
        pathGenerator = GaussianMultiPathGenerator(process, times, gaussianSequenceGenerator, False)        
        paths = Numpy.zeros(shape = (n, process.size(), len(timeGrid)))
        
        # loop through number of paths
        for i in range(n):
            # request multiPath, which contains the list of paths for each process
            multiPath = pathGenerator.next().value()
            # loop through number of processes
            for j in range(multiPath.assetNumber()):
                # request path, which contains the list of simulated prices for a process
                path = multiPath[j]
                # push prices to array
                paths[i, j, :] = Numpy.array([path[k] for k in range(len(path))])
        # resulting array dimension: n, process.size(), len(timeGrid)
        return paths

    # uncorrelated processes, use GaussianPathGenerator
    else:
        sequenceGenerator = UniformRandomSequenceGenerator(len(timeGrid), UniformRandomGenerator())
        gaussianSequenceGenerator = GaussianRandomSequenceGenerator(sequenceGenerator)
        maturity = timeGrid[len(timeGrid) - 1]
        pathGenerator = GaussianPathGenerator(process, maturity, len(timeGrid), gaussianSequenceGenerator, False)
        paths = Numpy.zeros(shape = (n, len(timeGrid)))
        for i in range(n):
            path = pathGenerator.next().value()
            paths[i, :] = Numpy.array([path[j] for j in range(len(timeGrid))])
        # resulting array dimension: n, len(timeGrid)
        return paths


# class for hosting schedule-related information (dates, times)
class Grid:
    def __init__(self, startDate, endDate, tenor):
        # create date schedule, ignore conventions and calendars
        self.schedule = Schedule(startDate, endDate, tenor, NullCalendar(), 
            Unadjusted, Unadjusted, DateGeneration.Forward, False)
        self.dayCounter = Actual365Fixed()
    def GetDates(self):
        # get list of scheduled dates
        dates = []
        [dates.append(self.schedule[i]) for i in range(self.GetSize())]
        return dates
    def GetTimes(self):
        # get list of scheduled times
        times = []
        [times.append(self.dayCounter.yearFraction(self.schedule[0], self.schedule[i])) 
            for i in range(self.GetSize())]
        return times
    def GetMaturity(self):
        # get maturity in time units
        return self.dayCounter.yearFraction(self.schedule[0], self.schedule[self.GetSteps()])
    def GetSteps(self):
        # get number of steps in schedule
        return self.GetSize() - 1    
    def GetSize(self):
        # get total number of items in schedule
        return len(self.schedule)    
    def GetTimeGrid(self):
        # get QuantLib TimeGrid object, constructed by using list of scheduled times
        return TimeGrid(self.GetTimes(), self.GetSize())
    def GetDt(self):
        # get constant time step
        return self.GetMaturity() / self.GetSteps()
    
main()

Sunday, June 3, 2018

QuantLib : Bootstrapping OIS curve

As we all are aware, all collateralized derivative contracts must be valued by using dual-curve discounting, by using separate curves for cash flow projection and cash flow discounting. While this post does not yet explain how to implement such dual-curve disconting scheme in QuantLib, it shows how to use QuantLib for bootstrapping required OIS discount curve from available market rates. Some further reading material, which nicely explains the justification for the use of dual-curve discounting, is available in here.

The program


First, for the sake of brevity, required market data (Eonia rates) has been hard-coded into map container. Then, container will be iterated through for constructing all required rate helpers. After this, piecewise yield term structure will be created by using constructed rate helpers. Finally, an equally-spaced time grid (monthly) will be created and discount factors are requested from bootstrapped Eonia curve for each grid time point (screenshot below). The curve itself (eurOisCurve) could now be used as a discounting curve in a dual-curve valuation scheme in QuantLib.

















Thanks for reading this blog.
-Mike

#include <iostream>
#include <ql\quantlib.hpp>
#include <algorithm>
#include <map>
using namespace QuantLib;

int main() {

 try {
  // create common data 
  Date today(7, Jul, 2017);
  DayCounter dayCounter = Actual360();
  Calendar calendar = TARGET();
  Date settlementDate = calendar.advance(today, Period(2, Days));
  Natural settlementDays = settlementDate - today;
  Settings::instance().evaluationDate() = today;

  // create overnight index for euro currency
  auto eoniaIndex = boost::make_shared<Eonia>();

  // create container for rate helpers
  std::vector<boost::shared_ptr<RateHelper>> rateHelpers;

  // create first 1d cash instrument for eonia curve - using deposit rate helper
  auto Q1D = boost::make_shared<SimpleQuote>(-0.0036);
  rateHelpers.push_back(boost::make_shared<DepositRateHelper>
   (Handle<Quote>(Q1D), Period(1, Days), eoniaIndex->fixingDays(),
   eoniaIndex->fixingCalendar(), eoniaIndex->businessDayConvention(), 
   eoniaIndex->endOfMonth(), eoniaIndex->dayCounter()));

  // create source data for eonia curve (period, rate)
  std::map<Period, Real> eoniaCurveData;
  eoniaCurveData.insert(std::make_pair(Period(1, Weeks), -0.00358));
  eoniaCurveData.insert(std::make_pair(Period(2, Weeks), -0.00357));
  eoniaCurveData.insert(std::make_pair(Period(3, Weeks), -0.00356));
  eoniaCurveData.insert(std::make_pair(Period(1, Months), -0.00358));
  eoniaCurveData.insert(std::make_pair(Period(2, Months), -0.00358));
  eoniaCurveData.insert(std::make_pair(Period(3, Months), -0.00358));
  eoniaCurveData.insert(std::make_pair(Period(4, Months), -0.00357));
  eoniaCurveData.insert(std::make_pair(Period(5, Months), -0.00356));
  eoniaCurveData.insert(std::make_pair(Period(6, Months), -0.00353));
  eoniaCurveData.insert(std::make_pair(Period(7, Months), -0.00351));
  eoniaCurveData.insert(std::make_pair(Period(8, Months), -0.00349));
  eoniaCurveData.insert(std::make_pair(Period(9, Months), -0.00344));
  eoniaCurveData.insert(std::make_pair(Period(10, Months), -0.0034));
  eoniaCurveData.insert(std::make_pair(Period(11, Months), -0.00336));
  eoniaCurveData.insert(std::make_pair(Period(1, Years), -0.00331));
  eoniaCurveData.insert(std::make_pair(Period(15, Months), -0.00314));
  eoniaCurveData.insert(std::make_pair(Period(18, Months), -0.00295));
  eoniaCurveData.insert(std::make_pair(Period(21, Months), -0.00273));
  eoniaCurveData.insert(std::make_pair(Period(2, Years), -0.00248));
  eoniaCurveData.insert(std::make_pair(Period(3, Years), -0.00138));
  eoniaCurveData.insert(std::make_pair(Period(4, Years), -0.0001245));
  eoniaCurveData.insert(std::make_pair(Period(5, Years), 0.0011945));
  eoniaCurveData.insert(std::make_pair(Period(6, Years), 0.00254));
  eoniaCurveData.insert(std::make_pair(Period(7, Years), 0.00387));
  eoniaCurveData.insert(std::make_pair(Period(8, Years), 0.0052));
  eoniaCurveData.insert(std::make_pair(Period(9, Years), 0.006474));
  eoniaCurveData.insert(std::make_pair(Period(10, Years), 0.007634));
  eoniaCurveData.insert(std::make_pair(Period(11, Years), 0.00868));
  eoniaCurveData.insert(std::make_pair(Period(12, Years), 0.0096045));
  eoniaCurveData.insert(std::make_pair(Period(15, Years), 0.0117245));
  eoniaCurveData.insert(std::make_pair(Period(17, Years), 0.0126797));
  eoniaCurveData.insert(std::make_pair(Period(20, Years), 0.0136245));
  eoniaCurveData.insert(std::make_pair(Period(25, Years), 0.01441));
  eoniaCurveData.insert(std::make_pair(Period(30, Years), 0.01479));

  // create other instruments for eonia curve - using ois rate helper
  std::for_each(eoniaCurveData.begin(), eoniaCurveData.end(),
   [settlementDays, &rateHelpers, &eoniaIndex](std::pair<Period, Real> p) -> void 
   { rateHelpers.push_back(boost::make_shared<OISRateHelper>(settlementDays,
   p.first, Handle<Quote>(boost::make_shared<SimpleQuote>(p.second)), eoniaIndex)); });

  // create piecewise term structure
  Handle<YieldTermStructure> eurOisCurve(boost::make_shared<PiecewiseYieldCurve<Discount, LogLinear>>
   (settlementDays, eoniaIndex->fixingCalendar(), rateHelpers, eoniaIndex->dayCounter()));

  // create equally-spaced (monthly) time grid
  Time maturity = 30.0;
  Size steps = 360;
  TimeGrid grid(maturity, steps);

  // print monthly discount factors
  std::for_each(grid.begin(), grid.end(), [&eurOisCurve](Time t) -> void
   { std::cout << eurOisCurve->discount(t) << std::endl; });

 }
 catch (std::exception& e) {
  std::cout << e.what() << std::endl;
 }
 return 0;
}

Saturday, June 4, 2016

Excel/VBA : Optimizing smooth OIS-adjusted Libor forward curve using Solver

Optimization for Libor forward curve has been presented in this blog post. This time, we will adjust the presented optimization procedure in such way, that OIS-adjusted Libor forward rates are going to be solved for a given fixed set of swap rates and OIS discount factors.

In a nutshell, justification for OIS-adjusted forward rates is the following :

  • In the "old world", we first bootstrapped Libor zero-coupon curve, from which we calculated Libor discount factors and Libor forward rates (for constructing floating leg coupons) at the same time. Only one curve was ever needed to accomplish this procedure. 
  • In the "new world", since all swap cash flows are discounted using OIS discount factors and par swap rates are still used for constructing swap fixed leg cash flows, forward rates (OIS-adjusted Libor forward rates) have to be slightly adjusted, in order to equate present value of all swap cash flows back to be zero.
All the relevant issues have been clearly explained in this research paper by Barclays.

Screenshots below are showing required Excel worksheet setups along with optimized OIS-adjusted Libor forward curve and required additions to existing VBA program needed to perform this optimization task. In order to validate the optimized OIS-adjusted Libor forward curve, a 10-year vanilla swap has been re-priced using optimized OIS-adjusted Libor forward rates and given set of fixed OIS discount factors.

When setting up VBA program, first implement the program presented in this blog post. After this, add two new functions (IRSwapOISPV, linearInterpolation) presented below into module XLSFunctions. Finally, remember to include references to Solver and Microsoft Scripting Runtime libraries.

Thanks for reading this blog.
-Mike





Public Function IRSwapOISPV(ByRef forwards As Excel.Range, ByRef OISDF As Excel.Range, ByVal swapRate As Double, _
ByVal yearsToMaturity As Integer, ByVal floatingLegPaymentFrequency As Integer, ByVal notional As Double) As Double
    '
    ' PV (fixed payer swap) = -swapRate * Q_fixed + Q_float
    ' where Q_fixed is sumproduct of fixed leg OIS discount factors and corresponding time accrual factors
    ' and Q_float is sumproduct of adjusted libor forward rates, OIS discount factors and corresponding time accrual factors
    ' assumption : fixed leg is always paid annually
    ' since fixed leg is paid annually and maturity is always integer, time accrual factor will always be exactly 1.0
    Dim f As Variant: f = forwards.Value2
    Dim DF As Variant: DF = OISDF.Value2
    Dim Q_fixed As Double
    Dim Q_float As Double
    Dim floatingLegTenor As Double: floatingLegTenor = 1 / CDbl(floatingLegPaymentFrequency)
    Dim nextFixedLegCouponDate As Integer: nextFixedLegCouponDate = 1
    Dim currentTimePoint As Double: currentTimePoint = CDbl(0)
    '
    Dim i As Integer
    For i = 1 To (yearsToMaturity * floatingLegPaymentFrequency)
        currentTimePoint = currentTimePoint + floatingLegTenor
        '
        ' update floating leg Q
        Q_float = Q_float + f(i, 1) * linearInterpolation(currentTimePoint, OISDF.Value2) * floatingLegTenor
        '
        ' update fixed leg Q, if current time point is coupon payment date
        If ((currentTimePoint - CDbl(nextFixedLegCouponDate)) = 0) Then
            Q_fixed = Q_fixed + linearInterpolation(currentTimePoint, OISDF.Value2)
            nextFixedLegCouponDate = nextFixedLegCouponDate + 1
        End If
    Next i
    IRSwapOISPV = (-swapRate * Q_fixed + Q_float) * notional
End Function
'
Public Function linearInterpolation(ByVal maturity As Double, ByRef curve As Variant) As Double
    '
    ' read range into Nx2 array
    Dim r As Variant: r = curve
    '
    Dim n As Integer: n = UBound(r, 1)
    '
    ' boundary checkings
    If ((r(LBound(r, 1), 1)) > maturity) Then linearInterpolation = r(LBound(r, 1), 2): Exit Function
    If ((r(UBound(r, 1), 1)) < maturity) Then linearInterpolation = r(UBound(r, 1), 2): Exit Function
    '
    Dim i As Long
    For i = 1 To n
        If ((r(i, 1) <= maturity) And (r(i + 1, 1) >= maturity)) Then
            '
            Dim y0 As Double: y0 = r(i, 2)
            Dim y1 As Double: y1 = r(i + 1, 2)
            Dim x0 As Double: x0 = r(i, 1)
            Dim x1 As Double: x1 = r(i + 1, 1)
            '
            linearInterpolation = y0 + (y1 - y0) * ((maturity - x0) / (x1 - x0))
            Exit For
        End If
    Next i
End Function
'

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