library(lossdistrib)
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) )
        }
    }
}

funded.tranche.pv <- function(L, R, cs, K1, K2, scaled = FALSE){
    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)])))
        interest <- crossprod(sizeadj * cs$coupons, cs$df)
        principal <- diff(c(0, trancherecov(R, K1, K2)))
        principal[length(principal)] <- principal[length(principal)] + size[length(size)]
        principal <- crossprod(cs$df, principal)
        if(scaled) {
            pv <- (interest + principal) / (K2-K1)
        }else{
            pv <- (interest + principal)
        }
        return(pv)
    }
}

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) )
}

BCtranche.legs <- function(index, K, rho, complement=FALSE) {
    ## computes the protection leg and couponleg of a 0-K tranche
    ## if complement==TRUE, computes the protection leg and coupon leg of a K-1 tranche
    if((K==0 && !complement) || (K==1 && complement)) {
        return(list(cl=0, pl=0))
    } else if((K==1 && !complement) || (K==0 && complement)) {
        return(BCindex.pv(index))
    } else {
        dist <- BClossdistC(index$defaultprob, index$issuerweights,
                            index$recov, rho, index$Z, index$w, index$N)
        if(complement) {
            return(list(cl=tranche.cl(dist$L, dist$R, index$cs, K, 1),
                        pl=tranche.pl(dist$L, index$cs, K, 1)))
        } else {
            return(list(cl=tranche.cl(dist$L, dist$R, index$cs, 0, K),
                        pl=tranche.pl(dist$L, index$cs, 0, K)))
        }
    }
}

BCtranche.pv <- function(index, protection=FALSE, complement=FALSE){
    ## computes the protection leg, couponleg, and bond price of a tranche
    ## in the base correlation setting
    ## if complement=FALSE compute the pvs starting from 0 (bottom-up skew)
    ## if complement=TRUE compute the pvs starting from 1 (top-down skew)
    pl <- rep(0, length(index$rho))
    cl <- rep(0, length(index$rho))
    for(i in seq_along(index$rho)) {
        temp <- BCtranche.legs(index, index$K[i], index$rho[i], complement)
        pl[i] <- temp$pl
        cl[i] <- temp$cl
    }
    dK <- diff(index$K)
    plvec <- diff(pl)/dK
    clvec <- diff(cl)/dK*index$tranches$running
    if(complement) {
        plvec <- -plvec
        clvec <- -clvec
    }
    if(protection) {
        bp <- - plvec - clvec
    }else{
        bp <- 1 + plvec + clvec
    }
    return(list(pl=plvec, cl=clvec, bp=bp))
}

MFtranche.pv <- function(index, dist, protection=FALSE) {
    n.tranches <- length(index$K) - 1
    pl <- rep(0, n.tranches)
    cl <- rep(0, n.tranches)
    for(i in seq.int(n.tranches)) {
        pl[i] <- tranche.pl(dist$L, index$cs, index$K[i], index$K[i+1])
        cl[i] <- tranche.cl(dist$L, dist$R, index$cs, index$K[i], index$K[i+1])
    }
    dK <- diff(index$K)
    plvec <- pl/dK
    clvec <- cl/dK*index$tranches$running
    if(protection) {
        bp <- -plvec-clvec
    } else {
        bp <- 1+plvec+clvec
    }
    return(list(pl=plvec, cl=clvec, bp=bp))
}

