More Copula stuff coming again. This time, I wanted to open my implementation for Basket Default Swap (BDS) pricer. It is assumed, that the reader is familiar with this exotic credit product.
PROGRAM DESIGN
The core of the program is
BasketCDS class, which is performing the actual Monte Carlo process. This class is using
NAGCopula class for creating correlated uniform random numbers needed for all simulations. It also uses CreditCurve class for retrieving credit-related data for calculating default times. Finally, it holds a delegate method for desired discounting method. For storing the calculated default and premium leg values, delegate methods are used to send these values to
Statistics class, which then stores the both leg values into its own data structures. In essence, BasketCDS class is doing the calculation and Statistics class is responsible for storing the calculated values and reporting its own set of desired (and configured) statistics when requested.
As already mentioned,
CreditCurve class is holding all the credit-related information (CDS curve, hazard rates, recovery rate).
DiscountCurve class holds the information on zero-coupon curve. In this class, discounting and interpolation methods are not hard-coded inside the class. Instead, the class is holding two delegate methods for this purpose. The actual implementations for the both methods are stored in static
MathTools class, which is nothing more than just collection of methods for different types of mathematical operations. Some matrix operations are needed for handling matrices when calculating default times in BasketCDS class. For this purpose, static
MatrixTools class is providing collection of methods for different types of matrix operations. Needless to say, new method implementations can be provided for concrete methods used by these delegates.
Source market data (co-variance matrix, CDS curve and zero-coupon curve) for this example program has been completely hard-coded. However, with the examples presented in
here or in
here, the user should be able to interface this program back to Excel, if so desired.
MONTE CARLO PROCESS
The client example program is calculating five fair spreads for all kth-to-default basket CDS for all five reference names in the basket. In a nutshell, Monte Carlo procedure to calculate kth-to-default BDS spread is the following.
- Simulate five correlated uniformly distributed random numbers by using Gaussian or Student Copula model.
- Transform simulated five uniform random numbers to default times by using inverse CDF of exponential distribution.
- By using term structure of hazard rates, define the exact default times for all five names.
- Calculate the value for default leg and premium leg.
- Store the both calculated leg values.
- Repeat steps from one to five N times.
- Finally, calculate BDS spread by dividing average of all stored default leg values with average of all stored premium leg values.
In order to run this program, one needs to create a new C# console project and copyPaste the program given below into a new cs file. It is also assumed, that the user is a registrated
NAG licence holder. Other users have to provide their own implementation for Copula needed in this program or use
my implementation. With the second approach, a bit more time will be spent in adjusting this program for all the changes required.
Results for one test run for the both Copula models is presented in the following console screenshot.
Observed results are behaving as theory is suggesting. BDS prices are decreasing as more reference names are needed for triggering the actual default. Also, tail-effect created by the use of Student Copula model is captured as higher prices after 1st-to-default BDS. Computational help provided by NAG G05 library is just amazing. After you outsource the handling of Copula model for external library, there is not so much to program either.
That is all I wanted to share with you this time. Thanks for reading my blog. Mike.
THE PROGRAM
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using NagLibrary;
//
namespace BasketCDSPricer
{
// delegates for updating leg values
public delegate void PremiumLegUpdate(double v);
public delegate void DefaultLegUpdate(double v);
//
// delegate methods for interpolation and discounting
public delegate double InterpolationAlgorithm(double t, ref double[,] curve);
public delegate double DiscountAlgorithm(double t, double r);
//
//
// *************************************************************
// CLIENT PROGRAM
class Program
{
static void Main(string[] args)
{
try
{
// create covariance matrix, discount curve, credit curve and NAG gaussian copula model
double[,] covariance = createCovarianceMatrix();
DiscountCurve zero = new DiscountCurve(createDiscountCurve(),
MathTools.LinearInterpolation, MathTools.ContinuousDiscountFactor);
CreditCurve cds = new CreditCurve(createCDSCurve(), zero.GetDF, 0.4);
NAGCopula copula = new NAGCopula(covariance, 1, 1, 2);
//
// create containers for basket pricers and statistics
// set number of simulations and reference assets
List<BasketCDS> engines = new List<BasketCDS>();
List<Statistics> statistics = new List<Statistics>();
int nSimulations = 250000;
int nReferences = 5;
//
// create basket pricers and statistics gatherers into containers
// 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, copula, cds, zero.GetDF));
engines[kth - 1].updateDefaultLeg += statistics[kth - 1].UpdateDefaultLeg;
engines[kth - 1].updatePremiumLeg += statistics[kth - 1].UpdatePremiumLeg;
}
// run basket pricers and print statistics
engines.ForEach(it => it.Process());
statistics.ForEach(it => it.PrettyPrint());
}
catch (Exception e)
{
Console.WriteLine(e.Message);
}
}
//
// hard-coded co-variance matrix data
static double[,] createCovarianceMatrix()
{
double[,] covariance = new double[5, 5];
covariance[0, 0] = 0.001370816; covariance[1, 0] = 0.001113856; covariance[2, 0] = 0.001003616;
covariance[3, 0] = 0.000937934; covariance[4, 0] = 0.00040375;
covariance[0, 1] = 0.001113856; covariance[1, 1] = 0.001615309; covariance[2, 1] = 0.000914043;
covariance[3, 1] = 0.000888594; covariance[4, 1] = 0.000486823;
covariance[0, 2] = 0.001003616; covariance[1, 2] = 0.000914043; covariance[2, 2] = 0.001302899;
covariance[3, 2] = 0.000845642; covariance[4, 2] = 0.00040185;
covariance[0, 3] = 0.000937934; covariance[1, 3] = 0.000888594; covariance[2, 3] = 0.000845642;
covariance[3, 3] = 0.000954866; covariance[4, 3] = 0.000325403;
covariance[0, 4] = 0.00040375; covariance[1, 4] = 0.000486823; covariance[2, 4] = 0.00040185;
covariance[3, 4] = 0.000325403; covariance[4, 4] = 0.001352532;
return covariance;
}
//
// hard-coded discount curve data
static double[,] createDiscountCurve()
{
double[,] curve = new double[5, 2];
curve[0, 0] = 1.0; curve[0, 1] = 0.00285;
curve[1, 0] = 2.0; curve[1, 1] = 0.00425;
curve[2, 0] = 3.0; curve[2, 1] = 0.00761;
curve[3, 0] = 4.0; curve[3, 1] = 0.0119;
curve[4, 0] = 5.0; curve[4, 1] = 0.01614;
return curve;
}
//
// hard-coded cds curve data
static double[,] createCDSCurve()
{
double[,] cds = new double[5, 5];
cds[0, 0] = 17.679; cds[0, 1] = 19.784;
cds[0, 2] = 11.242; cds[0, 3] = 23.5;
cds[0, 4] = 39.09; cds[1, 0] = 44.596;
cds[1, 1] = 44.921; cds[1, 2] = 27.737;
cds[1, 3] = 50.626; cds[1, 4] = 93.228;
cds[2, 0] = 54.804; cds[2, 1] = 64.873;
cds[2, 2] = 36.909; cds[2, 3] = 67.129;
cds[2, 4] = 107.847; cds[3, 0] = 83.526;
cds[3, 1] = 96.603; cds[3, 2] = 57.053;
cds[3, 3] = 104.853; cds[3, 4] = 164.815;
cds[4, 0] = 96.15; cds[4, 1] = 115.662;
cds[4, 2] = 67.82; cds[4, 3] = 118.643;
cds[4, 4] = 159.322;
return cds;
}
}
//
//
// *************************************************************
public class BasketCDS
{
public PremiumLegUpdate updatePremiumLeg; // delegate for sending value
public DefaultLegUpdate updateDefaultLeg; // delegate for sending value
private NAGCopula copula; // NAG copula model
private Func<double, double> df; // discounting method delegate
private CreditCurve cds; // credit 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, NAGCopula copula,
CreditCurve cds, Func<double, double> discountFactor)
{
this.k = kth;
this.n = simulations;
this.maturity = maturity;
this.copula = copula;
this.cds = cds;
this.df = discountFactor;
}
public void Process()
{
// request correlated random numbers from copula model
int[] seed = new int[1] { Math.Abs(Guid.NewGuid().GetHashCode()) };
double[,] randomArray = copula.Create(n, seed);
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 correct 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 += (df(t) * dt * p);
}
}
else
{
for (int i = 0; i < maturity; i++)
{
v += df(i + 1);
}
}
return v;
}
private double calculateDefaultLeg(double[,] defaultTimes, int nDefaults)
{
double v = 0.0;
if ((nDefaults > 0) && (nDefaults >= k))
{
v = (1 - cds.Recovery) * df(defaultTimes[0, k - 1]) * (1 / (double)m);
}
return v;
}
}
//
//
// *************************************************************
public class CreditCurve
{
private Func<double, double> df;
private double recovery;
private double[,] CDSSpreads;
private double[,] survivalMatrix;
private double[,] hazardMatrix;
private double[,] cumulativeHazardMatrix;
//
public CreditCurve(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]);
}
}
}
}
//
//
// *************************************************************
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;
}
}
}
//
//
// *************************************************************
// NAG G05 COPULAS WRAPPER
public class NAGCopula
{
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 double[,] Create(int n, int[] seed)
{
result = new double[n, m];
//
G05.G05State g05State = new G05.G05State(genID, subGenID, 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 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();
}
}
}
