#!/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"))

#figure out the workdate
filenames <- list.files(file.path(root.dir, "Scenarios", "Calibration"),
                               pattern = "hy19_singlenames*")

dates <- strtrim(sapply(strsplit(filenames, "_"), function(x)x[3]), 10)
workdate <- as.Date(sort(dates, decreasing = TRUE)[1])

#retreive yield curve data
MarkitData <- getMarkitIRData(workdate)
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_", workdate, ".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_", workdate, ".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(workdate), 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-", workdate, ".RData")))