big cleanup

1 files changed, 13 insertions, 575 deletions

diff --git a/R/tranche_functions.R b/R/tranche_functions.R
index 97fbc83..8b10394 100644
--- a/R/tranche_functions.R
+++ b/R/tranche_functions.R
@@ -355,22 +355,14 @@ lossdistC.prepay.jointZ <- function(dp, pp, w, S, N, defaultflag = FALSE, rho, Z
 }

 

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

+    lossdistrib2 <- if(useC) lossdistC

+    recovdist <- if(useC) recovdistC

     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)

-        }

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

-        }

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

+        R <- recovdist(defaultprob, prepayprob, w, S, N)

     }

     return(list(L=L, R=R))

 }

@@ -398,28 +390,17 @@ lossrecovdist.term <- function(defaultprob, prepayprob, w, S, N, defaultflag=FAL
 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.jointblas(defaultprob[,t], w, S[,t], N, defaultflag)

-            }

-        }else{

-            for(t in 1:ncol(defaultprob)){

-                Q[t,,] <- lossdistC.prepay.jointblas(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag)

-            }

+    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{

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

-            }

+        for(t in 1:ncol(defaultprob)){

+            Q[t,,] <- lossdist.prepay.jointblas(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag)

         }

     }

-    gc()

     return(Q)

 }

 

@@ -530,123 +511,6 @@ stochasticrecovC <- function(R, Rtilde, Z, w, rho, porig, pmod){
     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) )

-}

