fix up the fft computation

1 files changed, 629 insertions, 610 deletions

diff --git a/R/distrib.R b/R/distrib.R
index edc5081..8eb8c42 100644
--- a/R/distrib.R
+++ b/R/distrib.R
@@ -1,610 +1,629 @@
-## 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

-

-#' Gauss-Hermite quadrature weights

-#'

-#' \code{GHquad} computes the quadrature weights for integrating against a

-#' Gaussian distribution.

-#'

-#' if f is a function, then \eqn{\sum_{i=1}^n f(Z_i)w_i \approx

-#' \frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty f(x)e^{-\frac{x^2}{2}}\,dx}.

-#' @param n Integer, the number of nodes

-#' @return A list with two components:

-#'  \item{Z}{the list of nodes}

-#'  \item{w}{the corresponding weights}

-#'

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

-}

-

-#' Loss distribution of a portfolio

-#'

-#' \code{lossdistrib} computes the probability distribution of a sum

-#' of independent Bernouilli variables with unequal probabilities.

-#'

-#' We compute the probability distribution of \eqn{S = \sum_{i=1}^n X_i}

-#' where \eqn{X_i} is Bernouilli(\eqn{p_i}). This uses the simple recursive

-#' algorithm of Andersen, Sidenius and Basu

-#' @param p Numeric vector, the vector of success probabilities

-#' @return A vector q such that \eqn{q_k=\Pr(S=k)}

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

-}

-

-#' Loss distribution of a portfolio

-#'

-#' \code{lossdistrib.fft} computes the probability distribution of a sum

-#' of independent Bernouilli variables with unequal probabilities.

-#'

-#' We compute the probability distribution of \eqn{S = \sum_{i=1}^n X_i}

-#' where \eqn{X_i} is Bernouilli(\eqn{p_i}).

-#' This uses the FFT, thus omplexity is of order \eqn{O(n m) + O(m\log(m))}

-#' where \eqn{m} is the size of the grid and \eqn{n}, the number of probabilities.

-#' It is slower than the recursive algorithm in practice.

-#' @param p Numeric vector, the vector of success probabilities

-#' @return A vector such that \eqn{q_k=\Pr(S=k)}

-lossdistrib.fft <- function(p){

-    n <- length(p)

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

-}

-

-#' Loss distribution of a portfolio

-#'

-#' \code{lossdistrib2} computes the probability distribution of a sum

-#' of independent Bernouilli variables with unequal probabilities.

-#'

-#' We compute the probability distribution of \eqn{L = \sum_{i=1}^n w_i S_i X_i}

-#' where \eqn{X_i} is Bernouilli(\eqn{p_i}). If \code{defaultflag} is TRUE, we

-#' compute the distribution of \eqn{D = \sum_{i=1}^n w_i X_i} instead.

-#' This a recursive algorithm with first order correction for discretization.

-#' @param p Numeric, vector of default probabilities

-#' @param w Numeric, vector of weights

-#' @param S Numeric, vector of severities

-#' @param N Integer, number of ticks in the grid

-#' @param defaultflag Boolean, if TRUE, we compute the default distribution

-#' (instead of the loss distribution).

-#' @return a Numeric vector of size \code{N} computing the loss (resp.

-#' default) distribution if \code{defaultflag} is FALSE (resp. TRUE).

-lossdistrib2 <- function(p, w, S, N, defaultflag=FALSE){

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

-}

-

-#' Loss distribution truncated version

-#'

-#' \code{lossdistrib2.truncated} computes the probability distribution of a sum

-#' of independent Bernouilli variables with unequal probabilities up

-#' to a cutoff N.

-#'

-#' This is actually slower than \code{lossdistrib2}, but in C this is

-#' twice as fast. For high severities, M can become bigger than the cutoff, and

-#' there is some probability mass escaping.

-#' @param p Numeric, vector of default probabilities

-#' @param w Numeric, vector of weights

-#' @param S Numeric, vector of severities

-#' @param N Integer, number of ticks in the grid

-#' @param cutoff Integer, where to stop computing the exact probabilities

-#' @return a Numeric vector of size \code{N} computing the loss distribution

-lossdistrib2.truncated <- function(p, w, S, N, cutoff=N){

-    n <- length(p)

-    lu <- 1/(N-1)

-    q <- rep(0, cutoff)

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

-}

-

-#' Recovery distribution of a portfolio

-#'

