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

Configurable C# Monte Carlo zero-coupon bond pricer in Excel

Monte Carlo design, what was presented in my previous blog article could actually be used, not only for pricing options, but for all applications where one would like to simulate stochastic process. One such useful application would be to simulate short rate process for pricing zero-coupon bonds, for example.

In this article, we investigate how to implement such pricing scheme with the current design for one-factor short-rate models. As a result of this small project, we will have a framework for pricing zero-coupon bonds with any given one-factor short-rate model. Specifically, we implement zero-coupon bond pricing schemes for Vasicek model and CIR model. Since both of these models are nicely tractable Gaussian models, we also compare the simulated bond prices with the exact analytical model prices. Before going further with this article, make sure that you have clear understanding of the design presented in my previous blog article.

PROJECT OUTCOME

The outcome of this small project is fully configurable C# Monte Carlo pricer application in Excel. Application can be used to price zero-coupon bond prices with any different types of one-factor short-rate models. The application gets all the required input parameters directly from Excel, then performs calculations in C# and finally returns calculation results back to Excel. Excel and C# are interfaced with Excel-DNA and Excel itself is used only as data input/output platform, exactly like presented in this blog article.

PREPARATORY TASKS

Download and unzip Excel-DNA zip file to be ready when needed. There is also a step-by-step word documentation file available within the distribution folder. In this project, we are going to follow these instructions.

EXTENSIONS FOR CURRENT DESIGN

As a base design for this application, we will use the design presented in my previous blog article. There will be no changes made to core components of this base design (SDE, Discretization, RandomGenerator). We are going to implement the following extensions for this base design.

• Two new concrete class implementations for abstract class SDE (Vasicek, CIR)
• One new concrete class implementation for abstract class Pricer (BondPricer)

There will be some small modifications made to Builder component.

C# PROGRAM

Create a new C# Class Library project (BondPricer), save the project and copyPaste the following code blocks into separate cs files. These classes are all we need from our current base design. For the sake of clarity, I made a decision not to include any option-related classes to this design.

interface ISDE

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace BondPricer
{
    public interface ISDE
    { 
        // methods for calculating drift and diffusion term of stochastic differential equation
        double drift(double s, double t);
        double diffusion(double s, double t);
    }
}

abstract class SDE

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace BondPricer
{
    public abstract class SDE :ISDE
    { 
        // abstract class implementing ISDE interface
        public abstract double drift(double s, double t);
        public abstract double diffusion(double s, double t);
    }
}

interface IDiscretization

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace BondPricer
{
    public interface IDiscretization
    {
        // method for discretizing stochastic differential equation
        double next(double s, double t, double dt, double rnd);
    }
}

abstract class Discretization

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace BondPricer
{
    public abstract class Discretization : IDiscretization
    {
        // abstract class implementing IDiscretization interface
        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);
    }
}

concrete class EulerDiscretization

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace BondPricer
{
    public class EulerDiscretization : Discretization
    { 
        // concrete implementation for Euler discretization scheme
        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;
        }
    }
}

interface IRandomGenerator

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace BondPricer
{
    public interface IRandomGenerator
    {
        // method for generating normally distributed random variable
        double getRandom();
    }
}

abstract class RandomGenerator

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace BondPricer
{
    public abstract class RandomGenerator : IRandomGenerator
    {
        // abstract class implementing IRandomGenerator interface
        public abstract double getRandom();
    }
}

concrete class NormalApproximation

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace BondPricer
{
    public class NormalApproximation : RandomGenerator
    {
        // concrete implementation for normal random variable approximation
        // normRand = sum of 12 independent uniformly disctibuted random numbers, minus 6
        private Random random;
        public NormalApproximation()
        {
            random = new Random();
        }
        public override double getRandom()
        {
            // implementation uses C# uniform random generator
            double[] rnd = new double[12];
            Func<double> generator = () => { return random.NextDouble(); };
            return rnd.Select(r => generator()).Sum() - 6.0;
        }
    }
}

