path: root/tranche_functions.R

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){
    #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
    n <- length(p)
    lu <- 1/(N-1)
    q <- rep(0, N)
    q[1] <- 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 <- 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
        cat(dp1, dp2, pp1, pp2, "\n")
        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)
}

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
    ##   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)
    lu <- 1/(N-1)
    q <- matrix(0, N, N)
    q[1,1] <- 1
    for(k in 1:n){
        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)
}

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:
    ##   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
    ## 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){
        x <- S[k] * w[k]/lu
        y1 <- (1-S[k]) * w[k]/lu
        y2 <- w[k]/lu
        i <- floor(x)
        j1 <- floor(y1)
        j2 <- floor(y2)
        weights <- c((i+1-x)*(j1+1-y1), (i+1-x)*(y1-j1), (x-i)*(y1-j1), (j1+1-y1)*(x-i))
        dpsplit <- dp[k] * weights
        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)]
        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)
}

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.integer(N), q = double(N))$q
}

lossdistribC.truncated <- 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.integer(N), q = double(N))$q
}

lossdistribC.joint <- function(p, w, S, N){
    ## C version of lossdistrib.joint, roughly 20 times faster
    dyn.load("lossdistrib.dll")
    .C("lossdistrib_joint", as.double(p), as.integer(length(p)), as.double(w),
       as.double(S), as.integer(N), q = matrix(0, N, N))$q
}

lossdistribprepayC.joint <- function(dp, pp, w, S, N){
    ## C version of lossdistribprepay.joint
    dyn.load("lossdistrib.dll")
    .C("lossdistrib_prepay_joint", as.double(dp), as.double(pp), as.integer(length(dp)),
       as.double(w), as.double(S), as.integer(N), q=matrix(0, N, N))$q
}

lossrecovdist <- function(defaultprob, prepayprob, w, S, N, useC=TRUE){
    if(all(!prepayprob)){
        if(useC){
            L <- lossdistribC(defaultprob, w, S, N)
            R <- lossdistribC(defaultprob, w, 1-S, N)
        }else{
            L <- lossdistrib2(defaultprob, w, S, N)
            R <- lossdistrib2(defaultprob, w, 1-S, N)
        }
    }else{
        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, N, 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(all(!prepayprob)){
        for(t in 1:ncol(defaultprob)){
            temp <- lossrecovdist(defaultprob[,t], 0, w, S[,t], N, 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, 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.pl <- function(L, R, 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.pv <- function(L, R, cs, K1, K2, Ngrid=nrow(L)){
    return( tranche.pl(L, R, cs, K1, K2, Ngrid) + tranche.cl(L, R, cs, K1, K2, Ngrid))
}

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, T=N, n.int=100){
    dyn.load("lossdistrib.dll")
    quadrature <- gauss.quad.prob(n.int, "normal")
    Z <- quadrature$nodes
    w <- quadrature$weights
    L <- matrix(0, T, dim(SurvProb)[2])
    R <- matrix(0, T, 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.integer(T), 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, dist2$R, cs, 0, K[2])
    pl1 <- tranche.pl(dist1$L, dist1$R, 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) - indexpv(portfoliominus, index)
    if(K2==1){
        dPVtranche <- BCtranche.pv(portfolioplus, index, coupon, 0, K1, 0, rho1, lu)$bp -
            BCtranche.pv(portfoliominus, index, coupon, 0, K1, 0, rho1, lu)$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, lu)$bp -
            BCtranche.pv(portfoliominus, index, coupon, K1, K2, rho1, rho2, lu)$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, program$weight, rho, defaultprob[i,j])
        }
    }
    return( p )
}