aboutsummaryrefslogtreecommitdiffstats
path: root/tranche_functions.R
diff options
context:
space:
mode:
Diffstat (limited to 'tranche_functions.R')
-rw-r--r--tranche_functions.R184
1 files changed, 130 insertions, 54 deletions
diff --git a/tranche_functions.R b/tranche_functions.R
index a7262ea0..69b14511 100644
--- a/tranche_functions.R
+++ b/tranche_functions.R
@@ -6,7 +6,7 @@ lossdistrib <- function(p){
q <- rep(0, (n+1))
q[1] <- 1
for(i in 1:n){
- q[-(n+1)] <- p[i]*q[-(n+1)]+(1-p[i])*q[-1]
+ q[-1] <- p[i]*q[-(n+1)]+(1-p[i])*q[-1]
q[1] <- (1-p[i])*q[1]
}
return(q)
@@ -23,16 +23,15 @@ lossdistrib.fft <- function(p){
return(1/(n+1)*Re(fft(v)))
}
-lossdistrib2 <- function(p, w, S, lu){
+lossdistrib2 <- function(p, w, S, N){
#recursive algorithm with first order correction
#p vector of default probabilities
#w vector of weigths
#S vector of severities
- #lu loss unit
+ #N number of ticks in the grid
n <- length(p)
- N <- 1/lu + 1
- eps <- 1e-10
- q <- rep(0, N+eps)
+ lu <- 1/(N-1)
+ q <- rep(0, N)
q[1] <- 1
for(i in 1:n){
d <- S[i] * w[i] / lu
@@ -44,11 +43,75 @@ lossdistrib2 <- function(p, w, S, lu){
q2 <- c(rep(0,d2), p2*q[1:(N-d2)])
q <- q1 + q2 + (1-p[i])*q
}
- q[length(q)] <- q[length(q)]+1-sum(q)
+ ## q[length(q)] <- q[length(q)]+1-sum(q)
+ return(q)
+}
+
+lossdistrib2.fast <- function(p, w, S, N){
+ #recursive algorithm with first order correction
+ #p vector of default probabilities
+ #w vector of weigths
+ #S vector of severities
+ #N number of ticks in the grid
+ ## this is actually slower than lossdistrib2. But in C this is
+ ## twice as fast.
+ n <- length(p)
+ lu <- 1/(N-1)
+ q <- rep(0, N)
+ q[1] <- 1
+ M <- 1
+ for(i in 1:n){
+ d <- S[i] * w[i] / lu
+ d1 <- floor(d)
+ d2 <- ceiling(d)
+ p1 <- p[i]*(d2-d)
+ p2 <- p[i] - p1
+ q1 <- p1*q[1:M]
+ q2 <- p2*q[1:M]
+ q[1:M] <- (1-p[i])*q[1:M]
+ q[(d1+1):(M+d1)] <- q[(d1+1):(M+d1)]+q1
+ q[(d2+1):(M+d2)] <- q[(d2+1):(M+d2)]+q2
+ M <- M+d2
+ }
+ ## q[length(q)] <- q[length(q)]+1-sum(q)
return(q)
}
-recovdist <- function(dp, pp, w, S, lu){
+lossdistrib2.fast.trunc <- function(p, w, S, N, truncated=1){
+ ## recursive algorithm with first order correction
+ ## p vector of default probabilities
+ ## w vector of weigths
+ ## S vector of severities
+ ## N number of ticks in the grid
+ ## truncated: where to stop computing the exact probabilities
+ ## (useful for tranche computations)
+
+ ## this is actually slower than lossdistrib2. But in C this is
+ ## twice as fast.
+ n <- length(p)
+ lu <- 1/(N-1)
+ q <- rep(0, truncated)
+ q[1] <- 1
+ M <- 1
+ for(i in 1:n){
+ d <- S[i] * w[i] / lu
+ d1 <- floor(d)
+ d2 <- ceiling(d)
+ p1 <- p[i]*(d2-d)
+ p2 <- p[i] - p1
+ q1 <- p1*q[1:min(M, truncated-d1)]
+ q2 <- p2*q[1:min(M, truncated-d2)]
+ q[1:min(M, truncated)] <- (1-p[i])*q[1:min(M, truncated)]
+ q[(d1+1):min(M+d1, truncated)] <- q[(d1+1):min(M+d1, truncated)]+q1
+ q[(d2+1):min(M+d2, truncated)] <- q[(d2+1):min(M+d2, truncated)]+q2
+ M <- M+d2
+ }
+ ## q[length(q)] <- q[length(q)]+1-sum(q)
+ return(q)
+}
+
+
+recovdist <- function(dp, pp, w, S, N){
## computes the recovery distribution for a sum of independent variables
## R=\sum_{i=1}^n X_i
## where X_i = 0 w.p 1-dp[i]-pp[i]
@@ -58,9 +121,9 @@ recovdist <- function(dp, pp, w, S, lu){
## the pair of values floor(v/lu) and ceiling(v/lu) so that
## X_i has four non zero values
n <- length(dp)
- N <- 1/lu + 1
- q <- rep(0, ceiling(N))
+ q <- rep(0, N)
q[1] <- 1
+ lu <- 1/(N-1)
for(i in 1:n){
d1 <- w[i]*(1-S[i])/lu
d1l <- floor(d1)
@@ -82,20 +145,20 @@ recovdist <- function(dp, pp, w, S, lu){
return(q)
}
-lossdistrib.joint <- function(p, w, S, lu){
+lossdistrib.joint <- function(p, w, S, N){
## recursive algorithm with first order correction
## to compute the joint probability distribution of the loss and recovery
## inputs:
## p: vector of default probabilities
## w: vector of issuer weights
## S: vector of severities
- ## lu: loss unit
+ ## N: number of tick sizes on the grid
## output:
## q: matrix of joint loss, recovery probability
## colSums(q) is the recovery distribution marginal
## rowSums(q) is the loss distribution marginal
n <- length(p)
- N <- ceiling(1/lu + 1)
+ lu <- 1/(N-1)
q <- matrix(0, N, N)
q[1,1] <- 1
for(k in 1:n){
@@ -116,7 +179,7 @@ lossdistrib.joint <- function(p, w, S, lu){
return(q)
}
-lossdistribprepay.joint <- function(dp, pp, w, S, lu){
+lossdistribprepay.joint <- function(dp, pp, w, S, N){
## recursive algorithm with first order correction
## to compute the joint probability distribition of the loss and recovery
## inputs:
@@ -124,13 +187,13 @@ lossdistribprepay.joint <- function(dp, pp, w, S, lu){
## pp: vector of prepay probabilities
## w: vector of issuer weights
## S: vector of severities
- ## lu: loss unit
+ ## N: number of tick sizes on the grid
## outputs
## q: matrix of joint loss and recovery probability
## colSums(q) is the recovery distribution marginal
## rowSums(q) is the loss distribution marginal
n <- length(dp)
- N <- ceiling(1/lu + 1)
+ lu <- 1/(N-1)
q <- matrix(0, N, N)
q[1,1] <- 1
for(k in 1:n){
@@ -156,60 +219,64 @@ lossdistribprepay.joint <- function(dp, pp, w, S, lu){
return(q)
}
-
-lossdistribC <- function(p, w, S, lu){
+lossdistribC <- function(p, w, S, N){
## C version of lossdistrib2, roughly 50 times faster
dyn.load("lossdistrib.dll")
.C("lossdistrib", as.double(p), as.integer(length(p)),
- as.double(w), as.double(S), as.double(lu), q = double(1/lu+1))$q
+ as.double(w), as.double(S), as.integer(N), q = double(N))$q
}
-recovdistC <- function(dp, pp, w, S, lu){
+lossdistrib.fastC <- function(p, w, S, N, T=N){
+ ## C version of lossdistrib2, roughly 50 times faster
+ dyn.load("lossdistrib.dll")
+ .C("lossdistrib_fast", as.double(p), as.integer(length(p)),
+ as.double(w), as.double(S), as.integer(N), as.integer(T), q = double(T))$q
+}
+
+recovdistC <- function(dp, pp, w, S, N){
## C version of recovdist
dyn.load("lossdistrib.dll")
.C("recovdist", as.double(dp), as.double(pp), as.integer(length(dp)),
- as.double(w), as.double(S), as.double(lu), q = double(ceiling(1/lu+1)))$q
+ as.double(w), as.double(S), as.integer(N), q = double(N))$q
}
-lossdistribC.joint <- function(p, w, S, lu){
+lossdistribC.joint <- function(p, w, S, N){
## C version of lossdistrib.joint, roughly 20 times faster
- N <- ceiling(1/lu+1)
dyn.load("lossdistrib.dll")
.C("lossdistrib_joint", as.double(p), as.integer(length(p)), as.double(w),
- as.double(S), as.double(lu), q = matrix(0, N, N))$q
+ as.double(S), as.integer(N), q = matrix(0, N, N))$q
}
-lossdistribprepayC.joint <- function(dp, pp, w, S, lu){
+lossdistribprepayC.joint <- function(dp, pp, w, S, N){
## C version of lossdistribprepay.joint
dyn.load("lossdistrib.dll")
- N <- ceiling(1/lu+1)
.C("lossdistrib_prepay_joint", as.double(dp), as.double(pp), as.integer(length(dp)),
- as.double(w), as.double(S), as.double(lu), q=matrix(0, N, N))$q
+ as.double(w), as.double(S), as.integer(N), q=matrix(0, N, N))$q
}
-lossrecovdist <- function(defaultprob, prepayprob, w, S, lu, useC=TRUE){
+lossrecovdist <- function(defaultprob, prepayprob, w, S, N, useC=TRUE){
if(all(!prepayprob)){
if(useC){
- L <- lossdistribC(defaultprob, w, S, lu)
- R <- lossdistribC(defaultprob, w, 1-S, lu)
+ L <- lossdistribC(defaultprob, w, S, N)
+ R <- lossdistribC(defaultprob, w, 1-S, N)
}else{
- L <- lossdistrib2(defaultprob, w, S, lu)
- R <- lossdistrib2(defaultprob, w, 1-S, lu)
+ L <- lossdistrib2(defaultprob, w, S, N)
+ R <- lossdistrib2(defaultprob, w, 1-S, N)
}
}else{
- L <- lossdistribC(defaultprob, w, S, lu)
- R <- recovdistC(defaultprob, prepayprob, w, S, lu)
+ L <- lossdistribC(defaultprob, w, S, N)
+ R <- recovdistC(defaultprob, prepayprob, w, S, N)
}
return(list(L=L, R=R))
}
-lossrecovdist.term <- function(defaultprob, prepayprob, w, S, lu, useC=TRUE){
+lossrecovdist.term <- function(defaultprob, prepayprob, w, S, N, useC=TRUE){
## computes the loss and recovery distribution over time
- L <- array(0, dim=c(1/lu+1, ncol(defaultprob)))
- R <- array(0, dim=c(1/lu+1, ncol(defaultprob)))
+ L <- array(0, dim=c(N, ncol(defaultprob)))
+ R <- array(0, dim=c(N, ncol(defaultprob)))
if(all(!prepayprob)){
for(t in 1:ncol(defaultprob)){
- temp <- lossrecovdist(defaultprob[,t], 0, w, S[,t], lu, useC)
+ temp <- lossrecovdist(defaultprob[,t], 0, w, S[,t], N, useC)
L[,t] <- temp$L
R[,t] <- temp$R
}
@@ -324,21 +391,20 @@ tranche.bpvec <- function(K, L, R, cs){
return( r )
}
-BClossdist <- function(SurvProb, issuerweights, recov, rho, lu, n.int=100){
+BClossdist <- function(SurvProb, issuerweights, recov, rho, N, n.int=100){
quadrature <- gauss.quad.prob(n.int, "normal")
Z <- quadrature$nodes
w <- quadrature$weights
- n <- as.integer(1/lu+1)
- LZ <- matrix(0, n, n.int)
- RZ <- matrix(0, n, n.int)
- L <- matrix(0, n, ncol(SurvProb))
- R <- matrix(0, n, ncol(SurvProb))
+ LZ <- matrix(0, N, n.int)
+ RZ <- matrix(0, N, n.int)
+ L <- matrix(0, N, ncol(SurvProb))
+ R <- matrix(0, N, ncol(SurvProb))
for(t in 1:ncol(SurvProb)){
g <- 1 - SurvProb[, t]
for(i in 1:length(Z)){
g.shocked <- shockprob(g, rho, Z[i])
S.shocked <- shockseverity(1-recov, 1, Z[i], rho, g)
- temp <- lossrecovdist(g.shocked, 0, issuerweights, S.shocked, lu)
+ temp <- lossrecovdist(g.shocked, 0, issuerweights, S.shocked, N)
LZ[,i] <- temp$L
RZ[,i] <- temp$R
}
@@ -348,19 +414,18 @@ BClossdist <- function(SurvProb, issuerweights, recov, rho, lu, n.int=100){
list(L=L, R=R)
}
-BClossdistC <- function(SurvProb, issuerweights, recov, rho, lu, n.int=100){
+BClossdistC <- function(SurvProb, issuerweights, recov, rho, N, n.int=100){
dyn.load("lossdistrib.dll")
quadrature <- gauss.quad.prob(n.int, "normal")
Z <- quadrature$nodes
w <- quadrature$weights
- N <- as.integer(1/lu+1)
L <- matrix(0, N, dim(SurvProb)[2])
R <- matrix(0, N, dim(SurvProb)[2])
- r <- .C("BClossdist", SurvProb, as.integer(dim(SurvProb)[1]), as.integer(dim(SurvProb)[2]), as.double(issuerweights), as.double(recov), as.double(Z), as.double(w), as.integer(n.int), as.double(rho), as.double(lu), L=L, R=R)
+ r <- .C("BClossdist", SurvProb, as.integer(dim(SurvProb)[1]), as.integer(dim(SurvProb)[2]), as.double(issuerweights), as.double(recov), as.double(Z), as.double(w), as.integer(n.int), as.double(rho), as.integer(N), L=L, R=R)
return(list(L=r$L,R=r$R))
}
-BCtranche.pv <- function(portfolio, index, coupon, K1, K2, rho1, rho2, lu=0.01){
+BCtranche.pv <- function(portfolio, index, coupon, K1, K2, rho1, rho2, N=101){
## computes the protection leg, couponleg, and bond price of a tranche
## in the base correlation setting
if(K1==0){
@@ -374,9 +439,9 @@ BCtranche.pv <- function(portfolio, index, coupon, K1, K2, rho1, rho2, lu=0.01){
issuerweights <- rep(1/length(portfolio), length(portfolio))
K <- adjust.attachments(c(K1,K2), index$loss, index$factor)
dK <- K[2] - K[1]
- dist2 <- BClossdistC(SurvProb, issuerweights, recov, rho2, lu)
+ dist2 <- BClossdistC(SurvProb, issuerweights, recov, rho2, N)
if(rho1!=0){
- dist1 <- BClossdistC(SurvProb, issuerweights, recov, rho1, lu)
+ dist1 <- BClossdistC(SurvProb, issuerweights, recov, rho1, N)
}
cl2 <- tranche.cl(dist2$L, dist2$R, cs, 0, K[2])
cl1 <- tranche.cl(dist1$L, dist1$R, cs, 0, K[1])
@@ -386,7 +451,7 @@ BCtranche.pv <- function(portfolio, index, coupon, K1, K2, rho1, rho2, lu=0.01){
bp=100*(1+(pl2-pl1+cl2-cl1)/dK)))
}
-BCtranche.delta <- function(portfolio, index, coupon, K1, K2, rho1, rho2, lu=0.01){
+BCtranche.delta <- function(portfolio, index, coupon, K1, K2, rho1, rho2, N=101){
## computes the tranche delta (on current notional) by doing a proportional
## blip of all the curves
## if K2==1, then computes the delta using the lower attachment only
@@ -414,11 +479,11 @@ BCtranche.delta <- function(portfolio, index, coupon, K1, K2, rho1, rho2, lu=0.0
return( delta )
}
-BCstrikes <- function(portfolio, index, coupon, K, rho) {
+BCstrikes <- function(portfolio, index, coupon, K, rho, N=101) {
## computes the strikes as a percentage of expected loss
EL <- c()
for(i in 2:length(K)){
- EL <- c(EL, -BCtranche.pv(portfolio, index, coupon, K[i-1], K[i], rho[i-1], rho[i], lu)$pl)
+ EL <- c(EL, -BCtranche.pv(portfolio, index, coupon, K[i-1], K[i], rho[i-1], rho[i], N)$pl)
}
Kmodified <- adjust.attachments(K, index$loss, index$factor)
return(cumsum(EL*diff(Kmodified))/sum(EL*diff(Kmodified)))
@@ -431,3 +496,14 @@ delta.factor <- function(K1, K2, index){
-adjust.attachments(K1, index$loss, index$factor))/(K2-K1)
return( factor )
}
+
+rho <- 0.999
+N <- 101
+L <- matrix(0, N, 100)
+for(i in 1:100){
+ ptilde <- shockprob(p, rho, Z[i])
+ Sshocked <- shockseverity(S, 1, Z[i], rho, p)
+ L[,i] <- lossdistrib2(ptilde, w, Sshocked, N)
+}
+Lw <- L%*%quadrature$weights
+crossprod(Lw, seq(0, 1, length=N))
generated by cgit v1.2.3-70-g09d2 (git 2.52.0) at 2026-01-26 03:28:16 +0000