1 files changed, 131 insertions, 0 deletions
diff --git a/R/script_calibrate_tranches.R b/R/script_calibrate_tranches.R
new file mode 100644
index 00000000..60578c95
--- /dev/null
+++ b/R/script_calibrate_tranches.R
@@ -0,0 +1,131 @@
+#!/usr/bin/Rscript
+
+library("parallel")
+
+if(.Platform$OS.type == "unix"){
+  root.dir <- "/home/share/CorpCDOs"
+  cl <- makeForkCluster(4)
+}else{
+  root.dir <- "//WDSENTINEL/share/CorpCDOs"
+  cl <- makeCluster(6)
+}
+
+source(file.path(root.dir, "code", "R", "cds_utils.R"))
+source(file.path(root.dir, "code", "R", "cds_functions_generic.R"))
+source(file.path(root.dir, "code", "R", "index_definitions.R"))
+source(file.path(root.dir, "code", "R", "tranche_functions.R"))
+source(file.path(root.dir, "code", "R", "yieldcurve.R"))
+source(file.path(root.dir, "code", "R", "optimization.R"))
+
+MarkitData <- getMarkitIRData()
+L1m <- buildMarkitYC(MarkitData, dt = 1/12)
+L2m <- buildMarkitYC(MarkitData, dt = 1/6)
+L3m <-  buildMarkitYC(MarkitData)
+L6m <- buildMarkitYC(MarkitData, dt = 1/2)
+setEvaluationDate(as.Date(MarkitData$effectiveasof))
+
+## calibrate HY19
+## calibrate the single names curves
+singlenames.data <- read.csv(file.path(root.dir, "Scenarios", "Calibration",
+                                         paste0("hy19_singlenames_", Sys.Date(), ".csv")))
+nondefaulted <- singlenames.data[!singlenames.data$ticker %in% hy19$defaulted,]
+bps <- 1e-4
+
+cdsdates <- as.Date(character(0))
+for(tenor in paste0(1:5, "y")){
+    cdsdates <- c(cdsdates, cdsMaturity(tenor))
+}
+
+hy19portfolio <- c()
+for(i in 1:nrow(nondefaulted)){
+    SC <- new("creditcurve",
+              recovery=nondefaulted$recovery[i]/100,
+              startdate=today(),
+              issuer=as.character(nondefaulted$ticker[i]))
+    quotes <- data.frame(maturity=cdsdates, upfront = as.numeric(nondefaulted[i,4:8]) /100,
+                         running=rep(nondefaulted$running[i] * bps, 5))
+    SC@curve <- cdshazardrate(quotes, nondefaulted$recovery[i]/100)
+    hy19portfolio <- c(hy19portfolio, SC)
+}
+
+issuerweights <- rep(1/length(hy19portfolio), length(hy19portfolio))
+
+## load tranche data
+K <- c(0, 0.15, 0.25, 0.35, 1)
+Kmodified <- adjust.attachments(K, hy19$loss, hy19$factor)
+markit.data <- read.csv(file.path(root.dir, "Scenarios", "Calibration",
+                                  paste0("hy19_tranches_", Sys.Date(), ".csv")))
+
+tranche.upf <- markit.data$Mid
+tranche.running <- c(0.05, 0.05, 0.05, 0.05)
+# get the index ref
+hy19$indexref <- markit.data$bidRefPrice[1]/100
+hy19portfolio.tweaked <- tweakcurves(hy19portfolio, hy19)
+SurvProb <- SPmatrix(hy19portfolio.tweaked, hy19)
+
+Ngrid <- 2 * nrow(nondefaulted) + 1
+recov <- sapply(hy19portfolio.tweaked, attr, "recovery")
+cs <- couponSchedule(nextIMMDate(today()), hy19$maturity,"Q", "FIXED", 0.05, 0)
+
+##calibrate by modifying the factor distribution
+bottomup <- 1:3
+topdown <- 2:4
+n.int <- 500
+n.credit <- length(hy19portfolio)
+errvec <- c()
+quadrature <- gauss.quad.prob(n.int, "normal")
+w <- quadrature$weights
+Z <- quadrature$nodes
+w.mod <- w
+defaultprob <- 1 - SurvProb
+p <- defaultprob
+rho <- 0.45
+
+clusterExport(cl, list("root.dir", "shockprob", "issuerweights", "rho", "Z", "lossrecovdist.term",
+                       "lossrecovdist", "lossdistC", "Ngrid",
+                       "tranche.pvvec", "tranche.pv", "tranche.pl", "tranche.cl",
+                       "trancheloss", "trancherecov", "pos", "Kmodified", "cs"))
+
+## TODO: investigate if this is the right thing w.r.t recovery
+parf <- function(i){
+    pshocked <- apply(p, 2, shockprob, rho=rho, Z=Z[i])
+    S <- 1 - Rstoch[i,,]
+    dist <- lossrecovdist.term(pshocked, , issuerweights, S, Ngrid)
+    return( tranche.pvvec(Kmodified, dist$L, dist$R, cs))
+}
+
+for(l in 1:100){
+    Rstoch <- array(0, dim=c(n.int, n.credit, ncol(SurvProb)))
+    for(t in 1:ncol(SurvProb)){
+        for(i in 1:n.credit){
+            Rstoch[,i,t] <- stochasticrecov(recov[i], 0, Z, w.mod, rho, defaultprob[i,t], p[i,t])
+        }
+    }
+
+    clusterExport(cl, list("Rstoch", "p"))
+    result <- parSapply(cl, 1:n.int, parf)
+    ## solve the optimization problem
+    program <- KLfit(100*(result[bottomup,]+1), w, tranche.upf[bottomup])
+
+    err <- 0
+    for(i in 1:n.credit){
+        for(j in 1:ncol(p)){
+            err <- err + abs(crossprod(shockprob(p[i,j], rho, Z), program$weight) - defaultprob[i,j])
+        }
+    }
+    errvec <- c(errvec, err)
+
+    ## update the new probabilities
+    p <- MFupdate.prob(Z, program$weight, rho, defaultprob)
+
+    errvec <- c(errvec, err)
+    w.mod <- program$weight
+    cat(err,"\n")
+}
+
+write.table(data.frame(Z=Z, w=w.mod),
+            file=file.path(root.dir, "Scenarios", "Calibration",
+            paste0("calibration-", Sys.Date(), ".csv")),
+            col.names=T, row.names=F, sep=",")
+save(singlenames.data, hy19, tranche.upf,
+     file = file.path(root.dir, "Scenarios", "Calibration", paste0("marketdata-", Sys.Date(), ".RData")))