adjust.skew <- function(index1, index2, method=c("ATM", "TLP", "PM")) {
    #index1 is the index for which we already have computed the skew
    #index2 is the index we're mapping to
    # if method="ATM", do simple at the money mapping
    #    method="TLP", do tranche loss proportion mapping
    #    method="PM", do probability matching
    method <- match.arg(method)
    K1 <- index1$K[-c(1,length(index1$K))]
    K2 <- index2$K[-c(1,length(index2$K))]

    aux <- function(x, index1, el1, skew, index2, el2, K2) {
        newrho <- cap(skew(x))
        return(abs(ELtrunc(index1, x, newrho)/el1-
                   ELtrunc(index2, K2, newrho)/el2))
    }

    aux2 <- function(x, index1, skew, index2, K2) {
        newrho <- cap(skew(x))
        return(abs(log(Probtrunc(index1, x, newrho))-
                   log(Probtrunc(index2, K2, newrho))))
    }

    if(method %in% c("ATM", "TLP")) {
        el1 <- EL(index1)
        el2 <- EL(index2)
    }

    skew <- splinefun(K1, index1$rho[-c(1, length(index1$rho))], "natural")
    cap <- function(x) {
        pmax(pmin(x, 0.99), 0.01)
    }

    if(method == "ATM") {
        K1eq <- el1/el2 * K2
    } else if(method == "TLP") {
        K1eq <- c()
        m <- max(K2) + 0.3
        for(K2val in K2) {
            prog <- optimize(aux, interval=c(0,m),
                             index1=index1, el1=el1, skew=skew,
                             index2=index2, el2=el2, K2=K2val)
            K1eq <- c(K1eq, prog$minimum)
        }
    } else if (method == "PM") {
        K1eq <- c()
        m <- max(K2) + 0.25
        for(K2val in K2) {
            prog <- optimize(aux2, interval=c(K2val*0.1, K2val*1.8),
                             index1=index1, skew=skew, index2=index2, K2=K2val)
            K1eq <- c(K1eq, prog$minimum)
        }
     }
    return(c(NA, cap(skew(K1eq)), NA))
}

BCtranche.theta <- function(index, shortened=4, complement=FALSE,
                            method=c("ATM", "TLP", "PM")) {
    N <- nrow(index$cs) - shortened
    if(N < 0) {
        stop("Maturity too short for computing tranche thetas")
    }
    indexshort <- index
    indexshort$defaultprob <- index$defaultprob[,1:N]
    indexshort$cs <- index$cs[1:N,]
    indexshort$rho <- adjust.skew(index, indexshort, method)
    temp <- BCtranche.pv(index, complement=complement)
    temp2 <- BCtranche.pv(indexshort, complement=complement)
    temp3 <- BCtranche.delta(indexshort, complement=complement)
    return(data.frame(theta=temp2$bp-temp$bp+index$tranches$running,
                      fw.delta=temp3$delta))
}

BCtranche.delta <- function(index, complement=FALSE) {
    ## computes the tranche delta (on current notional) by doing a proportional
    ## blip of all the curves
    ## if complement is False, then computes deltas bottom-up
    ## if complement is True, then computes deltas top-down
    eps <- 1e-4
    index$N <- 301 ## for gamma computations we need all the precision we can get
    ## we build a lit of 4 indices with various shocks
    index.list <- lapply(c(0, eps, -eps, 2*eps), function(x) {
        if(x==0) {
            return(index)
        } else {
            newindex <- index
            newindex$portfolio <- tweakportfolio(newindex$portfolio, x)
            newindex$defaultprob <- 1 - SPmatrix(newindex$portfolio, index$cs$dates)
            return(newindex)
        }
    })
    bp <- matrix(0, length(index$K)-1, length(index.list))
    indexbp <- rep(0, length(index.list))
    for(j in seq_along(index.list)) {
        indexbp[j] <- BCindex.pv(index.list[[j]])$bp
        bp[,j] <- BCtranche.pv(index.list[[j]], complement=complement)$bp
    }

    deltas <- (bp[,2]-bp[,3])/(indexbp[2]-indexbp[3])*tranche.factor(index)/index$factor
    deltasplus <- (bp[,4]-bp[,1])/(indexbp[4]-indexbp[1])*tranche.factor(index)/index$factor
    gammas <- (deltasplus-deltas)/(indexbp[2]-indexbp[1])/100

    return( data.frame(delta=deltas, gamma=gammas) )
}

