## 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

GHquad <- function(n){
    n <- as.integer(n)
    Z <- double(n)
    w <- double(n)
    result <- .C("GHquad", n, Z=Z, w=w)
    result[[1]] <- NULL
    return(result)
}

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

lossdistCblas <- function(p, w, S, N, defaultflag=FALSE){
    ## C version of lossdistrib2, roughly 50 times faster
    .C("lossdistrib_blas", as.double(p), as.integer(length(p)),
       as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), q = double(N))$q
}

lossdistCZ <- function(p, w, S, N, defaultflag=FALSE, rho, Z, wZ){
    .C("lossdistrib_Z", as.double(p), as.integer(length(p)),
       as.double(w), as.double(S), as.integer(N), as.logical(defaultflag),
       as.double(rho), as.double(Z), as.integer(length(Z)),
       q = matrix(0, N, length(Z)))$q
}

lossdistC.truncated <- function(p, w, S, N, T=N){
    ## C version of lossdistrib2, roughly 50 times faster
    .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
    .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
    .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.jointblas <- function(p, w, S, N, defaultflag=FALSE){
    ## C version of lossdistrib.joint, roughly 20 times faster
    .C("lossdistrib_joint_blas", 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.jointZ <- function(dp, w, S, N, defaultflag = FALSE, rho, Z, wZ){
    ## N is the size of the grid
    ## dp is of size n.credits
    ## w is of size n.credits
    ## S is of size n.credits by nZ
    ## rho is a double
    ## Z is a vector of length nZ
    ## w is a vector if length wZ
    r <- .C("lossdistrib_joint_Z", as.double(dp), as.integer(length(dp)),
            as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), as.double(rho),
            as.double(Z), as.double(wZ), as.integer(length(Z)), q = matrix(0, N, N))$q
}

lossdistC.prepay.joint <- function(dp, pp, w, S, N, defaultflag=FALSE){
    ## C version of lossdist.prepay.joint
    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
    return(r)
}

lossdistC.prepay.jointZ <- function(dp, pp, w, S, N, defaultflag = FALSE, rho, Z, wZ){
  ## N is the size of the grid
  ## dp is of size n.credits
  ## pp is of size n.credits
  ## w is of size n.credits
  ## S is of size n.credits by nZ
  ## rho is a double
  ## Z is a vector of length nZ
  ## w is a vector if length wZ
    r <- .C("lossdistrib_prepay_joint_Z", as.double(dp), as.double(pp), as.integer(length(dp)),
            as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), as.double(rho),
            as.double(Z), as.double(wZ), as.integer(length(Z)), q = matrix(0, N, N))$q
}

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]
    for(t in 1:dim(dist.joint)[1]){
        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]
                }else{
                    distDR[t,Ngrid,j] <- distDR[t,Ngrid,j] +
                        dist.joint[t,i,j]
                }
            }
        }
        distDR[t,,] <- distDR[t,,]/sum(distDR[t,,])
    }
    return( distDR )
}

shockprob <- function(p, rho, Z, log.p=F){
  ## computes the shocked default probability as a function of the copula factor
  ## function is vectorized provided the below caveats:
  ## p and rho are vectors of same length n, Z is a scalar, returns vector of length n
  ## p and rho are scalars, Z is a vector of length n, returns vector of length n
  if(length(p)==1){
    if(rho==1){
      return(Z<=qnorm(p))
    }else{
      return(pnorm((qnorm(p)-sqrt(rho)*Z)/sqrt(1-rho), log.p=log.p))
    }
  }else{
    result <- double(length(p))
    result[rho==1] <- Z<=qnorm(p[rho==1])
    result[rho<1] <- pnorm((qnorm(p[rho<1])-sqrt(rho[rho<1])*Z)/sqrt(1-rho[rho<1]), log.p=log.p)
    return( result )
  }
}

