# n-th to default swap with five reference names

I would like to price a n-th to default swap on a basket of 5 assets or reference names. I started to code in R and I put the routine hereby. my doubt is how to use the m = {m1,m2,m3,m4,m5} series which represents the simulated defaults for each reference. Many Thanks for any tips you can give me if you have already worked for a similar project.

Many thanks for the support,

## Eugenio


library(moments)
library(psych)
library(plyr)
library(copula)
library(MASS)
# Importing the times series of the ref. names
# Calculating the weekly returns and the standardised returns
GS_r   = diff(log(GS[,6]))
BAC_r  = diff(log(BAC[,6]))
JPM_r  = diff(log(JPM[,6]))
CIT_r  = diff(log(CIT[,6]))
WFC_r  = diff(log(WFC[,6]))
Z_GS   = (GS_r -mean(GS_r))/sd(GS_r)
Z_BAC  = (BAC_r-mean(BAC_r))/sd(BAC_r)
Z_JPM  = (JPM_r-mean(JPM_r))/sd(JPM_r)
Z_CIT  = (CIT_r-mean(CIT_r))/sd(CIT_r)
Z_WFC  = (WFC_r-mean(WFC_r))/sd(WFC_r)
TAU    = cor(cbind(Z_GS, Z_BAC, Z_JPM, Z_CIT, Z_WFC), method = "kendall")
RHO    = sin((pi/2)*TAU)
# Credit spreads for the 5 Reference names from the market
CS_GS        = c(150, 155, 160, 170, 190)
CS_BAC       = c(20, 25, 32, 45, 57)
CS_JPM       = c(30, 45, 59, 67, 89)
CS_CIT       = c(100, 150, 155, 175, 210)
CS_WFC       = c(50, 59, 80, 120, 180)
CS           = matrix(c(CS_GS,CS_BAC,CS_JPM,CS_CIT,CS_WFC), nrow = 5, ncol = 5)

# Calculate the survival probabilities Matrix

ISP_mat  = function(CS) {
RR       = 0.40
df       = c(1,1.24,1.41,1.60,1.77,1.93) # Discounted Factors curve
nrows    = length(df)
ncols    = nrows-1
DF       = matrix(df,nrows,ncols)
ISP      = matrix(0,nrows,ncols)
h        = matrix(0,nrows,ncols)
terms    = matrix(0,nrows,ncols)
term     = matrix(0,nrows,ncols)
quotient = matrix(0,nrows,ncols)
firstTerm    = matrix(0,nrows,ncols)
lastTerm     = matrix(0,nrows,ncols)
for (i in 1:nrows) {
for (j in 1:ncols) {
if (i==1) { ISP[i,j] = 1 }
if (i==2) { ISP[i,j] = (1-RR) / ((1-RR) + CS[i-1,j]/10000) }
if (i> 2) {
for (k in 1:(i-2) ) {
term[i,j] = DF[k,j]*((1-RR)*ISP[k,j]-((1-RR)+(CS[i-1,j]/10000))*ISP[k+1,j])
terms[i,j] = terms[i,j] + term[i,j]
}
quotient[i,j]  = (DF[i-1,j]*((1-RR)+1*CS[i-1,j]/10000))
firstTerm[i,j] = (terms[i,j] / quotient[i,j])
lastTerm[i,j]  = ISP[i-1,j]*(1-RR)/(1-RR+1*CS[i-1,j]/10000)
ISP[i,j]       = (firstTerm[i,j] + lastTerm[i,j])
}
}
}
return(ISP)
}

surv.prob = ISP_mat(CS)

# Calculate hazard matrix and cumulative hazard rates matrix

nrows.h = nrow(surv.prob) - 1
ncols.h = ncol(surv.prob)
hazard  = matrix(0,nrows.h, ncols.h)
cumulativeHazard = matrix(0, nrows.h, ncols.h)
for (i in 1:nrows.h) {
for (j in 1:ncols.h) {
cumulativeHazard[i,j] = - log(surv.prob[i+1,j])
if (i==1) {
hazard[i,j] = cumulativeHazard[i,j]
}
if (i>1) {
hazard[i,j] = cumulativeHazard[i,j] - cumulativeHazard[i-1,j]
}
}
}

# Calculate the Default Times Matrix

nSim   = 10000
nrows  = nSim
ncols  = 5
time   = c(0,1,2,3,4,5)
df     = c(1,1.24,1.41,1.60,1.77,1.93)
e      = matrix(c(time,df), nrow = length(df), ncol = 2)
t      = matrix(0,nrows,ncols)
tenor  = matrix(0,ncols,ncols)
dt           = matrix(0,nrows,ncols)
defaultTenor = matrix(0,nrows,ncols)
defaultTime  = matrix(0,nrows,ncols)
u = pnorm(matrix(rnorm(nrows*ncols), ncol = ncols) %*% chol(RHO))
for (i in 1:nrows) {
for (j in 1:ncols) {
t[i,j] = abs(log(1-u[i,j]))
for (k in 1:ncols) {
if (cumulativeHazard[k,j] >= t[i,j] ) {
tenor[k,j] = k-1
}
if (tenor[k,j] >= 1) {
dt[i,j] = -(1/hazard[k,j])*log((1-u[i,j])/(exp(-cumulativeHazard[k-1,j])))
defaultTenor[i,j] = tenor[k,j] + dt[i,j]
}
}
if (defaultTenor[i,j] > 5 ) {
defaultTenor[i,j] = 0
}
if (defaultTenor[i,j] < 1 & defaultTenor[i,j] > 0) {
defaultTenor[i,j] = 0.50
}
}
}

