diff options
Diffstat (limited to 'R/tranche_functions.R')
| -rw-r--r-- | R/tranche_functions.R | 798 |
1 files changed, 798 insertions, 0 deletions
diff --git a/R/tranche_functions.R b/R/tranche_functions.R new file mode 100644 index 00000000..e1f0ab91 --- /dev/null +++ b/R/tranche_functions.R @@ -0,0 +1,798 @@ +library(statmod)
+
+## todo:
+## -investigate other ways to interpolate the random severities on the grid
+## I'm thinking that at eah severity that we add to the distribution, round it down
+## and keep track of the missing mass: namely if X_i=S_i w.p p_i, then add
+## X_i=lu*floor(S_i/lu) with probability p_i and propagate
+## X_{i+1}=S_{i+1}+(S_i-lu*floor(S_i/lu)) with the right probability
+## - investigate truncated distributions more (need to compute loss and recov distribution
+## separately, for the 0-10 equity tranche, we need the loss on the 0-0.1 support and
+## recovery with 0.1-1 support, so it's not clear that there is a big gain.
+## - do the correlation adjustments when computing the deltas since it seems to be
+## the market standard
+
+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[-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, w, S, N, defaultflag=FALSE){
+ ## 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
+ ## defaultflag if true computes the default distribution
+ n <- length(p)
+ lu <- 1/(N-1)
+ q <- rep(0, N)
+ q[1] <- 1
+ for(i in 1:n){
+ if(defaultflag){
+ d <- w[i] /lu
+ }else{
+ d <- S[i] * w[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
+ }
+ q[length(q)] <- q[length(q)]+1-sum(q)
+ return(q)
+}
+
+lossdistrib2.truncated <- function(p, w, S, N, cutoff=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 (for best accuracy should
+ ## be a multiple of the number of issuers)
+ ## cutoff 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.
+ ## for high severities, M can become bigger than N, and there is
+ ## some probability mass escaping.
+ 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, cutoff-d1)]
+ q2 <- p2*q[1:min(M, cutoff-d2)]
+ q[1:min(M, cutoff)] <- (1-p[i])*q[1:min(M, cutoff)]
+ q[(d1+1):min(M+d1, cutoff)] <- q[(d1+1):min(M+d1, cutoff)]+q1
+ q[(d2+1):min(M+d2, cutoff)] <- q[(d2+1):min(M+d2, cutoff)]+q2
+ M <- M+d2
+ }
+ 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]
+ ## = w[i]*(1-S[i]) w.p dp[i]
+ ## = w[i] w.p pp[i]
+ ## each non zero value v is interpolated on the grid as
+ ## the pair of values floor(v/lu) and ceiling(v/lu) so that
+ ## X_i has four non zero values
+ n <- length(dp)
+ 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)
+ d1u <- ceiling(d1)
+ d2 <- w[i] / lu
+ d2l <- floor(d2)
+ d2u <- ceiling(d2)
+ dp1 <- dp[i] * (d1u-d1)
+ dp2 <- dp[i] - dp1
+ pp1 <- pp[i] * (d2u - d2)
+ pp2 <- pp[i] - pp1
+ q1 <- c(rep(0, d1l), dp1 * q[1:(N-d1l)])
+ q2 <- c(rep(0, d1u), dp2 * q[1:(N-d1u)])
+ q3 <- c(rep(0, d2l), pp1 * q[1:(N-d2l)])
+ q4 <- c(rep(0, d2u), pp2 *q[1:(N-d2u)])
+ q <- q1+q2+q3+q4+(1-dp[i]-pp[i])*q
+ }
+ return(q)
+}
+
+lossdist.joint <- function(p, w, S, N, defaultflag=FALSE){
+ ## 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
+ ## N: number of tick sizes on the grid
+ ## defaultflag: if true computes the default distribution
+ ## 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)
+ lu <- 1/(N-1)
+ q <- matrix(0, N, N)
+ q[1,1] <- 1
+ for(k in 1:n){
+ if(defaultflag){
+ x <- w[k] / lu
+ }else{
+ x <- S[k] * w[k]/lu
+ }
+ y <- (1-S[k]) * w[k]/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)
+}
+
+lossdist.prepay.joint <- function(dp, pp, w, S, N, defaultflag=FALSE){
+ ## recursive algorithm with first order correction
+ ## to compute the joint probability distribition of the loss and recovery
+ ## inputs:
+ ## dp: vector of default probabilities
+ ## pp: vector of prepay probabilities
+ ## w: vector of issuer weights
+ ## S: vector of severities
+ ## N: number of tick sizes on the grid
+ ## defaultflag: if true computes the default
+ ## 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)
+ lu <- 1/(N-1)
+ q <- matrix(0, N, N)
+ q[1,1] <- 1
+ for(k in 1:n){
+ y1 <- (1-S[k]) * w[k]/lu
+ y2 <- w[k]/lu
+ j1 <- floor(y1)
+ j2 <- floor(y2)
+ if(defaultflag){
+ x <- y2
+ i <- j2
+ }else{
+ x <- y2-y1
+ i <- floor(x)
+ }
+
+ ## weights <- c((i+1-x)*(j1+1-y1), (i+1-x)*(y1-j1), (x-i)*(y1-j1), (j1+1-y1)*(x-i))
+ weights1 <- c((i+1-x)*(j1+1-y1), (i+1-x)*(y1-j1), (x-i)*(y1-j1), (j1+1-y1)*(x-i))
+ dpsplit <- dp[k] * weights1
+
+ if(defaultflag){
+ weights2 <- c((i+1-x)*(j2+1-y2), (i+1-x)*(y2-j2), (x-i)*(y2-j2), (j2+1-y2)*(x-i))
+ ppsplit <- pp[k] * weights2
+ }else{
+ ppsplit <- pp[k] * c(j2+1-y2, y2-j2)
+ }
+ qtemp <- matrix(0, N, N)
+ qtemp[(i+1):N,(j1+1):N] <- qtemp[(i+1):N,(j1+1):N] + dpsplit[1] * q[1:(N-i),1:(N-j1)]
+ qtemp[(i+1):N,(j1+2):N] <- qtemp[(i+1):N,(j1+2):N] + dpsplit[2] * q[1:(N-i), 1:(N-j1-1)]
+ qtemp[(i+2):N,(j1+2):N] <- qtemp[(i+2):N,(j1+2):N] + dpsplit[3] * q[1:(N-i-1), 1:(N-j1-1)]
+ qtemp[(i+2):N,(j1+1):N] <- qtemp[(i+2):N, (j1+1):N] + dpsplit[4] * q[1:(N-i-1), 1:(N-j1)]
+ if(defaultflag){
+ qtemp[(i+1):N,(j2+1):N] <- qtemp[(i+1):N,(j2+1):N] + ppsplit[1] * q[1:(N-i),1:(N-j2)]
+ qtemp[(i+1):N,(j2+2):N] <- qtemp[(i+1):N,(j2+2):N] + ppsplit[2] * q[1:(N-i), 1:(N-j2-1)]
+ qtemp[(i+2):N,(j2+2):N] <- qtemp[(i+2):N,(j2+2):N] + ppsplit[3] * q[1:(N-i-1), 1:(N-j2-1)]
+ qtemp[(i+2):N,(j2+1):N] <- qtemp[(i+2):N, (j2+1):N] + ppsplit[4] * q[1:(N-i-1), 1:(N-j2)]
+ }else{
+ qtemp[, (j2+1):N] <- qtemp[,(j2+1):N]+ppsplit[1]*q[,1:(N-j2)]
+ qtemp[, (j2+2):N] <- qtemp[,(j2+2):N]+ppsplit[2]*q[,1:(N-j2-1)]
+ }
+ q <- qtemp + (1-pp[k]-dp[k]) * q
+ }
+ q[length(q)] <- q[length(q)] + 1 - sum(q)
+ return(q)
+}
+
+lossdistC <- function(p, w, S, N, defaultflag=FALSE){
+ ## C version of lossdistrib2, roughly 50 times faster
+ if(!is.loaded("lossdistrib")){
+ dyn.load(paste0("lossdistrib", .Platform$dynlib.ext))
+ }
+ .C("lossdistrib", as.double(p), as.integer(length(p)),
+ as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), q = double(N))$q
+}
+
+lossdistC.truncated <- function(p, w, S, N, T=N){
+ ## C version of lossdistrib2, roughly 50 times faster
+ if(!is.loaded("lossdistrib_truncated")){
+ dyn.load(paste0("lossdistrib", .Platform$dynlib.ext))
+ }
+ dyn.load("lossdistrib.dll")
+ .C("lossdistrib_truncated", 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
+ if(!is.loaded("recovdist")){
+ dyn.load(paste0("lossdistrib", .Platform$dynlib.ext))
+ }
+ .C("recovdist", as.double(dp), as.double(pp), as.integer(length(dp)),
+ as.double(w), as.double(S), as.integer(N), q = double(N))$q
+}
+
+lossdistC.joint <- function(p, w, S, N, defaultflag=FALSE){
+ ## C version of lossdistrib.joint, roughly 20 times faster
+ if(!is.loaded("lossdistrib_joint")){
+ dyn.load(paste0("lossdistrib", .Platform$dynlib.ext))
+ }
+ .C("lossdistrib_joint", as.double(p), as.integer(length(p)), as.double(w),
+ as.double(S), as.integer(N), as.logical(defaultflag), q = matrix(0, N, N))$q
+}
+
+lossdistC.prepay.joint <- function(dp, pp, w, S, N, defaultflag=FALSE){
+ ## C version of lossdist.prepay.joint
+ if(!is.loaded("lossdistrib_prepay_joint")){
+ dyn.load(paste0("lossdistrib", .Platform$dynlib.ext))
+ }
+ r <- .C("lossdistrib_prepay_joint", as.double(dp), as.double(pp), as.integer(length(dp)),
+ as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), q=matrix(0, N, N))$q
+ gc()
+ return(r)
+}
+
+lossrecovdist <- function(defaultprob, prepayprob, w, S, N, defaultflag=FALSE, useC=TRUE){
+ if(missing(prepayprob)){
+ if(useC){
+ L <- lossdistC(defaultprob, w, S, N, defaultflag)
+ R <- lossdistC(defaultprob, w, 1-S, N)
+ }else{
+ L <- lossdistrib2(defaultprob, w, S, N, defaultflag)
+ R <- lossdistrib2(defaultprob, w, 1-S, N)
+ }
+ }else{
+ if(useC){
+ L <- lossdistC(defaultprob+defaultflag*prepayprob, w, S, N, defaultflag)
+ R <- recovdistC(defaultprob, prepayprob, w, S, N)
+ }else{
+ L <- lossdistrib2(defaultprob+defaultflag*prepayprob, w, S, N, defaultflag)
+ R <- recovdist(defaultprob, prepayprob, w, S, N)
+ }
+ }
+ return(list(L=L, R=R))
+}
+
+lossrecovdist.term <- function(defaultprob, prepayprob, w, S, N, defaultflag=FALSE, useC=TRUE){
+ ## computes the loss and recovery distribution over time
+ L <- array(0, dim=c(N, ncol(defaultprob)))
+ R <- array(0, dim=c(N, ncol(defaultprob)))
+ if(missing(prepayprob)){
+ for(t in 1:ncol(defaultprob)){
+ temp <- lossrecovdist(defaultprob[,t], , w, S[,t], N, defaultflag, useC)
+ L[,t] <- temp$L
+ R[,t] <- temp$R
+ }
+ }else{
+ for(t in 1:ncol(defaultprob)){
+ temp <- lossrecovdist(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag, useC)
+ L[,t] <- temp$L
+ R[,t] <- temp$R
+ }
+ }
+ return(list(L=L, R=R))
+}
+
+lossrecovdist.joint.term <- function(defaultprob, prepayprob, w, S, N, defaultflag=FALSE, useC=TRUE){
+ ## computes the joint loss and recovery distribution over time
+ Q <- array(0, dim=c(ncol(defaultprob), N, N))
+ if(useC){
+ if(missing(prepayprob)){
+ for(t in 1:ncol(defaultprob)){
+ Q[t,,] <- lossdistC.joint(defaultprob[,t], w, S[,t], N, defaultflag)
+ }
+ }else{
+ for(t in 1:ncol(defaultprob)){
+ Q[t,,] <- lossdistC.prepay.joint(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag)
+ }
+ }
+ }else{
+ if(missing(prepayprob)){
+ for(t in 1:ncol(defaultprob)){
+ Q[t,,] <- lossdist.joint(defaultprob[,t], w, S[,t], N, defaultflag)
+ }
+ }else{
+ for(t in 1:ncol(defaultprob)){
+ Q[t,,] <- lossdist.prepay.joint(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag)
+ }
+ }
+ }
+ gc()
+ return(Q)
+}
+
+dist.transform <- function(dist.joint){
+ ## compute the joint (D, R) distribution
+ ## from the (L, R) distribution using D = L+R
+ distDR <- array(0, dim=dim(dist.joint))
+ Ngrid <- dim(dist.joint)[2]
+ u <- seq(0, 1, length=Ngrid)
+ v <- seq(0, 1, length=Ngrid)
+ for(t in 1:ncol(dp)){
+ for(i in 1:Ngrid){
+ for(j in 1:Ngrid){
+ index <- i+j
+ if(index <= Ngrid){
+ distDR[t,index,j] <- distDR[t,index,j] + dist.joint[t,i,j]
+ }
+ }
+ }
+ }
+ return( distDR )
+}
+
+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
+ if(p0==0){
+ return(0)
+ }else{
+ 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, Ngrid=nrow(L), scaled=FALSE){
+ ## computes the couponleg of a tranche
+ ## if scaled is TRUE, scale it by the size of the tranche (K2-K1)
+ ## can make use of the fact that the loss and recov distribution are
+ ## truncated (in that case nrow(L) != Ngrid
+ if(K1==K2){
+ return( 0 )
+ }else{
+ support <- seq(0, 1, length=Ngrid)[1:nrow(L)]
+ size <- K2 - K1 - crossprod(trancheloss(support, K1, K2), L) -
+ crossprod(trancherecov(support, K1, K2), R)
+ sizeadj <- as.numeric(0.5 * (size + c(K2-K1, size[-length(size)])))
+ if(scaled){
+ return( 1/(K2-K1) * crossprod(sizeadj * cs$coupons, cs$df) )
+ }else{
+ return( crossprod(sizeadj * cs$coupons, cs$df) )
+ }
+ }
+}
+
+tranche.cl.scenarios <- function(l, r, cs, K1, K2, scaled=FALSE){
+ ## computes the couponleg of a tranche for one scenario
+ ## if scaled is TRUE, scale it by the size of the tranche (K2-K1)
+ ## can make use of the fact that the loss and recov distribution are
+ ## truncated (in that case nrow(L) != Ngrid
+ if(K1==K2){
+ return( 0 )
+ }else{
+ size <- K2 - K1 - trancheloss(l, K1, K2) - trancherecov(r, K1, K2)
+ sizeadj <- as.numeric(0.5 * (size + c(K2-K1, size[-length(size)])))
+ if(scaled){
+ return( 1/(K2-K1) * crossprod(sizeadj * cs$coupons, cs$df) )
+ }else{
+ return( crossprod(sizeadj * cs$coupons, cs$df) )
+ }
+ }
+}
+
+tranche.pl <- function(L, cs, K1, K2, Ngrid=nrow(L), scaled=FALSE){
+ ## computes the protection leg of a tranche
+ ## if scaled
+ if(K1==K2){
+ return(0)
+ }else{
+ support <- seq(0, 1, length=Ngrid)[1:nrow(L)]
+ cf <- K2 - K1 - crossprod(trancheloss(support, K1, K2), L)
+ cf <- c(K2 - K1, cf)
+ if(scaled){
+ return( 1/(K2-K1) * crossprod(diff(cf), cs$df))
+ }else{
+ return( crossprod(diff(cf), cs$df))
+ }
+ }
+}
+
+tranche.pl.scenarios <- function(l, cs, K1, K2, scaled=FALSE){
+ ## computes the protection leg of a tranche
+ ## if scaled
+ if(K1==K2){
+ return(0)
+ }else{
+ cf <- K2 - K1 - trancheloss(l, K1, K2)
+ cf <- c(K2 - K1, cf)
+ if(scaled){
+ return( 1/(K2-K1) * as.numeric(crossprod(diff(cf), cs$df)))
+ }else{
+ return( as.numeric(crossprod(diff(cf), cs$df)))
+ }
+ }
+}
+
+tranche.pv <- function(L, R, cs, K1, K2, Ngrid=nrow(L)){
+ return( tranche.pl(L, cs, K1, K2, Ngrid) + tranche.cl(L, R, cs, K1, K2, Ngrid))
+}
+
+tranche.pv.scenarios <- function(l, r, cs, K1, K2){
+ return( tranche.pl.scenarios(l, cs, K1, K2, TRUE) +
+ tranche.cl.scenarios(l, r, cs, K1, K2, TRUE))
+}
+
+
+adjust.attachments <- function(K, losstodate, factor){
+ ## computes the attachments adjusted for losses
+ ## on current notional
+ return( pmin(pmax((K-losstodate)/factor, 0),1) )
+}
+
+tranche.pvvec <- function(K, L, R, cs){
+ r <- rep(0, length(K)-1)
+ for(i in 1:(length(K)-1)){
+ r[i] <- 1/(K[i+1]-K[i]) * tranche.pv(L, R, cs, K[i], K[i+1])
+ }
+ return( r )
+}
+
+BClossdist <- function(SurvProb, issuerweights, recov, rho, N=length(recov)+1,
+ n.int=100){
+ quadrature <- gauss.quad.prob(n.int, "normal")
+ Z <- quadrature$nodes
+ w <- quadrature$weights
+ 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, N)
+ LZ[,i] <- temp$L
+ RZ[,i] <- temp$R
+ }
+ L[,t] <- LZ%*%w
+ R[,t] <- RZ%*%w
+ }
+ list(L=L, R=R)
+}
+
+BClossdistC <- function(SurvProb, issuerweights, recov, rho,
+ N=length(issuerweights)+1, n.int=100, defaultflag){
+ if(!is.loaded("BClossdist")){
+ dyn.load(paste0("lossdistrib", .Platform$dynlib.ext))
+ }
+ quadrature <- gauss.quad.prob(n.int, "normal")
+ Z <- quadrature$nodes
+ w <- quadrature$weights
+ 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.integer(N), as.logical(defaultflag), L=L, R=R)
+ return(list(L=r$L,R=r$R))
+}
+
+BCtranche.pv <- function(portfolio, index, coupon, K1, K2, rho1, rho2, N=length(portfolio)+1){
+ ## computes the protection leg, couponleg, and bond price of a tranche
+ ## in the base correlation setting
+ if(K1==0){
+ if(rho1!=0){
+ stop("equity tranche must have 0 lower correlation")
+ }
+ }
+ SurvProb <- SPmatrix(portfolio, index)
+ cs <- couponSchedule(nextIMMDate(Sys.Date()), index$maturity, "Q", "FIXED", coupon, 0)
+ recov <- sapply(portfolio, attr, "recovery")
+ 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, N)
+ if(rho1!=0){
+ 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])
+ pl2 <- tranche.pl(dist2$L, cs, 0, K[2])
+ pl1 <- tranche.pl(dist1$L, cs, 0, K[1])
+ return(list(pl=(pl2-pl1)/dK, cl=(cl2-cl1)/dK,
+ bp=100*(1+(pl2-pl1+cl2-cl1)/dK)))
+}
+
+BCtranche.delta <- function(portfolio, index, coupon, K1, K2, rho1, rho2, N=length(portolio)+1){
+ ## 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
+ ## this makes sense for bottom-up skews
+ eps <- 1e-4
+ portfolioplus <- portfolio
+ portfoliominus <- portfolio
+ for(i in 1:length(portfolio)){
+ portfolioplus[[i]]@curve@hazardrates <- portfolioplus[[i]]@curve@hazardrates * (1 + eps)
+ portfoliominus[[i]]@curve@hazardrates <- portfoliominus[[i]]@curve@hazardrates * (1- eps)
+ }
+ dPVindex <- indexpv(portfolioplus, index)$bp - indexpv(portfoliominus, index)$bp
+ if(K2==1){
+ dPVtranche <- BCtranche.pv(portfolioplus, index, coupon, 0, K1, 0, rho1, N)$bp -
+ BCtranche.pv(portfoliominus, index, coupon, 0, K1, 0, rho1, N)$bp
+ K1adj <- adjust.attachments(K1, index$loss, index$factor)
+ delta <- (1 - dPVtranche/(100*dPVindex) * K1adj)/(1-K1adj)
+ }else{
+ dPVtranche <- BCtranche.pv(portfolioplus, index, coupon, K1, K2, rho1, rho2, N)$bp -
+ BCtranche.pv(portfoliominus, index, coupon, K1, K2, rho1, rho2, N)$bp
+ delta <- dPVtranche/(100*dPVindex)
+ }
+ ## dPVindex <- BCtranche.pv(portfolioplus, index, coupon, 0, 1, 0, 0.5, lu)$bp-
+ ## BCtranche.pv(portfoliominus, index, coupon, 0, 1, 0, 0.5, lu)$bp
+ return( delta )
+}
+
+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], N)$pl)
+ }
+ Kmodified <- adjust.attachments(K, index$loss, index$factor)
+ return(cumsum(EL*diff(Kmodified))/sum(EL*diff(Kmodified)))
+}
+
+delta.factor <- function(K1, K2, index){
+ ## compute the factor to convert from delta on current notional to delta on original notional
+ ## K1 and K2 original strikes
+ factor <- (adjust.attachments(K2, index$loss, index$factor)
+ -adjust.attachments(K1, index$loss, index$factor))/(K2-K1)
+ return( factor )
+}
+
+MFupdate.prob <- function(Z, w, rho, defaultprob){
+ ## update the probabilites based on a non gaussian factor
+ ## distribution so that the pv of the cds stays the same.
+ p <- matrix(0, nrow(defaultprob), ncol(defaultprob))
+ for(i in 1:nrow(defaultprob)){
+ for(j in 1:ncol(defaultprob)){
+ p[i,j] <- fit.prob(Z, w, rho, defaultprob[i,j])
+ }
+ }
+ return( p )
+}
+
+MFlossrecovdist <- function(w, Z, rho, defaultprob, defaultprobmod, issuerweights, recov,
+ Ngrid=2*length(issuerweights)+1, defaultflag=FALSE){
+ ## computes the loss and recovery distribution using the modified factor distribution
+ n.credit <- length(issuerweights)
+ n.int <- length(w)
+ Rstoch <- array(0, dim=c(n.int, n.credit, ncol(defaultprob)))
+ for(t in 1:ncol(defaultprob)){
+ for(i in 1:n.credit){
+ Rstoch[,i,t] <- stochasticrecov(recov[i], 0, Z, w, rho, defaultprob[i,t], defaultprobmod[i,t])
+ }
+ }
+ parf <- function(i){
+ pshocked <- apply(defaultprobmod, 2, shockprob, rho=rho, Z=Z[i])
+ S <- 1 - Rstoch[i,,]
+ dist <- lossrecovdist.term(pshocked, 0, issuerweights, S, Ngrid, defaultflag)
+ }
+ L <- matrix(0, Ngrid, ncol(defaultprob))
+ R <- matrix(0, Ngrid, ncol(defaultprob))
+ for(i in 1:length(w)){
+ dist <- parf(i)
+ L <- L + dist$L * w[i]
+ R <- R + dist$R * w[i]
+ }
+ return( list(L=L, R=R) )
+}
+
+MFlossrecovdist.prepay <- function(w, Z, rho, defaultprob, defaultprobmod, prepayprob, prepayprobmod,
+ issuerweights, recov, Ngrid=2*length(issuerweights)+1, defaultflag=FALSE){
+ ## computes the loss and recovery distribution using the modified factor distribution
+ n.credit <- length(issuerweights)
+ n.int <- length(w)
+ Rstoch <- array(0, dim=c(n.int, n.credit, ncol(defaultprob)))
+ for(t in 1:ncol(defaultprob)){
+ for(i in 1:n.credit){
+ Rstoch[,i,t] <- stochasticrecov(recov[i], 0, Z, w, rho, defaultprob[i,t], defaultprobmod[i,t])
+ }
+ }
+ parf <- function(i){
+ dpshocked <- apply(defaultprobmod, 2, shockprob, rho=rho, Z=Z[i])
+ ppshocked <- apply(prepayprobmod, 2, shockprob, rho=rho, Z=-Z[i])
+ S <- 1 - Rstoch[i,,]
+ dist <- lossrecovdist.term(dpshocked, ppshocked, issuerweights, S, Ngrid, defaultflag)
+ }
+ L <- matrix(0, Ngrid, ncol(defaultprob))
+ R <- matrix(0, Ngrid, ncol(defaultprob))
+ for(i in 1:length(w)){
+ dist <- parf(i)
+ L <- L + dist$L * w[i]
+ R <- R + dist$R * w[i]
+ }
+ return( list(L=L, R=R) )
+}
+
+MFlossdist.joint <- function(cl, w, Z, rho, defaultprob, defaultprobmod, issuerweights, recov,
+ Ngrid=2*length(issuerweights)+1, defaultflag=FALSE){
+ ## rowSums(Q) is the loss/default distribution
+ ## colSums(Q) is the recovery distribution
+ ## so that recovery is the y axis and L/D is the x axis
+ ## if we use the persp function, losses is the axes facing us,
+ ## and R is the axis going away from us.
+ n.credit <- length(issuerweights)
+ n.int <- lenth(w)
+ Rstoch <- array(0, dim=c(n.int, n.credit, ncol(defaultprob)))
+ for(t in 1:ncol(defaultprob)){
+ for(i in 1:n.credit){
+ Rstoch[,i,t] <- stochasticrecov(recov[i], 0, Z, w, rho, defaultprob[i,t], defaultprobmod[i,t])
+ }
+ }
+ parf <- function(i){
+ pshocked <- apply(defaultprobmod, 2, shockprob, rho=rho, Z=Z[i])
+ S <- 1 - Rstoch[i,,]
+ dist <- lossrecovdist.joint.term(pshocked, 0, issuerweights, S, Ngrid, defaultflag)
+ gc()
+ return(dist)
+ }
+ temp <- parSapply(cl, 1:length(w), parf)
+ clusterCall(cl, gc)
+ Q <- array(0, dim=c(ncol(defaultprob), Ngrid, Ngrid))
+ for(i in 1:length(w)){
+ Q <- Q + w[i]*array(temp[,i], dim=c(ncol(defaultprob), Ngrid, Ngrid))
+ }
+ return( Q )
+}
+
+MFlossdist.prepay.joint <- function(cl, w, Z, rho, defaultprob, defaultprobmod,
+ prepayprob, prepayprobmod, issuerweights, recov,
+ Ngrid=2*length(issuerweights)+1, defaultflag=FALSE, n.chunks=4){
+ ## rowSums is the loss distribution
+ ## colSums is the recovery distribution
+ ## so that recovery is the y axis and L is the x axis
+ ## if we use the persp function, losses is the axes facing us,
+ ## and R is the axis going away from us.
+ n.credit <- length(issuerweights)
+ n.int <- length(w)
+ Rstoch <- array(0, dim=c(n.int, n.credit, ncol(defaultprob)))
+ for(t in 1:ncol(defaultprob)){
+ for(i in 1:n.credit){
+ Rstoch[,i,t] <- stochasticrecov(recov[i], 0, Z, w, rho, defaultprob[i,t], defaultprobmod[i,t])
+ }
+ }
+ parf <- function(i){
+ dpshocked <- apply(defaultprobmod, 2, shockprob, rho=rho, Z=Z[i])
+ ppshocked <- apply(prepayprobmod, 2, shockprob, rho=rho, Z=-Z[i])
+ S <- 1 - Rstoch[i,,]
+ dist <- lossrecovdist.joint.term(dpshocked, ppshocked, issuerweights, S, Ngrid, defaultflag)
+ return(dist)
+ }
+ Q <- array(0, dim=c(ncol(defaultprob), Ngrid, Ngrid))
+ ## when using the cluster this takes a crazy amount of memory
+ ## we cut it in chunks to keep it manageable
+ chunksize <- length(w)/n.chunks
+ for(chunk in 1:n.chunks){
+ if(length(cl)>0){
+ clusterExport(cl, list("Rstoch", "defaultprobmod", "prepayprobmod"), envir=environment())
+ temp <- parSapply(cl, (1+chunksize*(chunk-1)):(chunksize*chunk), parf)
+ ## need to call gc twice for some reason, otherwise doesn't release everything
+ clusterCall(cl, gc)
+ clusterCall(cl, gc)
+ }else{
+ temp <- sapply((1+chunksize*(chunk-1)):chunksize*chunk, parf)
+ }
+ for(i in 1:chunksize){
+ Q <- Q + w[i+chunksize*(chunk-1)]*array(temp[,i], dim=c(ncol(defaultprob), Ngrid, Ngrid))
+ }
+ ## same here
+ gc()
+ gc()
+ }
+ return( Q )
+}
+
+MFtranche.pv <- function(cl, cs, w, rho, defaultprob, defaultprobmod, issuerweights, recov,
+ Kmodified, Ngrid=length(issuerweights)+1){
+ ## computes the tranches pv using the modified factor distribution
+ ## p is the modified probability so that
+ n.credit <- length(issuerweights)
+ Rstoch <- array(0, dim=c(n.int, n.credit, ncol(defaultprob)))
+ for(t in 1:ncol(defaultprob)){
+ for(i in 1:n.credit){
+ Rstoch[,i,t] <- stochasticrecov(recov[i], 0, Z, w, rho, defaultprob[i,t], defaultprobmod[i,t])
+ }
+ }
+ parf <- function(i){
+ pshocked <- apply(defaultprobmod, 2, shockprob, rho=rho, Z=Z[i])
+ S <- 1 - Rstoch[i,,]
+ dist <- lossrecovdist.term(pshocked, 0, issuerweights, S, Ngrid)
+ return( tranche.pvvec(Kmodified, dist$L, dist$R, cs))
+ }
+ clusterExport(cl, list("Rstoch", "p"), envir=environment())
+ result <- parSapply(cl, 1:length(w), parf)
+ return( list(pv=100*(1+result%*%w), pv.w=result))
+}
|