1 files changed, 47 insertions, 34 deletions
diff --git a/tranche_functions.R b/tranche_functions.R
index 7329c971..479749a8 100644
--- a/tranche_functions.R
+++ b/tranche_functions.R
@@ -35,18 +35,23 @@ lossdistrib.fft <- function(p){
     return(1/(n+1)*Re(fft(v)))
 }
 
-lossdistrib2 <- function(p, w, S, 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
+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){
-        d <- S[i] * w[i] / lu
+        if(defaultflag){
+            d <- w[i] /lu
+        }else{
+            d <- S[i] * w[i] / lu
+        }
         d1 <- floor(d)
         d2 <- ceiling(d)
         p1 <- p[i]*(d2-d)
@@ -127,7 +132,7 @@ recovdist <- function(dp, pp, w, S, N){
     return(q)
 }
 
-lossdistrib.joint <- function(p, w, S, N){
+lossdistrib.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:
@@ -135,6 +140,7 @@ lossdistrib.joint <- function(p, w, S, N){
     ##   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
@@ -144,7 +150,11 @@ lossdistrib.joint <- function(p, w, S, N){
     q <- matrix(0, N, N)
     q[1,1] <- 1
     for(k in 1:n){
-        x <- S[k] * w[k]/lu
+        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)
@@ -161,7 +171,7 @@ lossdistrib.joint <- function(p, w, S, N){
     return(q)
 }
 
-lossdistribprepay.joint <- function(dp, pp, w, S, N){
+lossdistribprepay.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:
@@ -170,6 +180,7 @@ lossdistribprepay.joint <- function(dp, pp, w, S, N){
     ##   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
@@ -201,13 +212,13 @@ lossdistribprepay.joint <- function(dp, pp, w, S, N){
     return(q)
 }
 
-lossdistribC <- function(p, w, S, N){
+lossdistribC <- function(p, w, S, N, defaultflag){
     ## C version of lossdistrib2, roughly 50 times faster
     if(!is.loaded("lossdistrib")){
         dyn.load(paste0("lossdistrib", .Platform$dynlib.ext))
     }
     .C("lossdistrib", as.double(p), as.integer(length(p)),
-       as.double(w), as.double(S), as.integer(N), q = double(N))$q
+       as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), q = double(N))$q
 }
 
 lossdistribC.truncated <- function(p, w, S, N, T=N){
@@ -229,78 +240,78 @@ recovdistC <- function(dp, pp, w, S, N){
        as.double(w), as.double(S), as.integer(N), q = double(N))$q
 }
 
-lossdistribC.joint <- function(p, w, S, N){
+lossdistribC.joint <- function(p, w, S, N, defaultflag=FALSE){
     ## C version of lossdistrib.joint, roughly 20 times faster
     if(!is.loaded("lossdistrib_joint")){
         dyn.load(paste0("lossdistrib", .Platform$dynlib.ext))
     }
     .C("lossdistrib_joint", as.double(p), as.integer(length(p)), as.double(w),
-       as.double(S), as.integer(N), q = matrix(0, N, N))$q
+       as.double(S), as.integer(N), as.logical(defaultflag), q = matrix(0, N, N))$q
 }
 
-lossdistribprepayC.joint <- function(dp, pp, w, S, N){
+lossdistribprepayC.joint <- function(dp, pp, w, S, N, defaultflag=FALSE){
     ## C version of lossdistribprepay.joint
     if(!is.loaded("lossdistrib_prepay_joint")){
         dyn.load(paste0("lossdistrib", .Platform$dynlib.ext))
     }
     .C("lossdistrib_prepay_joint", as.double(dp), as.double(pp), as.integer(length(dp)),
-       as.double(w), as.double(S), as.integer(N), q=matrix(0, N, N))$q
+       as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), q=matrix(0, N, N))$q
 }
 
-lossrecovdist <- function(defaultprob, prepayprob, w, S, N, useC=TRUE){
+lossrecovdist <- function(defaultprob, prepayprob, w, S, N, defaultflag, useC=TRUE){
     if(all(!prepayprob)){
         if(useC){
-            L <- lossdistribC(defaultprob, w, S, N)
+            L <- lossdistribC(defaultprob, w, S, N, defaultflag)
             R <- lossdistribC(defaultprob, w, 1-S, N)
         }else{
-            L <- lossdistrib2(defaultprob, w, S, N)
+            L <- lossdistrib2(defaultprob, w, S, N, defaultflag)
             R <- lossdistrib2(defaultprob, w, 1-S, N)
         }
     }else{
         if(useC){
-            L <- lossdistribC(defaultprob, w, S, N)
+            L <- lossdistribC(defaultprob, w, S, N, defaultflag)
             R <- recovdistC(defaultprob, prepayprob, w, S, N)
         }else{
-            L <- lossdistrib2(defaultprob, w, S, N)
+            L <- lossdistrib2(defaultprob, 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, useC=TRUE){
+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)))
     for(t in 1:ncol(defaultprob)){
-        temp <- lossrecovdist(defaultprob[,t], prepayprob[,t], w, S[,t], N, useC)
+        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, useC=TRUE){
+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(all(!prepayprob)){
             for(t in 1:ncol(defaultprob)){
-                Q[t,,] <- lossdistribC.joint(defaultprob[,t], w, S[,t], N)
+                Q[t,,] <- lossdistribC.joint(defaultprob[,t], w, S[,t], N, defaultflag)
             }
         }else{
             for(t in 1:ncol(defaultprob)){
-                Q[t,,] <- lossdistribprepayC.joint(defaultprob[,t], prepayprob[,t], w, S[,t], N)
+                Q[t,,] <- lossdistribprepayC.joint(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag)
             }
         }
     }else{
         if(all(!prepayprob)){
             for(t in 1:ncol(defaultprob)){
-                Q[t,,] <- lossdistrib.joint(defaultprob[,t], w, S[,t], N)
+                Q[t,,] <- lossdistrib.joint(defaultprob[,t], w, S[,t], N, defaultflag)
             }
         }else{
             for(t in 1:ncol(defaultprob)){
-                Q[t,,] <- lossdistribprepay.joint(defaultprob[,t], prepayprob[,t], w, S[,t], N)
+                Q[t,,] <- lossdistribprepay.joint(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag)
             }
         }
     }
@@ -480,7 +491,7 @@ BClossdist <- function(SurvProb, issuerweights, recov, rho, N=length(recov)+1,
 }
 
 BClossdistC <- function(SurvProb, issuerweights, recov, rho,
-                        N=length(issuerweights)+1, n.int=100){
+                        N=length(issuerweights)+1, n.int=100, defaultflag){
     if(!is.loaded("BClossdist")){
         dyn.load(paste0("lossdistrib", .Platform$dynlib.ext))
     }
@@ -491,7 +502,7 @@ BClossdistC <- function(SurvProb, issuerweights, recov, rho,
     R <- matrix(0, N, dim(SurvProb)[2])
     r <- .C("BClossdist", SurvProb, as.integer(dim(SurvProb)[1]), as.integer(dim(SurvProb)[2]),
             as.double(issuerweights), as.double(recov), as.double(Z), as.double(w),
-            as.integer(n.int), as.double(rho), as.integer(N), L=L, R=R)
+            as.integer(n.int), as.double(rho), as.integer(N), as.logical(defaultflag), L=L, R=R)
     return(list(L=r$L,R=r$R))
 }
 
@@ -680,12 +691,14 @@ MFlossrecovdist2 <- function(cl, w, Z, rho, defaultprob, defaultprobmod,
         ppshocked <- apply(prepayprobmod, 2, shockprob, rho=rho, Z=-Z[i])
         S <- 1 - Rstoch[i,,]
         dist <- lossrecovdist.joint.term(dpshocked, ppshocked, issuerweights, S, Ngrid)
+        gc()
         return(dist)
     }
-    ## temp <- sapply(1:length(w), parf)
-    clusterExport(cl, list("Rstoch", "defaultprobmod", "prepayprobmod"), envir=environment())
-    temp <- parSapply(cl, 1:length(w), parf)
-    clusterCall(cl, gc)
+    temp <- sapply(1:length(w), parf)
+    ## clusterExport(cl, list("Rstoch", "defaultprobmod", "prepayprobmod"), envir=environment())
+    ## temp <- parSapply(cl, 1:length(w), parf)
+    ##gc()
+    ## 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))