MFtranche.delta <- function(index){
    ## computes the tranche delta (on current notional) by doing a proportional
    ## blip of all the curves
    ## if complement is False, then computes deltas bottom-up
    ## if complement is True, then computes deltas top-down
    eps <- 1e-4
    index$Ngrid <- 301 ## for gamma computations we need all the precision we can get
    ## we build a lit of 4 indices with various shocks
    index.list <- lapply(c(0, eps, -eps, 2*eps), function(x) {
        if(x==0) {
            return(index)
        } else {
            newindex <- index
            newindex$portfolio <- tweakportfolio(newindex$portfolio, x)
            newindex$defaultprob <- 1 - SPmatrix(newindex$portfolio, length(index$cs$dates))
            return(newindex)
        }
    })
    bp <- matrix(0, length(index$K)-1, length(index.list))
    indexbp <- rep(0, length(index.list))
    for(j in seq_along(index.list)) {
        indexbp[j] <- BCindex.pv(index.list[[j]])$bp
        dist <- MFdist(index.list[[j]])
        bp[,j] <- BCtranche.pv(index.list[[j]], dist)
    }

    deltas <- (bp[,2]-bp[,3])/(indexbp[2]-indexbp[3])*tranche.factor(index)/index$factor
    deltasplus <- (bp[,4]-bp[,1])/(indexbp[4]-indexbp[1])*tranche.factor(index)/index$factor
    gammas <- (deltasplus-deltas)/(indexbp[2]-indexbp[1])/100

    return( list(deltas=deltas, gammas=gammas) )
}

BCtranche.corr01 <- function(index, eps=0.01, complement=FALSE) {
    ##does a parallel shift of the skew and computes the change in pv
    before <- BCtranche.pv(index, complement=complement)
    index$rho[-1] <- index$rho[-1]^(1-eps)
    after <- BCtranche.pv(index, complement=complement)
    return(after$bp-before$bp)
}

BCtranche.VaR <- function(index, skew.shocks, complement=FALSE) {
    before <- BCtranche.pv(index, complement=complement)
    n <- length(index$rho)
    r <- matrix(0, nrow(skew.shocks), n-1)
    logskew <- log(index$rho[-c(1,n)])
    for(i in 1:nrow(skew.shocks)) {
        index$rho[-c(1,n)] <- exp(logskew * (1 + skew.shocks[i,]))
        after <- BCtranche.pv(index, complement=complement)
        r[i,] <- before$bp - after$bp
    }
    return( r )
}

EL <- function(index, discounted=TRUE, shortened=0) {
    ## computes the expected loss of a portfolio (time discounted if discounted is TRUE)
    ## given the default curves and recovery
    ## should be very close to the protection leg of the portfolio of cds
    ## index should be a list with issuerweights, recov, defaultprob and cs parameters
    ## shortened: number of quarters to shorten the maturity by
    Ncol <- ncol(index$defaultprob)
    if(is.null(Ncol) && shortened==0) {
        DP <- index$defaultprob
        df <- index$cs$df
    }else{
        DP <- index$defaultprob[,1:(Ncol-shortened)]
        df <- index$cs$df[1:(Ncol-shortened)]
    }
    ELvec <- as.numeric(crossprod(index$issuerweights * (1-index$recov), DP))
    if(!discounted) {
       return( ELvec[length(ELvec)] )
    } else {
        return( sum(df*diff(c(0, ELvec))) )
   }
}

BCindex.pv <- function(index, discounted=TRUE, shortened=0) {
    Ncol <- ncol(index$defaultprob)
    if(is.null(Ncol) && shortened == 0) {
        DP <- index$defaultprob
        Ncol <- 1
    } else {
        DP <- index$defaultprob[,1:(Ncol-shortened)]
    }
    ELvec <- as.numeric(crossprod(index$issuerweights * (1-index$recov), DP))
    size <- 1 - as.numeric(crossprod(index$issuerweights, DP))
    sizeadj <- 0.5*(c(1, size[-length(size)])+size)
    if(!discounted) {
        pl <- -ELvec[length(ELvec)]
        cl <- as.numeric(crossprod(index$cs$coupons[1:Ncol], sizeadj))
    } else {
        pl <-  -sum(index$cs$df[1:Ncol]* diff(c(0, ELvec)))
        cl <- as.numeric(crossprod(index$cs$coupons[1:Ncol], sizeadj * index$cs$df[1:Ncol] ))
    }
    bp <- 1 + cl*index$quotes$spread+pl
    return(list(pl=pl, cl=cl, bp=bp))
}

ELtrunc <- function(index, K, rho) {
    ## computes the expected loss of a portfolio below strike K
    ## could be written faster by using a truncated version of lossdist
    ## index should be a list with issuerweights, recov, defaultprob and cs parameters
    dist <- BClossdistC(index$defaultprob, index$issuerweights,
                        index$recov, rho, index$Z, index$w, index$N)
    return( -tranche.pl(dist$L, index$cs, 0, K))
}

