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

Builder Pattern for C# Basket Default Swap Pricer

Keep your options open is a phrase, we often hear when making decisions concerning our lives. Needless to say, this great advice is directly applicable to any programming scheme. Now, I have to confess that many times (too many times) I feel a big hurry just to finish something quickly and that kind of an approach then leads to programs having no options for changes or modifications. Good example of this is BDS pricer program presented in the previous post.

PROBLEM

So, what is then wrong with that program? There should be nothing wrong with the program, but experienced programmers should catch a lot of different issues to be improved. Program is actually not very flexible as it is at the moment. However, I will concentrate on just one specific problem: the object creation part of the program. When looking at its client (class Program), we can see how it is polluted with object creation procedures (DiscountCurve, CreditCurve, NAGCopula) and hard-coded data. What if I would like to create my objects and data from text file or from Excel workbook? The mess (main code size explosion) would be even greater. Needless to say, the program as it is now, is not leaving much of those options open.

SOLUTION

Solution for this problem is Builder design pattern. Gang of Four (GOF) definition for Builder is the following.

"The intent of the Builder design pattern is to separate the construction of a complex object from its representation. By doing so the same construction process can create different representations."

In the new program design, BasketCDS class needs three different objects: Copula, CreditMatrix and DiscountCurve. Ultimately, client will feed BasketCDS class with these three objects wrapped inside Tuple. Before this, client is using concrete ExcelBuilder (implementation of abstract ComponentBuilder) for creating all three objects. So, client is effectively outsourcing the creation part of this complex object (a tuple, consisting of three objects) for ExcelBuilder. Just for example, instead of ExcelBuilder we could have TextFileBuilder or ConsoleBuilder. So, the source from which we are creating these objects (Copula, CreditMatrix and DiscountCurve) can be changed. Moreover, our client has no part in the object creation process. As far as I see it, this scheme is satisfying that GOF definition for Builder pattern. With this small change in design, we are now having some important options open.

PROGRAM FLOW

In order to understand the program data flow a bit better, a simple class UML has been presented in the picture below. I would like to stress the fact, that in no way this presentation is formally correct. However, the point is just to open up how this program is working on a very general level.

Client is requesting for creation of several objects from ExcelBuilder object, which is an implementation of abstract ComponentBuilder class. ExcelBuilder is then creating three central objects (CreditMatrix, DiscountCurve and Copula) into generic Tuple.

Client is then feeding BasketCDS with Tuple, which is consisting all required objects. BasketCDS object is starting Monte Carlo process and sending calculated leg values to Statistics object by using delegate methods. Finally, client is requesting the result (basket CDS spread) from Statistics object.

MathTools and MatrixTools classes are static utility classes, used by several objects for different mathematical or matrix operations. 

EXCEL BUILDER

For this example program, Excel is the most convenient platform for feeding parameters for our C# program. At this posting, I am only going to present the Excel interface and the results. For this kind of scheme, interfacing C# to Excel with Excel-DNA has been thoroughly covered in this blog posting.

My current Excel worksheet interface is the following.

The idea is simple. Client (Main) will delegate the creation part of three core objects (Copula, CreditMatrix and DiscountCurve) for ExcelBuilder class. ExcelBuilder class has been carefully configured to get specific data from specific ranges in Excel workbook. For example, it will read co-variance matrix from the named Excel range (_COVARIANCE). Personally, I am always using named Excel ranges instead of range addresses, in order to maintain flexibility when changing location of some range in the worksheet. By using this kind of a scheme, no further modifications will be needed in the program, since it communicates only with named Excel ranges.

After creating all those named ranges, a Button (Form Control) has been created into Excel worksheet, having program name Execute in its Assign Macro Box. I have been using Form Controls in order to get rid of any VBA code completely. Any public static void method in .NET code will be registered automatically by Excel-DNA as a macro in Excel. Thanks for this tip goes to Govert Van Drimmelen (inventor, developer and author of Excel-DNA).

PROGRAM RUN

After pressing Execute button, my C# program will write the following Basket CDS prices and program running time in seconds back to Excel worksheet (yellow areas).

THE PROGRAM

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