defaultTime = defaultTenor
DL = matrix(0,nrow(defaultTime),ncol(defaultTime))
PL = matrix(0,nrow(defaultTime),ncol(defaultTime))

m  = defaultTime
# I extract the Default Times for each asset removing the zeros and sorting
m1 = sort(m[,1][!m[,1]==0])
m2 = sort(m[,2][!m[,2]==0])
m3 = sort(m[,3][!m[,3]==0])
m4 = sort(m[,4][!m[,4]==0])
m5 = sort(m[,5][!m[,5]==0])

# Linear Interpolation of the default times
getDF1 = rep(NA, length(m1))
getDF2 = rep(NA, length(m2))
getDF3 = rep(NA, length(m3))
getDF4 = rep(NA, length(m4))
getDF5 = rep(NA, length(m5))
r1 = rep(NA, length(m1))
r2 = rep(NA, length(m2))
r3 = rep(NA, length(m3))
r4 = rep(NA, length(m4))
r5 = rep(NA, length(m5))

for (i in 1:length(m1)) {
for (k in 1:5) {
if (m1[i] >= e[k,1] & m1[i] <= e[k+1,1]) {
r1[i] = e[k,2] + (e[k+1,2]-e[k,2])*(m1[i]-k)
getDF1[i] = exp(-r1[i]*m1[i])
}
}
}
for (i in 1:length(m2)) {
for (k in 1:5) {
if (m2[i] >= e[k,1] & m2[i] <= e[k+1,1]) {
r2[i] = e[k,2] + (e[k+1,2]-e[k,2])*(m2[i]-k)
getDF2[i] = exp(-r2[i]*m2[i])
}
}
}
for (i in 1:length(m3)) {
for (k in 1:5) {
if (m3[i] >= e[k,1] & m3[i] <= e[k+1,1]) {
r3[i] = e[k,2] + (e[k+1,2]-e[k,2])*(m3[i]-k)
getDF3[i] = exp(-r3[i]*m3[i])
}
}
}
for (i in 1:length(m4)) {
for (k in 1:5) {
if (m4[i] >= e[k,1] & m4[i] <= e[k+1,1]) {
r4[i] = e[k,2] + (e[k+1,2]-e[k,2])*(m4[i]-k)
getDF4[i] = exp(-r4[i]*m4[i])
}
}
}
for (i in 1:length(m5)) {
for (k in 1:5) {
if (m5[i] >= e[k,1] & m5[i] <= e[k+1,1]) {
r5[i] = e[k,2] + (e[k+1,2]-e[k,2])*(m5[i]-k)
getDF5[i] = exp(-r5[i]*m5[i])
}
}
}

# I calculate the Default Leg for each asset

DL1 = 1/5*(1-RR)*getDF1
DL2 = 1/5*(1-RR)*getDF2
DL3 = 1/5*(1-RR)*getDF3
DL4 = 1/5*(1-RR)*getDF4
DL5 = 1/5*(1-RR)*getDF5

# now I calculate the Protection Leg or Premium Leg

PL1 = getDF1*m1*5/5
PL2 = PL1 + getDF2*(m2-m1)*4/5
PL3 = PL2 + getDF3*(m3-m2)*3/5
PL4 = PL3 + getDF4*(m4-m3)*2/5
PL5 = PL4 + getDF5*(m5-m4)*1/5

# I calculate the spread as follows:

s1 = mean(DL1)/mean(PL1)
s2 = mean(DL2)/mean(PL2)
s3 = mean(DL3)/mean(PL3)
s4 = mean(DL4)/mean(PL4)
s5 = mean(DL5)/mean(PL5)



as you can see I extracted the default times for every single asset (in this case 5 names). Now I have to calculate the spread of each instrument and in a continuous time. the key concept is to calculate the Default Leg and the Protection Leg. Default Leg is calculated for each reference as following:

DL = 1/5 x (1-RR) x (Discounted Factors interpolated)

the protection leg is calculated in a different for each reference name.

For the 1st default, is

PL(1) = (DF_Interpolated_1 x m1 ) x (5/5).

For the 2nd Default I have:

PL(2) = PL(1) + (DF_Interpolated_2 x (m2 - m1) x (4/5)).

for the 3rd Default I have

PL(3) = PL(2) + (DF_Interpolated_3 x (m3 - m2) x (3/5)).

for the 4th Default I have

PL(4) = PL(3) + (DF_Interpolated_4 x (m4 - m3) x (2/5)).

for the 5th Default I have:

PL(5) = PL(4) + (DF_Interpolated_5 x (m5 - m4) x (1/5)).

What it is not clear is how to calculate the difference between two vectors having different lengths, since in my first round of 10,000 simulations, length(m1) = 1388, length(m2) = 406, length(m3) = 689, length(m4) = 1493, length(m5) = 1255.

if you need more details about the pricing model, please let me know.

• Please ensure that all data necessary to run this code is provided. Are you sure that tenor[k, j] is not NA the moment the error occurs? Apr 4 '20 at 17:38
• Hi Bob, yes now the code works out correctly because as you stated tenor was set as a matrix of NA , whereas it has to be a matrix of Zeroes. Yes I can provide more data, but at this stage , I am pretty sure I can continue with the code. thanks a lot Apr 6 '20 at 21:41
• I got the point. But my If condition is based on an inequality between a matrix c and another matrix t. Anyway, in the original code the tenor is not a matrix but it is defined as an Integer , hence A Number . I will post the new code with more data and comments so you can have a big picture. Many thanks in advance. Apr 8 '20 at 12:07

You have to initialize the tenor matrix to zero instead of NA. If you compare a number with NA it will yield an NA. Evaluating NA like this in an if-statement gives the error described.