Probtrunc <- function(index, K, rho) {
    dist <- BClossdistC(index$defaultprob, index$issuerweights, index$recov,
                        rho, index$Z, index$w, index$N)
    p <- cumsum(dist$L[,ncol(dist$L)])
    support <- seq(0, 1, length=index$N)
    probfun <- splinefun(support, p, method="hyman")
    return(probfun(K))
}


BCstrikes <- function(index, K, rho) {
    ## computes the strikes as a percentage of expected loss
    ## Kmodified is the current attachment points (adjusted for losses)
    el <- EL(index)
    ELvec <- c()
    for(i in 2:length(K)){
        ELvec <- c(ELvec, ELtrunc(index, K[i], rho[i]))
    }
    return( ELvec/el )
}

tranche.factor <- function(index){
    ## compute the factor to convert from delta on current notional to delta on original notional
    ## K1 and K2 original strikes
    return( diff(index$K)/diff(index$K.orig)*index$factor )
}

MFupdate.prob <- function(Z, w, rho, defaultprob, useC = TRUE){
    ## 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))
    fit.prob <- if(useC) fit.probC else fit.prob
    for(i in 1:nrow(defaultprob)) {
        for(j in 1:ncol(defaultprob)) {
            p[i,j] <- fit.prob(Z, w, rho[i], defaultprob[i,j])
        }
    }
    return( p )
}

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 default/loss distribution depending if
    ## defaultflag is TRUE or FALSE (default setting is FALSE)
    ## 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 axis 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(w, Z, rho, defaultprob, defaultprobmod,
                                    prepayprob, prepayprobmod, issuerweights, recov,
                                    Ngrid=2*length(issuerweights)+1, defaultflag=FALSE) {
  ## 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 axis 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.credit, n.int, ncol(defaultprob)))

  for(t in 1:ncol(defaultprob)) {
    for(i in 1:n.credit){
      Rstoch[i,,t] <- stochasticrecovC(recov[i], 0, Z, w, rho[i],
                                       defaultprob[i,t], defaultprobmod[i,t])
    }
  }
  Q <- array(0, dim=c(ncol(defaultprob), Ngrid, Ngrid))
  for(t in 1:ncol(defaultprob)) {
      S <- 1 - Rstoch[,,t]
      Q[t,,] <- lossdistC.prepay.jointZ(defaultprobmod[,t], prepayprobmod[,t],
                                        issuerweights, S, Ngrid, defaultflag, rho, Z, w)
  }
  return( Q )
}

MFrecovery <- function(index, defaultprobmod) {
    n.credit <- length(index$issuerweights)
    n.int <- length(index$Z)
    Rstoch <- array(0, dim=c(n.credit, n.int, ncol(index$defaultprob)))
    rho <- rep(0.45, n.credit)
    for(t in 1:ncol(index$defaultprob)) {
        for(i in 1:n.credit) {
            Rstoch[i,,t] <- stochasticrecovC(index$recov[i], 0, index$Z, index$w.mod,
                                             rho[i], index$defaultprob[i,t], defaultprobmod[i,t])
        }
    }
    return( Rstoch )
}

MFlossdist <- function(index) {
    n.credit <- length(index$issuerweights)
    rho <- rep(0.45, n.credit)
    defaultprobmod <- MFupdate.prob(index$Z, index$w.mod, rho, index$defaultprob)
    n.credit <- length(index$issuerweights)
    n.int <- length(index$Z)
    Rstoch <- MFrecovery(index, defaultprobmod)
    Lw <- matrix(0, index$N, n.int)
    Rw <- matrix(0, index$N, n.int)
    L <- matrix(0, index$N, ncol(index$defaultprob))
    R <- matrix(0, index$N, ncol(index$defaultprob))
    for(t in 1:ncol(index$defaultprob)) {
        S <- 1 - Rstoch[,,t]
        Lw <- lossdistCZ(defaultprobmod[,t], index$issuerweights, S, index$N, 0, rho, index$Z)
        Rw <- lossdistCZ(defaultprobmod[,t], index$issuerweights, 1-S, Ngrid, 0, rho, index$Z)
        L[,t] <- Lw%*%index$w.mod
        R[,t] <- Rw%*%index$w.mod
    }
    return(list(L=L, R=R))
}