-#' \code{recovdist} computes the recovery distribution of portfolio

-#' described by a vector of default probabilities, and prepay probabilities.

-#' \eqn{R=\sum_{i=1}^n w_i X_i} where \eqn{X_i=0} w.p. \eqn{1-dp_i-pp_i},

-#' \eqn{X_i=1-S_i} with probability \eqn{dp_i}, and \eqn{X_i=1} w.p. \eqn{pp_i}

-#'

-#' It is a recursive algorithm with first-order correction. For a unit of loss

-#' \eqn{lu}, each non-zero value \eqn{v} is interpolated on the grid

-#' as the pair of values

-#' \eqn{\left\lfloor\frac{v}{lu}\right\rfloor}  and

-#' \eqn{\left\lceil\frac{v}{lu}\right\rceil} so that \eqn{X_i} has

-#' four non zero values.

-#' @param dp Numeric, vector of default probabilities

-#' @param pp Numeric, vector of prepay probabilities

-#' @param w Numeric, vector of weights

-#' @param S Numeric, vector of severities

-#' @param N Integer, number of ticks in the grid

-#' @return a Numeric vector of size \code{N} computing the recovery distribution

-recovdist <- function(dp, pp, w, S, N){

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

-}

-

-#' Joint loss-recovery distributionrecursive algorithm with first order correction to compute the joint

-#' probability distribution of the loss and recovery.

-#'

-#' For high severities, M can become bigger than N, and there is

-#' some probability mass escaping.

-#' @param p Numeric, vector of default probabilities

-#' @param w Numeric, vector of weights

-#' @param S Numeric, vector of severities

-#' @param N Integer, number of ticks in the grid

-#' @param defaultflab Logical, whether to return the loss or default distribution

-#' @return q Matrix of joint loss, recovery probability distribution

-#' colSums(q) is the recovery distribution marginal

-#' rowSums(q) is the loss distribution marginal

-lossdist.joint <- function(p, w, S, N, defaultflag=FALSE){

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

-}

-

-lossdistCZ <- function(p, w, S, N, defaultflag=FALSE, rho, Z){

-    ##S is of size (length(p), length(Z))

-    stopifnot(length(rho)==length(p),

-              length(rho)==length(w),

-              nrow(S)==length(p),

-              ncol(S)==length(Z))

-    .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, defaultflag=FALSE){

-    ## truncated version of lossdistrib

-    ## q[i] is 0 for i>=T

-    .C("lossdistrib_truncated", as.double(p), as.integer(length(p)),

-       as.double(w), as.double(S), as.integer(N), as.integer(T), as.logical(defaultflag),

-       q = double(N))$q

-}

-

-exp.trunc <- function(p, w, S, N, K){

-    ## computes E[(K-L)^+]

-    r <- 0

-    .C("exp_trunc", as.double(p), as.integer(length(p)),

-       as.double(w), as.double(S), as.integer(N), as.double(K), res = r)$res

-}

-

-rec.trunc <- function(p, w, S, N, K){

-    ## computes E[(K-(1-R))^+] = E[(\tilde K- \bar R)]

-    ## where \tilde K = K-sum_i w_i S_i and \bar R=\sum_i w_i R_i (1-X_i)

-    Ktilde <- K-crossprod(w, S)

-    if(Ktilde < 0){

-        return( 0 )

-    }else{

-        return( exp.trunc(1-p, w, 1-S, N, Ktilde) )

-    }

-}

-

-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.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 vector of doubles of size n.credits

-    ## Z is a vector of length nZ

-    ## w is a vector if length wZ

-

-    r <- .C("lossdistrib_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)), output = matrix(0,N,N))

-    return(r$output)

-}

-

-lossrecovdist <- function(defaultprob, prepayprob, w, S, N, defaultflag=FALSE, useC=TRUE){

-    lossdistrib2 <- if(useC) lossdistC

-    recovdist <- if(useC) recovdistC

-    if(missing(prepayprob)){

-        L <- lossdistrib2(defaultprob, w, S, N, defaultflag)

-        R <- lossdistrib2(defaultprob, w, 1-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))

-    lossdist.joint <- if(useC) lossdistC.jointblas

-    lossdist.prepay.joint <- if(useC) lossdistC.prepay.jointblas

-    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.jointblas(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag)

-        }

-    }

-    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

-    if(p0==0){

-        return(0)

-    }

