R/t38.R

In YafeiXu/CDO: This package is designed to reproduce the iTraxx Europe Index Series 8 tranche spreads.

t38 <-
function(theta1Input,theta2Input,dF,MInput,dateInput){
# t38(.1,.3,13,10,c("2008-04-04"))
M=MInput # Monte Carlo runs 
d=125 # # of entities
#rho=rhoInput # correlation parameter in copula 
U=matrix(NA,d,M)
norm.cop=list()
tao=matrix(NA,d,M)
# Settlement Date  
CompDate=list()
CompDate=as.Date(dateInput) # Settlement Date 
#lambda=lambdaInput
lambda=numeric()

lambdaVector=read.csv("C:/defIntensity.csv")
lambda=lambdaVector[which(as.Date(lambdaVector[,1])==as.Date(CompDate)),3]

diffTime=numeric()
timeTable=data.frame()
Time=numeric()
TT=numeric()
BETA=numeric()
payDay=read.csv("C:/payday.csv")[,1]
BETA=numeric()

# Flat Recovery Rate		40% 
R=numeric()
R=0.4#recoveryInput
# Risk-free interest rate		5.0%  
IntRate=numeric()
IntRate=0.03#intRateInput

################################################################################
#---check needed packages  
is.installed=list()
is.installed <- function(mypkg) is.element(mypkg, installed.packages()[,1])

if(is.installed("copula")){library(copula)
}else{
install.packages("copula")
library(copula)}

if(is.installed("matrixcalc")){library(matrixcalc)
}else{
install.packages("matrixcalc")
library(matrixcalc)}


#---1) 2)Compute tao 
rho1=theta1Input
rho2=theta2Input
s1=matrix(rho2,30,30)+(1-rho2)*diag(30)
s2=matrix(rho2,25,25)+(1-rho2)*diag(25)
s3=matrix(rho2,20,20)+(1-rho2)*diag(20)
s4=matrix(rho2,20,20)+(1-rho2)*diag(20)
s5=matrix(rho2,20,20)+(1-rho2)*diag(20)
s6=matrix(rho2,10,10)+(1-rho2)*diag(10)
vsigma=matrix(rho1,125,125)
vsigma[1:30,1:30]=s1
vsigma[31:55,31:55]=s2
vsigma[56:75,56:75]=s3
vsigma[76:95,76:95]=s4
vsigma[96:115,96:115]=s5
vsigma[116:125,116:125]=s6
vsigma
dF=floor(dF)
#t1.cop=tCopula(dim = d, dispstr = "ex",df = dF, df.fixed = TRUE)
#---1) 2)Compute tao  
t.cop <- tCopula(P2p(vsigma), dim = 125, df = dF, disp = "un", df.fixed = TRUE)  
#U <- rCopula(M, norm.cop) 
#U=t(rCopula(M, norm.cop))

#tao=(-1/lambda)*log(t(rCopula(M, t.cop)))
#iTau(t.cop,rho)
#t.cop <- tCopula(rho, dim = d, dispstr = "ex",df = dF, df.fixed = TRUE)
#U <- rCopula(M, t.cop)  
#U=t(rCopula(M, t.cop))

tao=(-1/lambda)*log(t(rCopula(M, t.cop)))

#---3) 4)Compute Time  
diffTime=as.numeric(as.Date(payDay)-as.Date(CompDate))
timeTable=data.frame(as.Date(payDay),diffTime)
if(timeTable[1,2]<0){timeTable=timeTable[-which(timeTable[,2]<0),]}else{timeTable=timeTable}

Time=timeTable[,1] # term structure  
#Time=as.numeric(as.Date(pa)-as.Date(tStStart  
TT=length(timeTable[,1])  
tPercent=numeric()
tPercent=cumsum(as.numeric(as.Date(c(CompDate,Time))[-1]-as.Date(c(CompDate,Time))[-length(c(CompDate,Time))])/365)
ONE=matrix(NA,125*M,TT)



#---6) Compute ONE matrix   
taoVector=matrix(NA,125*M,1)
taoVector=vec(tao)
ONE=apply(taoVector,1,function(x)as.numeric(x<=tPercent))
ONE=t(ONE)

#---5) Compute BETA  
BETA=exp(-IntRate*(tPercent))

#---6) L matrix  
L=matrix(NA,M,TT)
for(i in 1:M){
   for(j in 1:TT){
      L[i,j]=sum(ONE[(i*125-124):(i*125),j])
   }
}
L=((1-R)/d)*L

#---7) Lj matrix 
l1=numeric()
u1=numeric()
l2=numeric()
u2=numeric()
l3=numeric()
u3=numeric()
l4=numeric()
u4=numeric()
l5=numeric()
u5=numeric()
l1=0
u1=0.03
l2=0.03
u2=0.06
l3=0.06
u3=0.09
l4=0.09
u4=0.12
l5=0.12
u5=0.22
fct=function(x){

L1=numeric()
L2=numeric()
L3=numeric()
L4=numeric()
L5=numeric()

L1=0.03*as.numeric(x>u1)+(x-l1)*as.numeric(x>l1 & x<=u1)+0*as.numeric(x<=l1)
L2=0.03*as.numeric(x>u2)+(x-l2)*as.numeric(x>l2 & x<=u2)+0*as.numeric(x<=l2)
L3=0.03*as.numeric(x>u3)+(x-l3)*as.numeric(x>l3 & x<=u3)+0*as.numeric(x<=l3)
L4=0.03*as.numeric(x>u4)+(x-l4)*as.numeric(x>l4 & x<=u4)+0*as.numeric(x<=l4)
L5=0.1*as.numeric(x>u5)+(x-l5)*as.numeric(x>l5 & x<=u5)+0*as.numeric(x<=l5)
return(c(L1,L2,L3,L4,L5))
}

LjResult=matrix(NA,5,M*TT)
LjResult=sapply(t(L),fct)
LjResult=t(LjResult)

#---8) LjResultMinus1
minusIndex=numeric()
LjResultTemp=matrix(NA,M*TT,5)
minusIndex=c(1:M)*TT
LjResultTemp=LjResult
LjResultTemp[minusIndex,]=0
LjResultMinus1=LjResultTemp
LjResultMinus1=LjResultMinus1[-length(LjResultMinus1[,1]),]
LjResultMinus1=rbind(seq(0,0,length.out=5),LjResultMinus1)

#---9) LjDiff
LjDiff=matrix(NA,M*TT,5)
LjDiff=LjResult-LjResultMinus1

#---10) BETAMatrix
betaMatrix=matrix(NA,M*TT,5)
betaMatrix=matrix(BETA,M*TT,5)

#---11) 12)upStar, upSumStar, eUpSumStar
upStar=matrix(NA,M*TT,5)
upSumStar=numeric()
eUpSumStar=numeric()
upStar=betaMatrix*LjDiff
upSumStar=colSums(upStar)
eUpSumStar=upSumStar*(1/M)

#---13) FjResult
FjResult=matrix(NA,M*TT,5)
FjResult=cbind((u1-l1)-LjResult[,1],(u2-l2)-LjResult[,2],(u3-l3)-LjResult[,3],(u4-l4)-LjResult[,4],(u5-l5)-LjResult[,5])

#---14) FjResultMinus1
FjResultMinus1=matrix(NA,M*TT,5)
FjResultTemp=matrix(NA,M*TT,5)
FjResultTemp=FjResult
FjResultTemp[minusIndex,]=0
FjResultMinus1=FjResultTemp[-length(FjResult[,1]),]
FjResultMinus1=rbind(seq(0,0,length.out=5),FjResultMinus1)

#---15) FjAdd
FjAdd=matrix(NA,M*TT,5)
FjAddHalf=matrix(NA,M*TT,5)
FjAdd=FjResult+FjResultMinus1
FjAddHalf=0.5*FjAdd

#---16) deltaMatrix
deltaT=numeric()
deltaT=as.numeric(as.Date(c(CompDate,Time))[-1]-as.Date(c(CompDate,Time))[-length(c(CompDate,Time))])/365
deltaMatrix=matrix(NA,M*TT,5)
deltaMatrix=matrix(deltaT,M*TT,5)

#---17) lowStar, lowSumStar, eLowSumStar
lowStar=matrix(NA,M*TT,5)
lowSumStar=numeric()
eLowSumStar=numeric()
lowStar=FjAddHalf*betaMatrix*deltaMatrix
lowSumStar= colSums(lowStar)
eLowSumStar=(1/M)*lowSumStar

#---18) tranches' spreads
trancheSpread1=(eUpSumStar[1]-0.05*eLowSumStar[1])/0.03
trancheSpread2=eUpSumStar[2]/eLowSumStar[2]
trancheSpread3=eUpSumStar[3]/eLowSumStar[3]
trancheSpread4=eUpSumStar[4]/eLowSumStar[4]
trancheSpread5=eUpSumStar[5]/eLowSumStar[5]

#---19) return value
return(c(trancheSpread1,trancheSpread2,trancheSpread3,trancheSpread4,trancheSpread5))
}