diff options
Diffstat (limited to 'tranche_functions.R')
| -rw-r--r-- | tranche_functions.R | 280 |
1 files changed, 280 insertions, 0 deletions
diff --git a/tranche_functions.R b/tranche_functions.R new file mode 100644 index 00000000..54290561 --- /dev/null +++ b/tranche_functions.R @@ -0,0 +1,280 @@ +library(statmod)
+
+lossdistrib <- function(p){
+ #basic recursive algorithm of Andersen, Sidenius and Basu
+ n <- length(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] <- (1-p[i])*q[1]
+ }
+ return(q)
+}
+
+lossdistrib.fft <- function(p){
+ #computes loss distribution using the fft
+ #complexity is of order O(n*m)+O(m*log(m))
+ #where m is the size of the grid and n the number of probabilities.
+ #this is slower than the recursive algorithm
+ theta <- 2*pi*1i*(0:n)/(n+1)
+ Phi <- 1 - p + p%o%exp(theta)
+ v <- apply(Phi, 2, prod)
+ return(1/(n+1)*Re(fft(v)))
+}
+
+lossdistrib2 <- function(p, S, lu){
+ #recursive algorithm with first order correction
+ #p vector of default probabilities
+ #S vector of severities
+ #lu loss unit
+ n <- length(p)
+ N <- 1/lu + 1
+ eps <- 1e-10
+ q <- rep(0, N+eps)
+ q[1] <- 1
+ for(i in 1:n){
+ d <- S[i]/(n*lu)
+ d1 <- floor(d)
+ d2 <- ceiling(d)
+ p1 <- p[i]*(d2-d)
+ p2 <- p[i] - p1
+ q1 <- c(rep(0,d1), p1*q[1:(N-d1)])
+ 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)
+ return(q)
+}
+
+lossdistrib.joint <- function(p, S, lu){
+ #recursive algorithm with first order correction
+ #to compute the joint probability distribition of the loss and recovery
+ #p vector of default probabilities
+ #S vector of severities
+ #lu loss unit
+ n <- length(p)
+ N <- ceiling(1/lu + 1)
+ q <- matrix(0, N, N)
+ q[1,1] <- 1
+ for(k in 1:n){
+ x <- S[k]/(n*lu)
+ y <- (1-S[k])/(n*lu)
+ i <- floor(x)
+ j <- floor(y)
+ weights <- c((i+1-x)*(j+1-y), (i+1-x)*(y-j), (x-i)*(y-j), (j+1-y)*(x-i))
+ psplit <- p[k] * weights
+ qtemp <- matrix(0, N, N)
+ qtemp[(i+1):N,(j+1):N] <- qtemp[(i+1):N,(j+1):N] + psplit[1] * q[1:(N-i),1:(N-j)]
+ qtemp[(i+1):N,(j+2):N] <- qtemp[(i+1):N,(j+2):N] + psplit[2] * q[1:(N-i), 1:(N-j-1)]
+ qtemp[(i+2):N,(j+2):N] <- qtemp[(i+2):N,(j+2):N] + psplit[3] * q[1:(N-i-1), 1:(N-j-1)]
+ qtemp[(i+2):N, (j+1):N] <- qtemp[(i+2):N, (j+1):N] + psplit[4] * q[1:(N-i-1), 1:(N-j)]
+ q <- qtemp + (1-p[k])*q
+ }
+ q[length(q)] <- q[length(q)]+1-sum(q)
+ return(q)
+}
+
+lossdistribC <- function(p, S, lu){
+ #C version of lossdistrib2, roughly 50 times faster
+ dyn.load("lossdistrib.dll")
+ .C("lossdistrib", as.double(p), as.integer(length(p)),
+ as.double(S), as.double(lu), q = double(1/lu+1))$q
+}
+
+lossdistribC.joint <- function(p, S, lu){
+ #C version of lossdistrib.joint, roughly 20 times faster
+ N <- ceiling(1/lu+1)
+ dyn.load("lossdistrib.dll")
+ .C("lossdistrib2", as.double(p), as.integer(length(p)), as.double(S),
+ as.double(lu), q = matrix(0, N, N))$q
+}
+
+lossdistrib3 <- function(p, S, lu){
+ #p vector of default probabilities
+ #S vector of severities times notional weights
+ #lu loss unit
+ n <- length(p)
+ N <- 1/lu + 1
+ eps <- 1e-10
+ q <- rep(0, N+eps)
+ q[1] <- 1
+ for(i in 1:n){
+ d <- S[i]/lu
+ d1 <- floor(d)
+ d2 <- ceiling(d)
+ p1 <- p[i]*(d2-d)
+ p2 <- p[i] - p1
+ q1 <- c(rep(0,d1), p1*q[1:(N-d1)])
+ q2 <- c(rep(0,d2), p2*q[1:(N-d2)])
+ q <- q1 + q2 + (1-p[i])*q
+ }
+ return(q)
+}
+
+lossrecovdist <- function(defaultprob, prepayprob, S, lu, useC=TRUE){
+ if(prepayprob==0){
+ if(useC){
+ L <- lossdistribC(defaultprob, S, lu)
+ R <- lossdistribC(defaultprob, 1-S, lu)
+ }else{
+ L <- lossdistrib2(defaultprob, S, lu)
+ R <- lossdistrib2(defaultprob, 1-S, lu)
+ }
+ }else{
+ u <- prepayprob/defaultprob
+ L <- lossdistribC(defaultprob+prepayprob, S/(1+u), lu)
+ R <- lossdistribC(defaultprob+prepayprob, (u+1-S)/(1+u), lu)
+ }
+ return(list(L=L, R=R))
+}
+
+lossrecovdist.term <- function(defaultprob, prepayprob, S, lu, 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)))
+ if(prepayprob==0){
+ for(t in 1:ncol(defaultprob)){
+ temp <- lossrecovdist(defaultprob[,t], 0, S[,t], lu, useC)
+ L[,t] <- temp$L
+ R[,t] <- temp$R
+ }
+ }
+ return(list(L=L, R=R))
+}
+
+shockprob <- function(p, rho, Z, log.p=F){
+ ## computes the shocked default probability as a function of the copula factor
+ pnorm((qnorm(p)-sqrt(rho)*Z)/sqrt(1-rho), log.p=log.p)
+}
+
+shockseverity <- function(S, Stilde=1, Z, rho, p){
+ #computes the severity as a function of the copula factor Z
+ Stilde * exp( shockprob(S/Stilde*p, rho, Z, TRUE) - shockprob(p, rho, Z, TRUE))
+}
+
+dshockprob <- function(p,rho,Z){
+ dnorm((qnorm(p)-sqrt(rho)*Z)/sqrt(1-rho))*dqnorm(p)/sqrt(1-rho)
+}
+
+dqnorm <- function(x){
+ 1/dnorm(qnorm(x))
+}
+
+fit.prob <- function(Z, w, rho, p0){
+ #if the weights are not perfectly gaussian, find the probability p such
+ # E_w(shockprob(p, rho, Z)) = p0
+ eps <- 1e-12
+ dp <- (crossprod(shockprob(p0,rho,Z),w)-p0)/crossprod(dshockprob(p0,rho,Z),w)
+ p <- p0
+ while(abs(dp) > eps){
+ dp <- (crossprod(shockprob(p,rho,Z),w)-p0)/crossprod(dshockprob(p,rho,Z),w)
+ phi <- 1
+ while ((p-phi*dp)<0 || (p-phi*dp)>1){
+ phi <- 0.8*phi
+ }
+ p <- p - phi*dp
+ }
+ return(p)
+}
+
+stochasticrecov <- function(R, Rtilde, Z, w, rho, porig, pmod){
+ ptilde <- fit.prob(Z, w, rho, (1-R)/(1-Rtilde) * porig)
+ return(abs(1-(1-Rtilde) * exp(shockprob(ptilde, rho, Z, TRUE) - shockprob(pmod, rho, Z, TRUE))))
+}
+
+pos <- function(x){
+ pmax(x, 0)
+}
+
+trancheloss <- function(L, K1, K2){
+ pos(L - K1) - pos(L - K2)
+}
+
+trancherecov <- function(R, K1, K2){
+ pos(R - 1 + K2) - pos(R - 1 +K1)
+}
+
+tranche.cl <- function(L, R, cs, K1, K2){
+ ## computes the couponleg of a tranche
+ if(K1==K2){
+ return( 0 )
+ }else{
+ n <- nrow(L)
+ support <- seq(0, 1, length=n)
+ accrued <- rep(0, length(cs$df))
+ for(t in 1:length(cs$df)){
+ accrued[t] <- (K2 - K1 - crossprod(trancheloss(support, K1, K2), L[,t]) +
+ crossprod(trancherecov(support, K1, K2), R[,t]))
+ }
+ return( crossprod(accrued * cs$coupons, cs$df)/(K2 - K1) )
+ }
+}
+
+tranche.pl <- function(L, R, cs, K1, K2){
+ ## computes the protection leg of a tranche
+ if(K1==K2){
+ return(0)
+ }else{
+ n <- nrow(L)
+ support <- seq(0, 1, length=n)
+ cf <- rep(0, length(cs$df))
+ for(t in 1:length(cs$df)){
+ cf[t] <- K2 - K1 - crossprod(trancheloss(support, K1, K2), L[,t])
+ }
+ cf <- c(1, cf / (K2-K1))
+ return( crossprod(diff(cf), cs$df))
+ }
+}
+
+tranche.bp <- function(L, R, cs, K1, K2){
+ return( 100 * ( 1 + tranche.pl(L, R, cs, K1, K2) + tranche.cl(L, R, cs, K1, K2)) )
+}
+
+adjust.attachments <- function(K, losstodate, factor){
+ return( pmin(pmax((K-losstodate)/factor, 0),1) )
+}
+
+tranche.bpvec <- function(K, L, R, cs){
+ r <- rep(0, length(K)-1)
+ for(i in 1:(length(K)-1)){
+ r[i] <- tranche.bp(L, R, cs, K[i], K[i+1])
+ }
+ return( r )
+}
+
+BClossdist <- function(SurvProb, recov, rho, lu, 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))
+ 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, S.shocked, lu)
+ LZ[,i] <- temp$L
+ RZ[,i] <- temp$R
+ }
+ L[,t] <- LZ%*%w
+ R[,t] <- RZ%*%w
+ }
+ list(L=L, R=R)
+}
+
+BClossdistC <- function(SurvProb, recov, rho, lu, 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(recov), as.double(Z), as.double(w), as.integer(n.int), as.double(rho), as.double(lu), L=L, R=R)
+ return(list(L=r$L,R=r$R))
+}
|