concrete class MonteCarloEngine

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace BondPricer
{
    public delegate void PathSender(ref double[] path);
    public delegate void ProcessStopper();
    //
    public class MonteCarloEngine
    {
        private SDE sde;
        private Discretization discretization;
        private RandomGenerator randomGenerator;
        private long paths;
        private int steps;
        public event PathSender sendPath;
        public event ProcessStopper stopProcess;
        //
        public MonteCarloEngine(Tuple<SDE, Discretization, RandomGenerator> parts, long paths, int steps)
        {
            sde = parts.Item1;
            discretization = parts.Item2;
            randomGenerator = parts.Item3;
            this.paths = paths;
            this.steps = steps;
        }
        public MonteCarloEngine(Builder builder, long paths, int steps)
        {
            Tuple<SDE, Discretization, RandomGenerator> parts = builder.build();
            sde = parts.Item1;
            discretization = parts.Item2;
            randomGenerator = parts.Item3;
            this.paths = paths;
            this.steps = steps;
        }
        public void run()
        {
            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++)
            {
                path[0] = vOld = discretization.InitialPrice;
                //
                for (int j = 1; j <= steps; j++)
                {
                    // get next value using discretization scheme
                    vNew = discretization.next(vOld, (dt * j), dt, randomGenerator.getRandom());
                    path[j] = vNew; vOld = vNew;
                }
                sendPath(ref path); // send one simulated path to pricer to be processed
            }
            stopProcess(); // simulation ends - notify pricer
        }
    }
}

interface IBuilder

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace BondPricer
{
    public interface IBuilder
    {
        // method for creating all the needed objects for asset price simulations
        Tuple<SDE, Discretization, RandomGenerator> build();
    }
}

abstract class Builder

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace BondPricer
{
    public abstract class Builder : IBuilder
    {
        // abstract class implementing IBuilder interface
        public abstract Tuple<SDE, Discretization, RandomGenerator> build();
    }
}

interface IPricer

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace BondPricer
{
    public interface IPricer
    {
        void processPath(ref double[] path);
        void calculate();
        double price();
    }
}

abstract class Pricer

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace BondPricer
{
    public delegate double OneFactorPayoff(double spot, double strike);
    //
    public abstract class Pricer : IPricer
    {
        protected OneFactorPayoff payoff; // delegate function for payoff calculation
        protected Func<double> discountFactor; // generic delegate function for discount factor
        protected double v; // option price
        protected long paths; // running counter
        //
        public Pricer() { }
        public Pricer(OneFactorPayoff payoff, Func<double> discountFactor)
        {
            this.payoff = payoff; this.discountFactor = discountFactor;
        }
        public abstract void processPath(ref double[] path);
        public abstract void calculate();
        public abstract double price();
    }
}

After implementing all the interfaces and classes described above, add the following new implementations to this design.

concrete class ExcelBuilder

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

