2 files changed, 50 insertions, 47 deletions
diff --git a/R/calibrate_tranches_MF.R b/R/calibrate_tranches_MF.R
index 95a9082e..32396a88 100644
--- a/R/calibrate_tranches_MF.R
+++ b/R/calibrate_tranches_MF.R
@@ -30,59 +30,16 @@ if(length(args) >= 1){
 exportYC(tradedate)
 ## calibrate HY21
 ## calibrate the single names curves
-n.int <- 250
+n.int <- 100
 list2env(GHquad(n.int), envir=parent.frame())
 Ngrid <- 201
 index <- load.index("hy21", date=tradedate, "5yr", Z, w, Ngrid)
 index <- set.singlenamesdata(index, tradedate)
-index <- set.tranchedata(index, tradedate)
-
 ## load tranche data
-n.credit <- length(index$portfolio)
-index$cs$coupons <- index$cs$coupons*0.05
+index <- set.tranchedata(index, tradedate)
 ##calibrate by modifying the factor distribution
-bottomup <- 1:3
-topdown <- 2:4
-
-index$w.mod <- w
-p <- index$defaultprob
-rho <- rep(0.45, n.credit)
-result <- matrix(0, 4, n.int)
-err <- Inf
-errvec <- c()
-while(err >0.01){
-    Rstoch <- MFrecovery(index, p)
-    L <- array(0, dim=c(index$N, n.int, ncol(index$defaultprob)))
-    R <- array(0, dim=c(index$N, n.int, ncol(index$defaultprob)))
-    for(t in 1:ncol(index$defaultprob)){
-        S <- 1 - Rstoch[,,t]
-        L[,,t] <- lossdistCZ(p[,t], index$issuerweights, S, index$N, 0, rho, index$Z)
-        R[,,t] <- lossdistCZ(p[,t], index$issuerweights, 1-S, index$N, 0, rho, index$Z)
-    }
-    for(i in 1:n.int){
-        result[,i] <- tranche.pvvec(index$K, L[,i,], R[,i,], index$cs)
-    }
-    ## solve the optimization problem
-    program <- KLfit(-result[bottomup,], w, index$tranche.quotes[bottomup])
-
-    err <- 0
-    for(i in 1:n.credit){
-        for(j in 1:ncol(p)){
-            err <- err + abs(crossprod(shockprob(p[i,j], rho[i], Z), program$weight) - index$defaultprob[i,j])
-        }
-    }
-    errvec <- c(errvec, err)
-
-    ## update the new probabilities
-    p <- MFupdate.probC(Z, program$weight, rho, index$defaultprob)
-
-    errvec <- c(errvec, err)
-    index$w.mod <- program$weight
-    cat(err,"\n")
-}
-
-dist <- MFdist(index)
-
+index$w.mod <- build.MFdist(index)
+dist <- MFlossdist(index)
 write.table(data.frame(Z=Z, w=index$w.mod),
             file=file.path(root.dir, "Scenarios", "Calibration",
             paste0("calibration-", tradedate, ".csv")),
diff --git a/R/calibration.R b/R/calibration.R
index 1b802d99..b8dc9295 100644
--- a/R/calibration.R
+++ b/R/calibration.R
@@ -118,3 +118,49 @@ build.skew <- function(index, type="bottomup"){
     }
     return( rhovec )
 }
+
+build.MFdist <- function(index, type="bottomup", tol=1e-2){
+    index$w.mod <- w
+    p <- index$defaultprob
+    n.credit <- length(index$issuerweights)
+    rho <- rep(0.45, n.credit)
+    result <- matrix(0, 4, n.int)
+    n.tranches <- length(index$K)-1
+    err <- Inf
+    if(type=="bottomup"){
+        select <- 1:(n.tranches-1)
+    }else if(type=="topdown"){
+        select <- 2:n.tranches
+    }
+    while(err > tol){
+        Rstoch <- MFrecovery(index, p)
+        L <- array(0, dim=c(index$N, n.int, ncol(index$defaultprob)))
+        R <- array(0, dim=c(index$N, n.int, ncol(index$defaultprob)))
+        for(t in 1:ncol(index$defaultprob)){
+            S <- 1 - Rstoch[,,t]
+            L[,,t] <- lossdistCZ(p[,t], index$issuerweights, S, index$N, 0, rho, index$Z)
+            R[,,t] <- lossdistCZ(p[,t], index$issuerweights, 1-S, index$N, 0, rho, index$Z)
+        }
+        for(i in 1:n.int){
+            dist <- list(L=L[,i,], R=R[,i,])
+            result[,i] <- MFtranche.pv(index, dist, protection=TRUE)$bp
+        }
+        ## solve the optimization problem
+        program <- KLfit(result[select,], w, index$tranche.quotes[select])
+
+        err <- 0
+        for(i in 1:n.credit){
+            for(j in 1:ncol(p)){
+                err <- err + abs(crossprod(shockprob(p[i,j], rho[i], Z), program$weight) - index$defaultprob[i,j])
+            }
+        }
+
+        ## update the new probabilities
+        p <- MFupdate.probC(index$Z, program$weight, rho, index$defaultprob)
+
+        index$w.mod <- program$weight
+        cat("=")
+    }
+    cat("\n")
+    return(index$w.mod)
+}