-    if(rho == 1){

-        distw <- distr::DiscreteDistribution(Z, w)

-        return(distr::pnorm(distr::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){

-    stopifnot(length(Z)==length(w))

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

-}

-

-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

-    }

-    stopifnot(length(Z)==length(w),

-              nrow(defaultprob) == length(issuerweights))

-

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

-    if(is.null(dim(defaultprob))){

-        dim(defaultprob) <- c(length(defaultprob),1)

-    }

-    stopifnot(length(Z)==length(w),

-              nrow(defaultprob)==length(issuerweights),

-              nrow(defaultprob)==length(recov))

-    L <- matrix(0, N, ncol(defaultprob))

-    R <- matrix(0, N, ncol(defaultprob))

-    rho <- rep(rho, length(issuerweights))

-    r <- .C("BCloss_recov_dist", defaultprob, as.integer(nrow(defaultprob)),

-            as.integer(ncol(defaultprob)), 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))

-}

-

-BCER <- function(defaultprob, issuerweights, recov, K, rho, Z, w,

-                 N=length(issuerweights)+1, defaultflag=FALSE){

-    stopifnot(length(Z)==length(w),

-              nrow(defaultprob)==length(issuerweights))

-

-    rho <- rep(rho, length(issuerweights))

-    ELt <- numeric(ncol(defaultprob))

-    ERt <- numeric(ncol(defaultprob))

-    r <- .C("BCloss_recov_trunc", defaultprob, as.integer(nrow(defaultprob)),

-            as.integer(ncol(defaultprob)),

-            as.double(issuerweights), as.double(recov), as.double(Z), as.double(w),

-            as.integer(length(Z)), as.double(rho), as.integer(N), as.double(K),

-            as.logical(defaultflag), ELt=ELt, ERt=ERt)

-    return(list(ELt=r$ELt, ERt=r$ERt))

-}

