new function using truncated dists

1 files changed, 52 insertions, 28 deletions

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

+}