From 6528b3cf1b976b8392e4ae9d0f1857c22276ae9a Mon Sep 17 00:00:00 2001
From: Guillaume Horel <guillaume.horel@serenitascapital.com>
Date: Wed, 11 Mar 2015 16:38:53 -0400
Subject: new function using truncated dists

---
 R/distrib.R | 80 +++++++++++++++++++++++++++++++++++++++----------------------
 1 file changed, 52 insertions(+), 28 deletions(-)

(limited to 'R')

diff --git a/R/distrib.R b/R/distrib.R
index b2c028a..536ed20 100644
--- a/R/distrib.R
+++ b/R/distrib.R
@@ -122,7 +122,7 @@ lossdistrib2.truncated <- function(p, w, S, N, cutoff=N){
     ## some probability mass escaping.
     n <- length(p)
     lu <- 1/(N-1)
-    q <- rep(0, truncated)
+    q <- rep(0, cutoff)
     q[1] <- 1
     M <- 1
     for(i in 1:n){
@@ -508,6 +508,7 @@ fit.prob <- function(Z, w, rho, p0){
 }
 
 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)
@@ -533,38 +534,61 @@ stochasticrecovC <- function(R, Rtilde, Z, w, rho, porig, pmod){
 
 BClossdist <- function(defaultprob, issuerweights, recov, rho, Z, w,
                        N=length(recov)+1, defaultflag=FALSE, n.int=500){
-  if(missing(Z)){
-    quadrature <- GHquad(n.int)
-    Z <- quadrature$Z
-    w <- quadrature$w
-  }
-  ## do not use if weights are not gaussian, results would be incorrect
-  ## since shockseverity is invalid in that case (need to use stochasticrecov)
-  LZ <- matrix(0, N, length(Z))
-  RZ <- matrix(0, N, length(Z))
-  L <- matrix(0, N, ncol(defaultprob))
-  R <- matrix(0, N, ncol(defaultprob))
-  for(t in 1:ncol(defaultprob)){
-    for(i in 1:length(Z)){
-      g.shocked <- shockprob(defaultprob[,t], rho, Z[i])
-      S.shocked <- shockseverity(1-recov, 1, Z[i], rho, defaultprob[,t])
-      temp <- lossrecovdist(g.shocked, , issuerweights, S.shocked, N)
-      LZ[,i] <- temp$L
-      RZ[,i] <- temp$R
+
+    if(missing(Z)){
+        quadrature <- GHquad(n.int)
+        Z <- quadrature$Z
+        w <- quadrature$w
     }
-    L[,t] <- LZ%*%w
-    R[,t] <- RZ%*%w
-  }
-  list(L=L, R=R)
+    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){
-    L <- matrix(0, N, dim(defaultprob)[2])
-    R <- matrix(0, N, dim(defaultprob)[2])
+    stopifnot(length(Z)==length(w),
+              nrow(defaultprob)==length(issuerweights))
+    L <- matrix(0, N, ncol(defaultprob))
+    R <- matrix(0, N, ncol(defaultprob))
     rho <- rep(rho, length(issuerweights))
-    r <- .C("BCloss_recov_dist", defaultprob, dim(defaultprob)[1], dim(defaultprob)[2],
-            as.double(issuerweights), as.double(recov), as.double(Z), as.double(w),
-            as.integer(length(Z)), as.double(rho), as.integer(N), as.logical(defaultflag), L=L, R=R)
+    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(issuwerweights)+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))
+}
-- 
cgit v1.3