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
