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