+## 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
+
+#' Gauss-Hermite quadrature weights
+#'
+#' \code{GHquad} computes the quadrature weights for integrating against a
+#' Gaussian distribution.
+#'
+#' if f is a function, then \eqn{\sum_{i=1}^n f(Z_i)w_i \approx
+#' \frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty f(x)e^{-\frac{x^2}{2}}\,dx}.
+#' @param n Integer, the number of nodes
+#' @return A list with two components:
+#'  \item{Z}{the list of nodes}
+#'  \item{w}{the corresponding weights}
+#'
+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)
+}
+
+#' Loss distribution of a portfolio
+#'
+#' \code{lossdistrib} computes the probability distribution of a sum
+#' of independent Bernouilli variables with unequal probabilities
+#' (also called the Poisson-Binomial distribution).
+#'
+#' We compute the probability distribution of \eqn{S = \sum_{i=1}^n X_i}
+#' where \eqn{X_i} is Bernouilli(\eqn{p_i}). This uses the simple recursive
+#' algorithm of Andersen, Sidenius and Basu
+#' @param p Numeric vector, the vector of success probabilities
+#' @return A vector q such that \eqn{q_k=\Pr(S=k)}
+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)
+}
+
+#'Convolution of two discrete probability distributions
+#'
+#' This is the probability of the sum of the two random variables
+#' assuming they are independent.
+#' @param dist1 Numeric vector of probabilities
+#' @param dist2 Numeric vector of probabilities
+#' @return Convolution of the two vectors
+convolve <- function(dist1, dist2) {
+    n1 <- length(dist1)
+    n2 <- length(dist2)
+    dist1 <- c(dist1, rep(0, n2-1))
+    dist2 <- c(dist2, rep(0, n1-1))
+    Re(fft(fft(dist1) * fft(dist2), inverse=T))/length(dist1)
+}
+
+#' Loss distribution of a portfolio
+#'
+#' \code{lossdistrib.fft} computes the probability distribution of a sum
+#' of independent Bernouilli variables with unequal probabilities
+#' (also called the Poisson-Binomial distribution).
+#'
+#' We compute the probability distribution of \eqn{S = \sum_{i=1}^n X_i}
+#' where \eqn{X_i} is Bernouilli(\eqn{p_i}).
+#' This uses the FFT, thus complexity is of order \eqn{O(n \log(n))},
+#' compared to \eqn{O(n^2)} for the recurvise algotithm.
+#' @param p Numeric vector, the vector of success probabilities
+#' @return A vector such that \eqn{q_k=\Pr(S=k)}
+lossdistrib.fft <- function(p) {
+    ## haven't tested when p is not a poiwer of 2.
+    if(length(p) == 1){
+        c(1-p, p)
+    }else {
+        convolve(lossdistrib.fft(p[1:(length(p)/2)]),
+                 lossdistrib.fft(p[-(1:(length(p)/2))]))
+    }
+}
+
+#' Loss distribution of a portfolio
+#'
+#' \code{lossdistrib2} computes the probability distribution of a sum
+#' of independent Bernouilli variables with unequal probabilities.
+#'
+#' We compute the probability distribution of \eqn{L = \sum_{i=1}^n w_i S_i X_i}
+#' where \eqn{X_i} is Bernouilli(\eqn{p_i}). If \code{defaultflag} is TRUE, we
+#' compute the distribution of \eqn{D = \sum_{i=1}^n w_i X_i} instead.
+#' This a recursive algorithm with first order correction for discretization.
+#' Complexity is of ordier \eqn{O(nN)}, so linear for a given grid size.
+#' @param p Numeric, vector of default probabilities
+#' @param w Numeric, vector of weights
+#' @param S Numeric, vector of severities
+#' @param N Integer, number of ticks in the grid
+#' @param defaultflag Boolean, if TRUE, we compute the default distribution
+#' (instead of the loss distribution).
+#' @return a Numeric vector of size \code{N} computing the loss (resp.
+#' default) distribution if \code{defaultflag} is FALSE (resp. TRUE).
+lossdistrib2 <- function(p, w, S, N, defaultflag=FALSE){
+    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)
+}
+
+#' Loss distribution truncated version
+#'
+#' \code{lossdistrib2.truncated} computes the probability distribution of a sum
+#' of independent Bernouilli variables with unequal probabilities up
+#' to a cutoff N.
+#'
+#' This is actually slower than \code{lossdistrib2}, but in C this is
+#' twice as fast. For high severities, M can become bigger than the cutoff, and
+#' there is some probability mass escaping.
+#' @param p Numeric, vector of default probabilities
+#' @param w Numeric, vector of weights
+#' @param S Numeric, vector of severities
+#' @param N Integer, number of ticks in the grid
+#' @param cutoff Integer, where to stop computing the exact probabilities
+#' @return a Numeric vector of size \code{N} computing the loss distribution
+lossdistrib2.truncated <- function(p, w, S, N, cutoff=N){
+    n <- length(p)
+    lu <- 1/(N-1)
+    q <- rep(0, cutoff)
+    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)
+}
+
+#' Recovery distribution of a portfolio
+#'
+#' \code{recovdist} computes the recovery distribution of portfolio
+#' described by a vector of default probabilities, and prepay probabilities.
+#' \eqn{R=\sum_{i=1}^n w_i X_i} where \eqn{X_i=0} w.p. \eqn{1-dp_i-pp_i},
+#' \eqn{X_i=1-S_i} with probability \eqn{dp_i}, and \eqn{X_i=1} w.p. \eqn{pp_i}
+#'
+#' It is a recursive algorithm with first-order correction. For a unit of loss
+#' \eqn{lu}, each non-zero value \eqn{v} is interpolated on the grid
+#' as the pair of values
+#' \eqn{\left\lfloor\frac{v}{lu}\right\rfloor}  and
+#' \eqn{\left\lceil\frac{v}{lu}\right\rceil} so that \eqn{X_i} has
+#' four non zero values.
+#' @param dp Numeric, vector of default probabilities
+#' @param pp Numeric, vector of prepay probabilities
+#' @param w Numeric, vector of weights
+#' @param S Numeric, vector of severities
+#' @param N Integer, number of ticks in the grid
+#' @return a Numeric vector of size \code{N} computing the recovery distribution
+recovdist <- function(dp, pp, w, S, N){
+    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)
+}
+
+#' Joint loss-recovery distributionrecursive algorithm with first order correction to compute the joint
+#' probability distribution of the loss and recovery.
+#'
+#' For high severities, M can become bigger than N, and there is
+#' some probability mass escaping.
+#' @param p Numeric, vector of default probabilities
+#' @param w Numeric, vector of weights
+#' @param S Numeric, vector of severities
+#' @param N Integer, number of ticks in the grid
+#' @param defaultflab Logical, whether to return the loss or default distribution
+#' @return q Matrix of joint loss, recovery probability distribution
+#' colSums(q) is the recovery distribution marginal
+#' rowSums(q) is the loss distribution marginal
+lossdist.joint <- function(p, w, S, N, defaultflag=FALSE){
+    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.integer(N), as.logical(defaultflag), q = double(N))$q
+}
+
+lossdistCZ <- function(p, w, S, N, defaultflag=FALSE, rho, Z){
+    ##S is of size (length(p), length(Z))
+    stopifnot(length(rho)==length(p),
+              length(rho)==length(w),
+              nrow(S)==length(p),
+              ncol(S)==length(Z))
+    .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, defaultflag=FALSE){
+    ## truncated version of lossdistrib
+    ## q[i] is 0 for i>=T
+    .C("lossdistrib_truncated", as.double(p), as.integer(length(p)),
+       as.double(w), as.double(S), as.integer(N), as.integer(T), as.logical(defaultflag),
+       q = double(N))$q
+}
+
+exp.trunc <- function(p, w, S, N, K){
+    ## computes E[(K-L)^+]
+    r <- 0
+    .C("exp_trunc", as.double(p), as.integer(length(p)),
+       as.double(w), as.double(S), as.integer(N), as.double(K), res = r)$res
+}
+
+rec.trunc <- function(p, w, S, N, K){
+    ## computes E[(K-(1-R))^+] = E[(\tilde K- \bar R)]
+    ## where \tilde K = K-sum_i w_i S_i and \bar R=\sum_i w_i R_i (1-X_i)
+    Ktilde <- K-crossprod(w, S)
+    if(Ktilde < 0){
+        return( 0 )
+    }else{
+        return( exp.trunc(1-p, w, 1-S, N, Ktilde) )
+    }
+}
+
+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.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 vector of doubles of size n.credits
+    ## Z is a vector of length nZ
+    ## w is a vector if length wZ
+
+    r <- .C("lossdistrib_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)), output = matrix(0,N,N))
+    return(r$output)
+}
+
+lossrecovdist <- function(defaultprob, prepayprob, w, S, N, defaultflag=FALSE, useC=TRUE){
+    lossdistrib2 <- if(useC) lossdistC
+    recovdist <- if(useC) recovdistC
+    if(missing(prepayprob)){
+        L <- lossdistrib2(defaultprob, w, S, N, defaultflag)
+        R <- lossdistrib2(defaultprob, w, 1-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))
+    lossdist.joint <- if(useC) lossdistC.jointblas
+    lossdist.prepay.joint <- if(useC) lossdistC.prepay.jointblas
+    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.jointblas(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag)
+        }
+    }
+    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
+    if(p0==0){
+        return(0)
+    }
+    if(rho == 1){
+        distw <- distr::DiscreteDistribution(Z, w)
+        return(distr::pnorm(distr::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){
+    stopifnot(length(Z)==length(w))
+    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)
+}
+
+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
+    }
+    stopifnot(length(Z)==length(w),
+              nrow(defaultprob) == length(issuerweights))
+
+    ## 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){
+    if(is.null(dim(defaultprob))){
+        dim(defaultprob) <- c(length(defaultprob),1)
+    }
+    stopifnot(length(Z)==length(w),
+              nrow(defaultprob)==length(issuerweights),
+              nrow(defaultprob)==length(recov))
+    L <- matrix(0, N, ncol(defaultprob))
+    R <- matrix(0, N, ncol(defaultprob))
+    rho <- rep(rho, length(issuerweights))
+    r <- .C("BCloss_recov_dist", defaultprob, as.integer(nrow(defaultprob)),
+            as.integer(ncol(defaultprob)), 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))
+}
+
+BCER <- function(defaultprob, issuerweights, recov, K, rho, Z, w,
+                 N=length(issuerweights)+1, defaultflag=FALSE){
+    stopifnot(length(Z)==length(w),
+              nrow(defaultprob)==length(issuerweights))
+
+    rho <- rep(rho, length(issuerweights))
+    ELt <- numeric(ncol(defaultprob))
+    ERt <- numeric(ncol(defaultprob))
+    r <- .C("BCloss_recov_trunc", defaultprob, as.integer(nrow(defaultprob)),
+            as.integer(ncol(defaultprob)),
+            as.double(issuerweights), as.double(recov), as.double(Z), as.double(w),
+            as.integer(length(Z)), as.double(rho), as.integer(N), as.double(K),
+            as.logical(defaultflag), ELt=ELt, ERt=ERt)
+    return(list(ELt=r$ELt, ERt=r$ERt))
+}