namespace BondPricer
{
    public class ExcelBuilder : Builder
    {
        private double t;
        private dynamic Excel = ExcelDnaUtil.Application;
        //
        public ExcelBuilder(double t) { this.t = t; }
        public override Tuple<SDE, Discretization, RandomGenerator> build()
        {
            // build all objects needed for asset path simulations
            SDE sde = build_SDE();
            Discretization discretization = build_discretization(sde);
            RandomGenerator randomGenerator = build_randomGenerator();
            return new Tuple<SDE, Discretization, RandomGenerator>(sde, discretization, randomGenerator);
        }
        private SDE build_SDE()
        {
            SDE sde = null;
            string sdeType = (string)Excel.Range("_shortRateModel").Value;
            double longTermRate = (double)Excel.Range("_longTermRate").Value2;
            double reversionSpeed = (double)Excel.Range("_reversionSpeed").Value2;
            double rateVolatility = (double)Excel.Range("_rateVolatility").Value2;
            //
            if (sdeType == "VASICEK") sde = new Vasicek(longTermRate, reversionSpeed, rateVolatility);
            if (sdeType == "CIR") sde = new Vasicek(longTermRate, reversionSpeed, rateVolatility);
            // insert new stochastic model choices here
            return sde;
        }
        private Discretization build_discretization(SDE sde)
        {
            Discretization discretization = null;
            string discretizationType = (string)Excel.Range("_discretization").Value;
            //
            if (discretizationType == "EULER")
            {
                double initialPrice = (double)Excel.Range("_spotRate").Value2;
                discretization = new EulerDiscretization(sde, initialPrice, t);
            }
            // insert new discretization scheme choices here
            return discretization;
        }
        private RandomGenerator build_randomGenerator()
        {
            RandomGenerator randomGenerator = null;
            string randomGeneratorType = (string)Excel.Range("_randomGenerator").Value;
            //
            if (randomGeneratorType == "CLT")
            {
                randomGenerator = new NormalApproximation();
            }
            // insert new random generator choices here
            return randomGenerator;
        }
    }
}

concrete class Vasicek

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace BondPricer
{
    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;
        }
    }
}

concrete class CIR

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace BondPricer
{
    public class CIR : SDE
    {
        private double longTermRate;
        private double reversionSpeed;
        private double rateVolatility;
        //
        public CIR(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 * Math.Sqrt(s);
        }
    }
}

concrete class ZCBPricer

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace BondPricer
{
    public class ZCBPricer : Pricer
    {
        private double t;
        private double dt;
        //
        public ZCBPricer(double t, double dt) : base() { this.t = t; this.dt = dt; }
        //
        public override void processPath(ref double[] path)
        {
            int n = (int)(t / dt);
            double integral = 0.0;
            for (int i = 0; i < n; i++) integral += path[i] * dt;
            v += integral;
            paths++;
        }
        public override void calculate()
        {
            v /= paths; // calculate expectation
        }
        public override double price()
        {
            return Math.Exp(-v); // return discount factor
        }
    }
}

concerete class BondPricer (VBA program will call run method of this class).

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

namespace BondPricer
{
    public static class BondPricer
    {
        private static dynamic Excel;
        private static MonteCarloEngine engine;
        private static List<Pricer> bonds;
        //
        public static void run()
        {
            try
            {
                // create Excel application
                Excel = ExcelDnaUtil.Application;
                //
                // fetch bond pricing parameters from named Excel ranges
                int steps = (int)Excel.Range("_steps").Value2;
                long paths = (long)Excel.Range("_paths").Value2;
                //
                // read maturities from named range
                dynamic maturityArray = Excel.Range["_maturity"].Value2;
                int rows = maturityArray.GetUpperBound(0);
                double maxMaturity = (double)maturityArray.GetValue(rows, 1);
                double dt = (double)maxMaturity / steps;
                //
                // create Monte Carlo engine
                engine = new MonteCarloEngine(new ExcelBuilder(maxMaturity), paths, steps);
                bonds = new List<Pricer>();
                //
                // add pricers to list, order path updates and stop notification from engine
                for (int i = 0; i < rows; i++)
                {
                    double t = (double)maturityArray.GetValue(i + 1, 1);
                    bonds.Add(new ZCBPricer(t, dt));
                    engine.sendPath += bonds[i].processPath;
                    engine.stopProcess += bonds[i].calculate;
                }
                //
                // run Monte Carlo engine
                engine.run();
                //
                // create output array and add zero-coupon prices
                double[] zeroCouponCurve = new double[rows];
                for (int i = 0; i < rows; i++) zeroCouponCurve[i] = bonds[i].price();
                //
                // print curve back to Excel
                Excel.Range["_zcb"] = Excel.WorksheetFunction.Transpose(zeroCouponCurve);
            }
            catch (Exception e)
            {
                MessageBox.Show(e.Message.ToString());
            }
        }
    }
}

