1 files changed, 90 insertions, 37 deletions
diff --git a/tranche_functions.R b/tranche_functions.R
index 9fdd3dcb..1b20fdc7 100644
--- a/tranche_functions.R
+++ b/tranche_functions.R
@@ -1,7 +1,7 @@
 library(statmod)
 
 lossdistrib <- function(p){
-    #basic recursive algorithm of Andersen, Sidenius and Basu
+    ## basic recursive algorithm of Andersen, Sidenius and Basu
     n <- length(p)
     q <- rep(0, (n+1))
     q[1] <- 1
@@ -13,19 +13,20 @@ lossdistrib <- function(p){
 }
 
 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)))
+    ## computes loss distribution using the fft
+    ## complexity is of order O(n*m)+O(m*log(m))
+    ## where m is the size of the grid and n the number of probabilities.
+    ## this is slower than the recursive algorithm
+    theta <- 2*pi*1i*(0:n)/(n+1)
+    Phi <- 1 - p + p%o%exp(theta)
+    v <- apply(Phi, 2, prod)
+    return(1/(n+1)*Re(fft(v)))
 }
 
-lossdistrib2 <- function(p, S, lu){
+lossdistrib2 <- function(p, w, S, lu){
     #recursive algorithm with first order correction
     #p vector of default probabilities
+    #w vector of weigths
     #S vector of severities
     #lu loss unit
     n <- length(p)
@@ -34,7 +35,7 @@ lossdistrib2 <- function(p, S, lu){
     q <- rep(0, N+eps)
     q[1] <- 1
     for(i in 1:n){
-        d <- S[i]/(n*lu)
+        d <- S[i] * w[i] / lu
         d1 <- floor(d)
         d2 <- ceiling(d)
         p1 <- p[i]*(d2-d)
@@ -47,7 +48,7 @@ lossdistrib2 <- function(p, S, lu){
     return(q)
 }
 
-lossdistrib3 <- function(dp, pp, w, S, lu){
+recovdist <- function(dp, pp, w, S, lu){
     ## 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]
@@ -64,13 +65,14 @@ lossdistrib3 <- function(dp, pp, w, S, lu){
         d1 <- w[i]*(1-S[i])/lu
         d1l <- floor(d1)
         d1u <- ceiling(d1)
-        d2 <- w[i]/lu
+        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
+        cat(dp1, dp2, pp1, pp2, "\n")
         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)])
@@ -80,19 +82,20 @@ lossdistrib3 <- function(dp, pp, w, S, lu){
     return(q)
 }
 
-lossdistrib.joint <- function(p, S, lu){
-    #recursive algorithm with first order correction
-    #to compute the joint probability distribition of the loss and recovery
-    #p vector of default probabilities
-    #S vector of severities
-    #lu loss unit
+lossdistrib.joint <- function(p, w, S, lu){
+    ## recursive algorithm with first order correction
+    ## to compute the joint probability distribution of the loss and recovery
+    ## p vector of default probabilities
+    ## w vector of issuer weights
+    ## S vector of severities
+    ## lu loss unit
     n <- length(p)
     N <- ceiling(1/lu + 1)
     q <- matrix(0, N, N)
     q[1,1] <- 1
     for(k in 1:n){
-        x <- S[k]/(n*lu)
-        y <- (1-S[k])/(n*lu)
+        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))
@@ -108,34 +111,84 @@ lossdistrib.joint <- function(p, S, lu){
     return(q)
 }
 
-lossdistribC <- function(p, S, lu){
-    #C version of lossdistrib2, roughly 50 times faster
+lossdistribprepay.joint <- function(dp, pp, w, S, lu){
+    ## recursive algorithm with first order correction
+    ## to compute the joint probability distribition of the loss and recovery
+    ## dp vector of default probabilities
+    ## pp vector of prepay probabilities
+    ## w vector of issuer weights
+    ## S vector of severities
+    ## lu loss unit
+    n <- length(dp)
+    N <- ceiling(1/lu + 1)
+    q <- matrix(0, N, N)
+    q[1,1] <- 1
+    for(k in 1:n){
+        x <- S[k] * w[k]/lu
+        y1 <- (1-S[k]) * w[k]/lu
+        y2 <- w[k]/lu
+        i <- floor(x)
+        j1 <- floor(y1)
+        j2 <- floor(y2)
+        weights <- c((i+1-x)*(j1+1-y1), (i+1-x)*(y1-j1), (x-i)*(y1-j1), (j1+1-y1)*(x-i))
+        dpsplit <- dp[k] * weights
+        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)]
+        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)
+}
+
+
+lossdistribC <- function(p, w, S, lu){
+    ## C version of lossdistrib2, roughly 50 times faster
     dyn.load("lossdistrib.dll")
     .C("lossdistrib", as.double(p), as.integer(length(p)),
-       as.double(S), as.double(lu), q = double(1/lu+1))$q
+       as.double(w), as.double(S), as.double(lu), q = double(1/lu+1))$q
 }
 
-lossdistribC.joint <- function(p, S, lu){
-    #C version of lossdistrib.joint, roughly 20 times faster
+recovdistC <- function(dp, pp, w, S, lu){
+    ## C version of recovdist
+    dyn.load("lossdistrib2.dll")
+    .C("recovdist", as.double(dp), as.double(pp), as.integer(length(dp)),
+       as.double(w), as.double(S), as.double(lu), q = double(ceiling(1/lu+1)))$q
+}
+
+lossdistribC.joint <- function(p, w, S, lu){
+    ## C version of lossdistrib.joint, roughly 20 times faster
     N <- ceiling(1/lu+1)
-    dyn.load("lossdistrib.dll")
-    .C("lossdistrib2", as.double(p), as.integer(length(p)), as.double(S),
-       as.double(lu), q = matrix(0, N, N))$q
+    dyn.load("lossdistrib2.dll")
+    .C("lossdistrib_joint", as.double(p), as.integer(length(p)), as.double(w),
+       as.double(S), as.double(lu), q = matrix(0, N, N))$q
+}
+
+lossdistribprepayC.joint <- function(dp, pp, w, S, lu){
+    ## C version of lossdistribprepay.joint
+    dyn.load("lossdistrib2.dll")
+    N <- ceiling(1/lu+1)
+    .C("lossdistrib_prepay_joint", as.double(dp), as.double(pp), as.integer(length(dp)),
+       as.double(w), as.double(S), as.double(lu), q=matrix(0, N, N))$q
 }
 
-lossrecovdist <- function(defaultprob, prepayprob, S, lu, useC=TRUE){
+lossrecovdist <- function(defaultprob, prepayprob, w, S, lu, useC=TRUE){
     if(prepayprob==0){
         if(useC){
-            L <- lossdistribC(defaultprob, S, lu)
-            R <- lossdistribC(defaultprob, 1-S, lu)
+            L <- lossdistribC(defaultprob, w, S, lu)
+            R <- lossdistribC(defaultprob, w, 1-S, lu)
         }else{
-            L <- lossdistrib2(defaultprob, S, lu)
-            R <- lossdistrib2(defaultprob, 1-S, lu)
+            L <- lossdistrib2(defaultprob, w, S, lu)
+            R <- lossdistrib2(defaultprob, w, 1-S, lu)
         }
     }else{
-        u <- prepayprob/defaultprob
-        L <- lossdistribC(defaultprob+prepayprob, S/(1+u), lu)
-        R <- lossdistribC(defaultprob+prepayprob, (u+1-S)/(1+u), lu)
+        L <- lossdistribC(defaultprob, w, S, lu)
+        R <- recovdistC(defaultprob, prepayprob, w, S, lu)
     }
     return(list(L=L, R=R))
 }