shockseverity <- function(S, Stilde=1, Z, rho, p){
  ## computes the severity as a function of the copula factor Z
  result <- double(length(S))
  result[p==0] <- 0
  result[p!=0] <- Stilde * exp( shockprob(S[p!=0]/Stilde*p[p!=0], rho[p!=0], Z, TRUE) -
           shockprob(p[p!=0], rho[p!=0], Z, TRUE))
  return(result)
}

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
    require(distr)
    if(p0==0){
        return(0)
    }
    if(rho == 1){
        distw <- DiscreteDistribution(Z, w)
        return(pnorm(q(distw)(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)
}

fit.probC <- function(Z, w, rho, p0){
    r <- .C("fitprob", as.double(Z), as.double(w), as.integer(length(Z)),
            as.double(rho), as.double(p0), q = double(1))
    return(r$q)
}

stochasticrecov <- function(R, Rtilde, Z, w, rho, porig, pmod){
    ## if porig == 0 (probably matured asset) then return orginal recovery
    ## it shouldn't matter anyway since we never default that asset
    if(porig == 0){
        return(rep(R, length(Z)))
    }else{
        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))))
    }
}

stochasticrecovC <- function(R, Rtilde, Z, w, rho, porig, pmod){
    r <- .C("stochasticrecov", as.double(R), as.double(Rtilde), as.double(Z),
            as.double(w), as.integer(length(Z)), as.double(rho), as.double(porig),
            as.double(pmod), q = double(length(Z)))
    return(r$q)
}

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

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(defaultprob, issuerweights, recov, rho, Z, w,
                       N=length(recov)+1, defaultflag=FALSE, n.int=500){
  if(missing(Z)){
    quadrature <- GHquad(n.int)
    Z <- quadrature$Z
    w <- quadrature$w
  }
  ## do not use if weights are not gaussian, results would be incorrect
  ## since shockseverity is invalid in that case (need to use stochasticrecov)
  LZ <- matrix(0, N, length(Z))
  RZ <- matrix(0, N, length(Z))
  L <- matrix(0, N, ncol(defaultprob))
  R <- matrix(0, N, ncol(defaultprob))
  for(t in 1:ncol(defaultprob)){
    for(i in 1:length(Z)){
      g.shocked <- shockprob(defaultprob[,t], rho, Z[i])
      S.shocked <- shockseverity(1-recov, 1, Z[i], rho, defaultprob[,t])
      temp <- lossrecovdist(g.shocked, , 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(defaultprob, issuerweights, recov, rho, Z, w,
                        N=length(issuerweights)+1, defaultflag=FALSE){
    L <- matrix(0, N, dim(defaultprob)[2])
    R <- matrix(0, N, dim(defaultprob)[2])
    rho <- rep(rho, length(issuerweights))
    r <- .C("BCloss_recov_dist", defaultprob, dim(defaultprob)[1], dim(defaultprob)[2],
            as.double(issuerweights), as.double(recov), as.double(Z), as.double(w),
            as.integer(length(Z)), as.double(rho), as.integer(N), as.logical(defaultflag), L=L, R=R)
    return(list(L=r$L,R=r$R))
}

BCtranche.legs <- function(index, K, rho){
    ## computes the protection leg and couponleg of a 0-K tranche
    dist <- BClossdistC(index$defaultprob, index$issuerweights, index$recov, rho, index$Z, index$w, index$N)
    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){
    ## computes the protection leg, couponleg, and bond price of a tranche
    ## in the base correlation setting

    pl <- rep(0, length(index$rho)-1)
    cl <- rep(0, length(index$rho)-1)
    for(i in 2:(length(index$rho)-1)){
        temp <- BCtranche.legs(index, index$K[i], index$rho[i])
        pl[i] <- temp$pl
        cl[i] <- temp$cl
    }
    dK <- diff(index$K)
    plvec <- diff(pl)/dK[-length(dK)]
    clvec <- diff(cl)/dK[-length(dK)]*index$tranche.running[-length(index$tranche.running)]
    ## compute the supersenior as the left over
    temp <- indexpv(index, tradedate = tradedate, clean = FALSE)
    indexpl <- temp$pl
    indexcl <- temp$cl/index$quotes$spread
    plvec <- c(plvec, -(indexpl+pl[4])/dK[length(dK)])
    clvec <- c(clvec, (indexcl-cl[4])*index$tranche.running[length(index$tranche.running)]/dK[length(dK)])
    if(protection){
        bp <- -plvec-clvec
    }else{
        bp <- 1+plvec+clvec
    }
    return(list(pl=plvec, cl=clvec, bp=bp))
}

adjust.skew <- function(index1, index2, method="ATM"){
    #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

    K1 <- index2$K[-c(1,length(index1$K))]
    K2 <- index2$K[-c(1,length(index1$K))]
    aux <- function(x, index1, el1, skew, index2, el2, K2){
        return(abs(ELtrunc(index1, x, skew(x))/el1-
                   ELtrunc(index2, K2, skew(x))/el2))
    }
    aux2 <- function(x, index1, skew, index2, K2){
        return(abs(log(Probtrunc(index1, x, skew(x)))-
                   log(Probtrunc(index2, K2, skew(x)))))
    }
    if(method %in% c("ATM", "TLP")){
        el1 <- EL(index1)
        el2 <- EL(index2)
    }
    skew <- function(x){
        #we cap the correlation at 0.99
        f <- splinefun(K1, index1$rho[-c(1, length(index1$rho))], "natural")
        return(min(f(x), 0.99))
    }
    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.3
        for(K2val in K2){
            prog <- optimize(aux2, interval=c(0, m),
                             index1=index1, skew=skew, index2=index2, K2=K2val)
            K1eq <- c(K1eq, prog$minimum)
        }
     }
    return(c(0, skew(K1eq), NA))
}

theta.adjust.skew <- function(index, shortened=4, method="ATM"){
    #ajust the correlation skew by doing ATM mapping on the expected loss
    indexshort <- index
    N <- nrow(index$cs)-shortened
    indexshort$defaultprob <- indexshort$defaultprob[,1:N]
    indexshort$cs <- indexshort$cs[1:N,]
    return(adjust.skew(index, indexshort, method))
}

BCtranche.theta <- function(index, tradedate=Sys.Date(), shortened=4){
    temp <- BCtranche.pv(index)
    rho.adj <- theta.adjust.skew(index, shortened)
    if(any(rho.adj[-c(1, length(rho.adj))]<=0)){
        print("probable inverted skew: no adjustment")
    }else{
        index$rho <- rho.adj
    }
    N <- nrow(index$cs)-shortened
    index$cs <- index$cs[1:N,]
    index$defaultprob <- index$defaultprob[,1:N]
    temp2 <- BCtranche.pv(index)
    temp3 <- BCtranche.delta(index, tradedate)
    return(list(theta=temp2$bp-temp$bp+index$tranche.running,
                delta=temp3$delta))
}

BCtranche.delta <- function(index, tradedate = Sys.Date()){
    ## 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
    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)){
        temp <- BCtranche.pv(index.list[[j]])
        indexbp[j] <- indexpv(index.list[[j]], tradedate = tradedate, clean = FALSE)$bp
        bp[,j] <- temp$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( list(deltas=deltas, gammas=gammas) )
}

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

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)-shortened
    ELvec <- as.numeric(crossprod(index$issuerweights * (1-index$recov), index$defaultprob[,1:Ncol]))
    if(!discounted){
       return( diff(c(0, ELvec)) )
   }else{
       return( sum(index$cs$df[1:Ncol]*diff(c(0, ELvec))) )
   }
}

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(defaultprob, issuerweights, recov, cs, K, rho, Z, w, N=101) {
    ## 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){
    ## 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[i], defaultprob[i,j])
        }
    }
    return( p )
}

MFupdate.probC <- function(Z, w, rho, defaultprob){
    ## update the probabilities 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] <- .C("fitprob", as.double(Z), as.double(w), as.integer(length(Z)),
                     as.double(rho[i]), as.double(defaultprob[i,j]), q = double(1))$q
      }
    }
    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 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(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 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.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 )
}

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