One essential input for this application is the array of zero-coupon bond maturities, fetched directly from Excel worksheet. Time in years for the longest bond maturity is saved into local variable maxMaturity and time step length (dt) is calculated as maxMaturity divided by the number of steps. In the Excel example below, maxMaturity is 10 years, steps per simulation is 2500 and dt is 0.004. During the simulation process, MonteCarloEngine will then always simulate a path having maturity equal to maxMaturity (10.0) and having time step length equal to dt (0.004).

Now, for each simulated short-rate path, generated by MonteCarloEngine, ZCBPricer will get a path update (array of doubles) by its processPath method. ZCBPricer will then calculate zero-coupon bond price as integral of this short rate path. However, time length of short-rate path in years, provided by MonteCarloEngine (10.0), could be longer than the maturity of zero-coupon bond defined in ZCBPricer (double t). In order to process only the rates where time is less or equal to maturity of zero-coupon bond, local variable n (t / dt) is used in loop condition of pricer's processPath method. So, the pricer will always receive the whole path, but it picks only the items it needs.

In essence, most of the time our MonteCarloEngine is "working more than it should", since most of the bond maturities are less than maxMaturity and MonteCarloEngine will always simulate a path having maturity equal to maxMaturity. However, this design has been deliberate choice for this particular application, in order to maintain a scheme with only one engine plus n pricers. Needless to say, the price we pay for this particular scheme is increased processing time.

EXCEL-DNA INTEGRATION

After implementing all the previous cs files into C# Class Library project, we are receiving a lot of errors. However, all errors should be related to missing references to Excel-DNA integration library and Windows Forms library. Next, carefully follow the instructions described here in step two.

In a nutshell

• add reference to Excel-DNA library (ExcelDna.Integration.dll). From the properties of this reference, set Copy Local to be False.
• add reference to Windows Forms library (System.Windows.Forms)
• create BondPricer.dna file (consisting XML tags). DnaLibrary Name="BondPricer" and Path="BondPricer.dll". From the properties of this dna file, set Copy to Output Directory to be Copy if newer.
• copy ExcelDna.xll file into your project folder and rename it to be BondPricer.xll. From the properties of this xll file, set Copy to Output Directory to be Copy if newer.
  

Make sure, that all properties for these references and files are exactly the same as described in this post. After adding all the required references, files and building this program once again, my project folder has the following four files.

USER INTERFACE AND C#-VBA INTEGRATION

The essence of this part of the process has been described here in step three. Open a new Excel workbook. Create the following source data into worksheet. In the picture below, I have also included Excel formulas for calculating bond yield from zero-coupon bond price and calculating exact analytical price for Vasicek model (coefficients A, B).

From Excel Name Manager (Formulas - Name Manager), set the following range names.

In VB editor, create the following event handling program for CommandButton (run).

TEST RUN

At this point, our zero-coupon bond pricing application is ready for test run. While the current Excel workbook is open, doubleClick BondPricer.xll file in your \\BondPricer\bin\Release folder. After this, xll file content can be used by Excel and our C# program is available to be called from VBA program (Application.Run). VBA program will call and start C# program run, which then reads all input data, performs calculations and sends result data back to Excel worksheet.

After pressing command button in Excel interface, my C# MC zero-coupon bond pricer application simulated the following prices for all requested maturities from 1 year to 10 years (Monte Carlo simulation).

AFTERTHOUGHTS

This small project was presenting one possible Monte Carlo design for pricing zero-coupon bonds with one-factor short-rate processes. For this application, we basically just extended our existing flexible Monte Carlo design. Comparing our simulated bond prices with the exact analytical counterparties, we can conclude that our application is pricing bonds correctly, related to short-rate model used.

Thanks for reading and have a nice summer!
-Mike Juniperhill