-

 BClossdist <- function(defaultprob, issuerweights, recov, rho, Z, w,

                        N=length(recov)+1, defaultflag=FALSE, n.int=500){

   if(missing(Z)){

@@ -684,429 +548,3 @@ BClossdistC <- function(defaultprob, issuerweights, recov, rho, Z, 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, complement=FALSE){

-    ## computes the protection leg and couponleg of a 0-K tranche

-    ## if complement==TRUE, computes the protection leg and coupon leg of a K-1 tranche

-    if((K==0 && !complement) || (K==1 && complement)){

-        return(list(cl=0, pl=0))

-    }else if((K==1 && !complement) || (K==0 && complement)){

-        return(BCindex.pv(index))

-    }else{

-        dist <- BClossdistC(index$defaultprob, index$issuerweights, index$recov, rho, index$Z, index$w, index$N)

-        if(complement){

-            return(list(cl=tranche.cl(dist$L, dist$R, index$cs, K, 1),

-                        pl=tranche.pl(dist$L, index$cs, K, 1)))

-        }else{

-            return(list(cl=tranche.cl(dist$L, dist$R, index$cs, 0, K),

-                        pl=tranche.pl(dist$L, index$cs, 0, K)))

-        }

-    }

-}

-

-BCtranche.pv <- function(index, protection=FALSE, complement=FALSE){

-    ## computes the protection leg, couponleg, and bond price of a tranche

-    ## in the base correlation setting

-    ## if complement=FALSE compute the pvs starting from 0 (bottom-up skew)

-    ## if complement=TRUE compute the pvs starting from 1 (top-down skew)

-    pl <- rep(0, length(index$rho))

-    cl <- rep(0, length(index$rho))

-    for(i in seq_along(index$rho)){

-        temp <- BCtranche.legs(index, index$K[i], index$rho[i], complement)

-        pl[i] <- temp$pl

-        cl[i] <- temp$cl

-    }

-    dK <- diff(index$K)

-    plvec <- diff(pl)/dK

-    clvec <- diff(cl)/dK*index$tranches$running

-    if(complement){

-        plvec <- -plvec

-        clvec <- -clvec

-    }

-    if(protection){

-        bp <- -plvec-clvec

-    }else{

-        bp <- 1+plvec+clvec

-    }

-    return(list(pl=plvec, cl=clvec, bp=bp))

-}

-

-MFtranche.pv <- function(index, dist, protection=FALSE){

-    n.tranches <- length(index$K)-1

-    pl <- rep(0, n.tranches)

-    cl <- rep(0, n.tranches)

-    for(i in seq.int(n.tranches)){

-        pl[i] <- tranche.pl(dist$L, index$cs, index$K[i], index$K[i+1])

-        cl[i] <- tranche.cl(dist$L, dist$R, index$cs, index$K[i], index$K[i+1])

-    }

-    dK <- diff(index$K)

-    plvec <- pl/dK

-    clvec <- cl/dK*index$tranches$running

-    if(protection){

-        bp <- -plvec-clvec

-    }else{

-        bp <- 1+plvec+clvec

-    }

-    return(list(pl=plvec, cl=clvec, bp=bp))

-}

-

-adjust.skew <- function(index1, index2, method="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 <- index1$K[-c(1,length(index1$K))]

-    K2 <- index2$K[-c(1,length(index2$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 and 0.01

-        f <- splinefun(K1, index1$rho[-c(1, length(index1$rho))], "natural")

-        return(pmax(pmin(f(x), 0.99), 0.01))

-    }

-    if(method=="ATM"){

-        K1eq <- el1/el2 * K2

-    }else if(method == "TLP"){

-        K1eq <- c()

-        m <- max(K2) + 0.3

-        for(K2val in K2){

-            prog <- optimize(aux, interval=c(0,m),

-                             index1=index1, el1=el1, skew=skew,

-                             index2=index2, el2=el2, K2=K2val)

-            K1eq <- c(K1eq, prog$minimum)

-        }

-    }else if (method=="PM"){

-        K1eq <- c()

-        m <- max(K2) + 0.25

-        for(K2val in K2){

-            prog <- optimize(aux2, interval=c(0, m),

-                             index1=index1, skew=skew, index2=index2, K2=K2val)

-            K1eq <- c(K1eq, prog$minimum)

-        }

-     }

-    return(c(NA, 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, shortened=4, complement=FALSE, method="ATM"){

-    temp <- BCtranche.pv(index, complement=complement)

-    rho.adj <- theta.adjust.skew(index, shortened, method)

-    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, complement=complement)

-    temp3 <- BCtranche.delta(index, complement=complement)

-    return(data.frame(theta=temp2$bp-temp$bp+index$tranches$running,

-                      fw.delta=temp3$delta))

-}

-

-BCtranche.delta <- function(index, complement=FALSE){

-    ## computes the tranche delta (on current notional) by doing a proportional

-    ## blip of all the curves

-    ## if complement is False, then computes deltas bottom-up

-    ## if complement is True, then computes deltas top-down

-    eps <- 1e-4

-    index$N <- 301 ## for gamma computations we need all the precision we can get

-    ## we build a lit of 4 indices with various shocks

-    index.list <- lapply(c(0, eps, -eps, 2*eps), function(x){

-        if(x==0){

-            return(index)

-        }else{

-            newindex <- index

-            newindex$portfolio <- tweakportfolio(newindex$portfolio, x)

-            newindex$defaultprob <- 1 - SPmatrix(newindex$portfolio, length(index$cs$dates))

-            return(newindex)

-        }

-    })

-    bp <- matrix(0, length(index$K)-1, length(index.list))

-    indexbp <- rep(0, length(index.list))

-    for(j in seq_along(index.list)){

-        indexbp[j] <- BCindex.pv(index.list[[j]])$bp

-        bp[,j] <- BCtranche.pv(index.list[[j]], complement=complement)$bp

-    }

-

-    deltas <- (bp[,2]-bp[,3])/(indexbp[2]-indexbp[3])*tranche.factor(index)/index$factor

-    deltasplus <- (bp[,4]-bp[,1])/(indexbp[4]-indexbp[1])*tranche.factor(index)/index$factor

-    gammas <- (deltasplus-deltas)/(indexbp[2]-indexbp[1])/100

-

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

-}

-

-MFtranche.delta <- function(index){

-    ## computes the tranche delta (on current notional) by doing a proportional

-    ## blip of all the curves

-    ## if complement is False, then computes deltas bottom-up

-    ## if complement is True, then computes deltas top-down

-    eps <- 1e-4

-    index$Ngrid <- 301 ## for gamma computations we need all the precision we can get

-    ## we build a lit of 4 indices with various shocks

-    index.list <- lapply(c(0, eps, -eps, 2*eps), function(x){

-        if(x==0){

-            return(index)

-        }else{

-            newindex <- index

-            newindex$portfolio <- tweakportfolio(newindex$portfolio, x)

-            newindex$defaultprob <- 1 - SPmatrix(newindex$portfolio, length(index$cs$dates))

-            return(newindex)

-        }

-    })

-    bp <- matrix(0, length(index$K)-1, length(index.list))

-    indexbp <- rep(0, length(index.list))

-    for(j in seq_along(index.list)){

-        indexbp[j] <- BCindex.pv(index.list[[j]])$bp

-        dist <- MFdist(index.list[[j]])

-        bp[,j] <- BCtranche.pv(index.list[[j]], dist)

-    }

-

-    deltas <- (bp[,2]-bp[,3])/(indexbp[2]-indexbp[3])*tranche.factor(index)/index$factor

-    deltasplus <- (bp[,4]-bp[,1])/(indexbp[4]-indexbp[1])*tranche.factor(index)/index$factor

-    gammas <- (deltasplus-deltas)/(indexbp[2]-indexbp[1])/100

-

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

-}

-

-BCtranche.corr01 <- function(index, eps=0.01, complement=FALSE){

-    ##does a parallel shift of the skew and computes the change in pv

-    before <- BCtranche.pv(index, complement=complement)

-    index$rho[-1] <- index$rho[-1]+eps

-    after <- BCtranche.pv(index, complement=complement)

-    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( ELvec[length(ELvec)] )

-   }else{

-       return( sum(index$cs$df[1:Ncol]*diff(c(0, ELvec))) )

-   }

-}

-

-BCindex.pv <- function(index, discounted=TRUE, shortened=0){

-    Ncol <- ncol(index$defaultprob)-shortened

-    ELvec <- as.numeric(crossprod(index$issuerweights * (1-index$recov), index$defaultprob[,1:Ncol]))

-    size <- 1-as.numeric(crossprod(index$issuerweights, index$defaultprob[,1:Ncol]))

-    sizeadj <- 0.5*(c(1, size[-length(size)])+size)

-    if(!discounted){

-        pl <- -ELvec[length(ELvec)]

-        cl <- as.numeric(crossprod(index$cs$coupons[1:Ncol], sizeadj))

-        bp <- 1+cl+pl

-    }else{

-        pl <-  -sum(index$cs$df[1:Ncol]* diff(c(0, ELvec)))

-        cl <- as.numeric(crossprod(index$cs$coupons[1:Ncol], sizeadj * index$cs$df[1:Ncol] ))

-        bp <- 1+cl+pl

-    }

-    bp <- 1+cl*index$quotes$spread+pl

-    return(list(pl=pl, cl=cl, bp=bp))

-}

-

-ELtrunc <- function(index, K, rho){

-    ## computes the expected loss of a portfolio below strike K

-    ## could be written faster by using a truncated version of lossdist

-    ## index should be a list with issuerweights, recov, defaultprob and cs parameters

-    dist <- BClossdistC(index$defaultprob, index$issuerweights, index$recov, rho, index$Z, index$w, index$N)

-    return( -tranche.pl(dist$L, index$cs, 0, K))

-}

-

-Probtrunc <- function(index, K, rho){

-    dist <- BClossdistC(index$defaultprob, index$issuerweights, index$recov, rho, index$Z, index$w, index$N)

-    p <- cumsum(dist$L[,ncol(dist$L)])

-    support <- seq(0, 1, length=index$N)

-    probfun <- splinefun(support, p, method="hyman")

-    return(probfun(K))

-}

-

-

-BCstrikes <- function(index, K, rho) {

-    ## computes the strikes as a percentage of expected loss

-    ## Kmodified is the current attachment points (adjusted for losses)

-    el <- EL(index)

-    ELvec <- c()

-    for(i in 2:length(K)){

-        ELvec <- c(ELvec, ELtrunc(index, K[i], rho[i]))

-    }

-    return( ELvec/el )

-}

-

-tranche.factor <- function(index){

-    ## compute the factor to convert from delta on current notional to delta on original notional

-    ## K1 and K2 original strikes

-    return( diff(index$K)/diff(index$K.orig)*index$factor )

-}

-

-MFupdate.prob <- function(Z, w, rho, defaultprob){

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

-}

-

-MFrecovery <- function(index, defaultprobmod){

-    n.credit <- length(index$issuerweights)

-    n.int <- length(index$Z)

-    Rstoch <- array(0, dim=c(n.credit, n.int, ncol(index$defaultprob)))

-    rho <- rep(0.45, n.credit)

-    for(t in 1:ncol(index$defaultprob)){

-        for(i in 1:n.credit){

-            Rstoch[i,,t] <- stochasticrecovC(index$recov[i], 0, index$Z, index$w.mod,

-                                             rho[i], index$defaultprob[i,t], defaultprobmod[i,t])

-        }

-    }

-    return( Rstoch )

-}

-

-MFlossdist <- function(index){

-    n.credit <- length(index$issuerweights)

-    rho <- rep(0.45, n.credit)

-    defaultprobmod <- MFupdate.probC(index$Z, index$w.mod, rho, index$defaultprob)

-    n.credit <- length(index$issuerweights)

-    n.int <- length(index$Z)

-    Rstoch <- MFrecovery(index, defaultprobmod)

-    Lw <- matrix(0, index$N, n.int)

-    Rw <- matrix(0, index$N, n.int)

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

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

-    for(t in 1:ncol(index$defaultprob)){

-        S <- 1 - Rstoch[,,t]

-        Lw <- lossdistCZ(defaultprobmod[,t], index$issuerweights, S, index$N, 0, rho, index$Z)

-        Rw <- lossdistCZ(defaultprobmod[,t], index$issuerweights, 1-S, Ngrid, 0, rho, index$Z)

-        L[,t] <- Lw%*%index$w.mod

-        R[,t] <- Rw%*%index$w.mod

-    }

-    return(list(L=L, R=R))

-}