namespace ExcelBasketCDSPricer
{
    // CLIENT
    public static class Main
    {
        private static Stopwatch timer;
        private static dynamic Excel;
        private static ComponentBuilder builder;
        private static Tuple<Copula, CreditMatrix, DiscountCurve> components;
        private static List<BasketCDS> engines;
        private static List<Statistics> statistics;
        //
        public static void Execute()
        {
            try
            {
                // create timer and Excel objects
                timer = new Stopwatch();
                timer.Start();
                Excel = ExcelDnaUtil.Application;
                //
                // create Excel builder objects and request for required components for all pricers
                builder = new ExcelBuilder();
                components = builder.GetComponents();
                int nSimulations = (int)Excel.Range["_N"].Value2;
                int nReferences = (int)Excel.Range["_REFERENCES"].Value2;
                //
                // create basket pricers and statistics gatherers into containers
                engines = new List<BasketCDS>();
                statistics = new List<Statistics>();
                //
                // order updating for statistics gatherers from basket pricers
                for (int kth = 1; kth <= nReferences; kth++)
                {
                    statistics.Add(new Statistics(String.Concat(kth.ToString(), "-to-default")));
                    engines.Add(new BasketCDS(kth, nSimulations, nReferences, components));
                    engines[kth - 1].updateDefaultLeg += statistics[kth - 1].UpdateDefaultLeg;
                    engines[kth - 1].updatePremiumLeg += statistics[kth - 1].UpdatePremiumLeg;
                }
                //
                // run basket pricers
                engines.ForEach(it => it.Process());
                double[] prices =
                    statistics.Select(it => Math.Round(it.Spread() * 10000, 3)).ToArray<double>();
                timer.Stop();
                //
                // print statistics back to Excel
                Excel.Range["_PRICES"] = Excel.WorksheetFunction.Transpose(prices);
                Excel.Range["_RUNTIME"] = Math.Round(timer.Elapsed.TotalSeconds, 3);
            }
            catch (Exception e)
            {
                MessageBox.Show(e.Message);
            }
        }
    }
    //
    // *******************************************************************************
    // abstract base class for all component builders
    public abstract class ComponentBuilder
    {
        public abstract Copula CreateCopula();
        public abstract CreditMatrix CreateCreditMatrix();
        public abstract DiscountCurve CreateDiscountCurve();
        public abstract Tuple<Copula, CreditMatrix, DiscountCurve> GetComponents();
    }
    //
    // implementation for builder reading parameters from Excel
    public class ExcelBuilder : ComponentBuilder
    {
        dynamic Excel = null;
        private Copula copula = null;
        private CreditMatrix credit = null;
        private DiscountCurve curve = null;
        //
        public ExcelBuilder()
        {
            Excel = ExcelDnaUtil.Application;
        }
        public override Copula CreateCopula()
        {
            // create covariance matrix
            dynamic matrix = Excel.Range["_COVARIANCE"].Value2;
            int rows = matrix.GetUpperBound(0);
            int columns = matrix.GetUpperBound(1);
            double[,] covariance = new double[rows, columns];
            //
            for (int i = 0; i < rows; i++)
            {
                for (int j = 0; j < columns; j++)
                {
                    covariance[i, j] = (double)matrix.GetValue(i + 1, j + 1);
                }
            }
            //
            // create copula source information
            string source = (string)Excel.Range["_COPULA"].Value2;
            //
            // special parameters for NAG copula
            if (source == "NAG")
            {
                // generator id
                int genid = (int)Excel.Range["_GEN"].Value2;
                //
                // sub generator id
                int subgenid = (int)Excel.Range["_SUBGEN"].Value2;
                //
                // mode
                int mode = (int)Excel.Range["_MODE"].Value2;
                //
                // degrees of freedom
                int df = (int)Excel.Range["_DF"].Value2;
                //
                // create NAG copula
                copula = new NAGCopula(covariance, genid, subgenid, mode, df);
            }
            return copula;
        }
        public override CreditMatrix CreateCreditMatrix()
        {
            // create cds matrix
            dynamic matrix = Excel.Range["_CDS"].Value2;
            int rows = matrix.GetUpperBound(0);
            int columns = matrix.GetUpperBound(1);
            double[,] cds = new double[rows, columns];
            //
            for (int i = 0; i < rows; i++)
            {
                for (int j = 0; j < columns; j++)
                {
                    cds[i, j] = (double)matrix.GetValue(i + 1, j + 1);
                }
            }
            //
            // recovery estimate and discounting method as generic delegate
            double recovery = (double)Excel.Range["_RECOVERY"].Value2;
            Func<double, double> df = curve.GetDF;
            //
            CreditMatrix credit = new CreditMatrix(cds, df, recovery);
            return credit;
        }
        public override DiscountCurve CreateDiscountCurve()
        {
            // create zero-coupon curve
            dynamic matrix = Excel.Range["_ZERO"].Value2;
            int rows = matrix.GetUpperBound(0);
            int columns = matrix.GetUpperBound(1);
            double[,] zeroCouponCurve = new double[rows, columns];
            //
            for (int i = 0; i < rows; i++)
            {
                for (int j = 0; j < columns; j++)
                {
                    zeroCouponCurve[i, j] = (double)matrix.GetValue(i + 1, j + 1);
                }
            }
            //
            // create interpolation method
            int interpolationMethodID = (int)Excel.Range["_INTERPOLATION"].Value2;
            InterpolationAlgorithm interpolation = null;
            switch (interpolationMethodID)
            {
                case 1:
                    interpolation = MathTools.LinearInterpolation;
                    break;
                default:
                    throw new Exception("interpolation algorithm not defined");
            }
            //
            // create discounting method
            int discountingMethodID = (int)Excel.Range["_DISCOUNTING"].Value2;
            DiscountAlgorithm discounting = null;
            switch (discountingMethodID)
            {
                case 1:
                    discounting = MathTools.ContinuousDiscountFactor;
                    break;
                default:
                    throw new Exception("discounting algorithm not defined");
            }
            DiscountCurve curve = new DiscountCurve(zeroCouponCurve, interpolation, discounting);
            return curve;
        }
        public override Tuple<Copula, CreditMatrix, DiscountCurve> GetComponents()
        {
            copula = CreateCopula();
            curve = CreateDiscountCurve();
            credit = CreateCreditMatrix();
            Tuple<Copula, CreditMatrix, DiscountCurve> components =
                new Tuple<Copula, CreditMatrix, DiscountCurve>(copula, credit, curve);
            return components;
        }
    }
    //
    // *******************************************************************************
    public class BasketCDS
    {
        public PremiumLegUpdate updatePremiumLeg; // delegate for sending value
        public DefaultLegUpdate updateDefaultLeg; // delegate for sending value
        //
        private Copula copula; // NAG copula model
        private CreditMatrix cds; // credit matrix
        private DiscountCurve curve; // discount curve
        private int k; // kth-to-default
        private int m; // number of reference assets
        private int n; // number of simulations
        private int maturity; // basket cds maturity in years (integer)
        //
        public BasketCDS(int kth, int simulations, int maturity,
            Tuple<Copula, CreditMatrix, DiscountCurve> components)
        {
            this.k = kth;
            this.n = simulations;
            this.maturity = maturity;
            this.copula = components.Item1;
            this.cds = components.Item2;
            this.curve = components.Item3;
        }
        public void Process()
        {
            // request correlated random numbers from copula model
            double[,] randomArray = copula.Create(n, Math.Abs(Guid.NewGuid().GetHashCode()));
            m = randomArray.GetLength(1);
            //
            // process n sets of m correlated random numbers
            for (int i = 0; i < n; i++)
            {
                // create a set of m random numbers needed for one simulation round
                double[,] set = new double[1, m];
                for (int j = 0; j < m; j++)
                {
                    set[0, j] = randomArray[i, j];
                }
                //
                // calculate default times for reference name set
                calculateDefaultTimes(set);
            }
        }
        private void calculateDefaultTimes(double[,] arr)
        {
            // convert uniform random numbers into default times
            int cols = arr.GetLength(1);
            double[,] defaultTimes = new double[1, cols];
            //
            for (int j = 0; j < cols; j++)
            {
                // iteratively, find out the default tenor bucket
                double u = arr[0, j];
                double t = Math.Abs(Math.Log(1 - u));
                //
                for (int k = 0; k < cds.CumulativeHazardMatrix.GetLength(1); k++)
                {
                    int tenor = 0;
                    double dt = 0.0; double defaultTenor = 0.0;
                    if (cds.CumulativeHazardMatrix[k, j] >= t)
                    {
                        tenor = k;
                        if (tenor >= 1)
                        {
                            // calculate the exact default time for a given reference name
                            dt = -(1 / cds.HazardMatrix[k, j]) * Math.Log((1 - u)
                                / (Math.Exp(-cds.CumulativeHazardMatrix[k - 1, j])));
                            defaultTenor = tenor + dt;
                            //
                            // default time after basket maturity
                            if (defaultTenor >= maturity) defaultTenor = 0.0;
                        }
                        else
                        {
                            // hard-coded assumption
                            defaultTenor = 0.5;
                        }
                        defaultTimes[0, j] = defaultTenor;
                        break;
                    }
                }
            }
            // proceed to calculate leg values
            updateLegValues(defaultTimes);
        }
        private void updateLegValues(double[,] defaultTimes)
        {
            // check for defaulted reference names, calculate leg values 
            // and send statistics updates for BasketCDSStatistics
            int nDefaults = getNumberOfDefaults(defaultTimes);
            if (nDefaults > 0)
            {
                // for calculation purposes, remove zeros and sort matrix
                MatrixTools.RowMatrix_removeZeroValues(ref defaultTimes);
                MatrixTools.RowMatrix_sort(ref defaultTimes);
            }
            // calculate and send values for statistics gatherer
            double dl = calculateDefaultLeg(defaultTimes, nDefaults);
            double pl = calculatePremiumLeg(defaultTimes, nDefaults);
            updateDefaultLeg(dl);
            updatePremiumLeg(pl);
        }
        private int getNumberOfDefaults(double[,] arr)
        {
            int nDefaults = 0;
            for (int i = 0; i < arr.GetLength(1); i++)
            {
                if (arr[0, i] > 0.0) nDefaults++;
            }
            return nDefaults;
        }
        private double calculatePremiumLeg(double[,] defaultTimes, int nDefaults)
        {
            double dt = 0.0; double t; double p = 0.0; double v = 0.0;
            if ((nDefaults > 0) && (nDefaults >= k))
            {
                for (int i = 0; i < k; i++)
                {
                    if (i == 0)
                    {
                        // premium components from 0 to t1
                        dt = defaultTimes[0, i] - 0.0;
                        t = dt;
                        p = 1.0;
                    }
                    else
                    {
                        // premium components from t1 to t2, etc.
                        dt = defaultTimes[0, i] - defaultTimes[0, i - 1];
                        t = defaultTimes[0, i];
                        p = (m - i) / (double)m;
                    }
                    v += (curve.GetDF(t) * dt * p);
                }
            }
            else
            {
                for (int i = 0; i < maturity; i++)
                {
                    v += curve.GetDF(i + 1);
                }
            }
            return v;
        }
        private double calculateDefaultLeg(double[,] defaultTimes, int nDefaults)
        {
            double v = 0.0;
            if ((nDefaults > 0) && (nDefaults >= k))
            {
                v = (1 - cds.Recovery) * curve.GetDF(defaultTimes[0, k - 1]) * (1 / (double)m);
            }
            return v;
        }
    }
    //
    // *******************************************************************************
    // abstract base class for all copula models
    public abstract class Copula
    {
        // request matrix of correlated uniform random numbers 
        // number of rows are given argument n
        // Number of columns are inferred from the size of covariance matrix and
        public abstract double[,] Create(int n, int seed);
    }
    //
    // NAG G05 COPULAS WRAPPER
    public class NAGCopula : Copula
    {
        private double[,] covariance;
        private int genID;
        private int subGenID;
        private int mode;
        private int df;
        private int m;
        private int errorNumber;
        private double[] r;
        private double[,] result;
        //
        public NAGCopula(double[,] covariance, int genID,
            int subGenID, int mode, int df = 0)
        {
            // ctor : create correlated uniform random numbers
            // degrees-of-freedom parameter (df), being greater than zero, 
            // will automatically trigger the use of student copula
            this.covariance = covariance;
            this.genID = genID;
            this.subGenID = subGenID;
            this.mode = mode;
            this.df = df;
            this.m = covariance.GetLength(1);
            r = new double[m * (m + 1) + (df > 0 ? 2 : 1)];
        }
        public override double[,] Create(int n, int seed)
        {
            result = new double[n, m];
            G05.G05State g05State = new G05.G05State(genID, subGenID, new int[1] { seed }, out errorNumber);
            if (errorNumber != 0) throw new Exception("G05 state failure");
            //
            if (this.df != 0)
            {
                // student copula
                G05.g05rc(mode, n, df, m, covariance, r, g05State, result, out errorNumber);
            }
            else
            {
                // gaussian copula
                G05.g05rd(mode, n, m, covariance, r, g05State, result, out errorNumber);
            }
            //
            if (errorNumber != 0) throw new Exception("G05 method failure");
            return result;
        }
    }
    //
    // *******************************************************************************
    public class CreditMatrix
    {
        private Func<double, double> df;
        private double recovery;
        private double[,] CDSSpreads;
        private double[,] survivalMatrix;
        private double[,] hazardMatrix;
        private double[,] cumulativeHazardMatrix;
        //
        public CreditMatrix(double[,] CDSSpreads,
            Func<double, double> discountFactor, double recovery)
        {
            this.df = discountFactor;
            this.CDSSpreads = CDSSpreads;
            this.recovery = recovery;
            createSurvivalMatrix();
            createHazardMatrices();
        }
        // public read-only accessors to class data
        public double[,] HazardMatrix { get { return this.hazardMatrix; } }
        public double[,] CumulativeHazardMatrix { get { return this.cumulativeHazardMatrix; } }
        public double Recovery { get { return this.recovery; } }
        //
        private void createSurvivalMatrix()
        {
            // bootstrap matrix of survival probabilities from given CDS data
            int rows = CDSSpreads.GetUpperBound(0) + 2;
            int cols = CDSSpreads.GetUpperBound(1) + 1;
            survivalMatrix = new double[rows, cols];
            //
            double term = 0.0; double firstTerm = 0.0; double lastTerm = 0.0;
            double terms = 0.0; double quotient = 0.0;
            int i = 0; int j = 0; int k = 0;
            //
            for (i = 0; i < rows; i++)
            {
                for (j = 0; j < cols; j++)
                {
                    if (i == 0) survivalMatrix[i, j] = 1.0;
                    if (i == 1) survivalMatrix[i, j] = (1 - recovery) / ((1 - recovery) + 1 * CDSSpreads[i - 1, j] / 10000);
                    if (i > 1)
                    {
                        terms = 0.0;
                        for (k = 0; k < (i - 1); k++)
                        {
                            term = df(k + 1) * ((1 - recovery) * survivalMatrix[k, j] -
                            (1 - recovery + 1 * CDSSpreads[i - 1, j] / 10000) * survivalMatrix[k + 1, j]);
                            terms += term;
                        }
                        quotient = (df(i) * ((1 - recovery) + 1 * CDSSpreads[i - 1, j] / 10000));
                        firstTerm = (terms / quotient);
                        lastTerm = survivalMatrix[i - 1, j] * (1 - recovery) / (1 - recovery + 1 * CDSSpreads[i - 1, j] / 10000);
                        survivalMatrix[i, j] = firstTerm + lastTerm;
                    }
                }
            }
        }
        private void createHazardMatrices()
        {
            // convert matrix of survival probabilities into two hazard rate matrices
            int rows = survivalMatrix.GetUpperBound(0);
            int cols = survivalMatrix.GetUpperBound(1) + 1;
            hazardMatrix = new double[rows, cols];
            cumulativeHazardMatrix = new double[rows, cols];
            int i = 0; int j = 0;
            //
            for (i = 0; i < rows; i++)
            {
                for (j = 0; j < cols; j++)
                {
                    cumulativeHazardMatrix[i, j] = -Math.Log(survivalMatrix[i + 1, j]);
                    if (i == 0) hazardMatrix[i, j] = cumulativeHazardMatrix[i, j];
                    if (i > 0) hazardMatrix[i, j] = (cumulativeHazardMatrix[i, j] - cumulativeHazardMatrix[i - 1, j]);
                }
            }
        }
    }
    //
    // *******************************************************************************
    // delegate methods for interpolation and discounting
    public delegate double InterpolationAlgorithm(double t, ref double[,] curve);
    public delegate double DiscountAlgorithm(double t, double r);
    //
    public class DiscountCurve
    {
        // specific algorithms for interpolation and discounting
        private InterpolationAlgorithm interpolationMethod;
        private DiscountAlgorithm discountMethod;
        private double[,] curve;
        //
        public DiscountCurve(double[,] curve,
            InterpolationAlgorithm interpolationMethod, DiscountAlgorithm discountMethod)
        {
            this.curve = curve;
            this.interpolationMethod = interpolationMethod;
            this.discountMethod = discountMethod;
        }
        public double GetDF(double t)
        {
            // get discount factor from discount curve
            return discountMethod(t, interpolation(t));
        }
        private double interpolation(double t)
        {
            // get interpolation from discount curve
            return interpolationMethod(t, ref this.curve);
        }
    }
    //
    // *******************************************************************************
    // collection of methods for different types of mathematical operations
    public static class MathTools
    {
        public static double LinearInterpolation(double t, ref double[,] curve)
        {
            int n = curve.GetUpperBound(0) + 1;
            double v = 0.0;
            //
            // boundary checkings
            if ((t < curve[0, 0]) || (t > curve[n - 1, 0]))
            {
                if (t < curve[0, 0]) v = curve[0, 1];
                if (t > curve[n - 1, 0]) v = curve[n - 1, 1];
            }
            else
            {
                // iteration through all given curve points
                for (int i = 0; i < n; i++)
                {
                    if ((t >= curve[i, 0]) && (t <= curve[i + 1, 0]))
                    {
                        v = curve[i, 1] + (curve[i + 1, 1] - curve[i, 1]) * (t - (i + 1));
                        break;
                    }
                }
            }
            return v;
        }
        public static double ContinuousDiscountFactor(double t, double r)
        {
            return Math.Exp(-r * t);
        }
    }
    //
    // *******************************************************************************
    // collection of methods for different types of matrix operations
    public static class MatrixTools
    {
        public static double[,] CorrelationToCovariance(double[,] corr, double[] stdev)
        {
            // transform correlation matrix to co-variance matrix
            double[,] cov = new double[corr.GetLength(0), corr.GetLength(1)];
            //
            for (int i = 0; i < cov.GetLength(0); i++)
            {
                for (int j = 0; j < cov.GetLength(1); j++)
                {
                    cov[i, j] = corr[i, j] * stdev[i] * stdev[j];
                }
            }
            //
            return cov;
        }
        public static void RowMatrix_sort(ref double[,] arr)
        {
            // sorting a given row matrix to ascending order
            // input must be 1 x M matrix
            // bubblesort algorithm implementation
            int cols = arr.GetUpperBound(1) + 1;
            double x = 0.0;
            //
            for (int i = 0; i < (cols - 1); i++)
            {
                for (int j = (i + 1); j < cols; j++)
                {
                    if (arr[0, i] > arr[0, j])
                    {
                        x = arr[0, i];
                        arr[0, i] = arr[0, j];
                        arr[0, j] = x;
                    }
                }
            }
        }
        public static void RowMatrix_removeZeroValues(ref double[,] arr)
        {
            // removes zero values from a given row matrix
            // input must be 1 x M matrix
            List<double> temp = new List<double>();
            int cols = arr.GetLength(1);
            int counter = 0;
            for (int i = 0; i < cols; i++)
            {
                if (arr[0, i] > 0)
                {
                    counter++;
                    temp.Add(arr[0, i]);
                }
            }
            if (counter > 0)
            {
                arr = new double[1, temp.Count];
                for (int i = 0; i < temp.Count; i++)
                {
                    arr[0, i] = temp[i];
                }
            }
            else
            {
                arr = null;
            }
        }
    }
    //
    // *******************************************************************************
    public delegate void PremiumLegUpdate(double v);
    public delegate void DefaultLegUpdate(double v);
    //
    public class Statistics
    {
        // data structures for storing leg values
        private List<double> premiumLeg;
        private List<double> defaultLeg;
        private string ID;
        //
        public Statistics(string ID)
        {
            this.ID = ID;
            premiumLeg = new List<double>();
            defaultLeg = new List<double>();
        }
        public void UpdatePremiumLeg(double v)
        {
            premiumLeg.Add(v);
        }
        public void UpdateDefaultLeg(double v)
        {
            defaultLeg.Add(v);
        }
        public void PrettyPrint()
        {
            // hard-coded 'report' output
            Console.WriteLine("{0} : {1} bps", ID, Math.Round(Spread() * 10000, 2));
        }
        public double Spread()
        {
            return defaultLeg.Average() / premiumLeg.Average();
        }
    }
}

CONCLUSIONS

In this blog posting, Builder design pattern was applied to solve some of the problems observed in hard-coded program. By outsourcing the object creation for Builder (ExcelBuilder), client code size explosion has been avoided completely. Moreover, the source for objects creation is now open for new implementations (ConsoleBuilder, TextFileBuilder, etc). The program is still not as flexible as it could be, but that is another story.

Again, Thanks for Govert Van Drimmelen for his amazing Excel-DNA. 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. Excel-DNA is an open-source project, and we (the happy users) can invest its future development by making a donation.

Last week, I spent three days in Datasim training course. During the course, we went systematically through most of the GOF design patterns in C#, using traditional object-oriented approach. Plus, the instructor was also presenting, how to use C# generics (delegates) for implementing some of the design patterns, using functional programming approach. The course was truly having a great balance between theory and programming. So, if anyone is looking for very practical hands-on training for this GOF stuff, this course is getting all five stars from me.

Thank You for reading my blog. I wish you have discovered some usage for Builder pattern in your programs. Mike.