#!/usr/bin/Rscript
require(methods)
options(warn=2)
args <- commandArgs(trailingOnly=TRUE)

if(.Platform$OS.type == "unix"){
  root.dir <- "/home/share/CorpCDOs"
}else{
  root.dir <- "//WDSENTINEL/share/CorpCDOs"
}

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
if(length(args) >= 1){
  workdate <- as.Date(args[1])
}else{
  workdate <- Sys.Date()
}

#retrieve 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 <- rep(0.45, n.credit)

result <- matrix(0, 4, n.int)
for(l in 1:100){
    Rstoch <- array(0, dim=c(n.credit, n.int, ncol(SurvProb)))
    for(t in 1:ncol(SurvProb)){
        for(i in 1:n.credit){
            Rstoch[i,,t] <- stochasticrecovC(recov[i], 0, Z, w.mod, rho[i], defaultprob[i,t], p[i,t])
        }
    }
    L <- array(0, dim=c(n.int, Ngrid, ncol(defaultprob)))
    R <- array(0, dim=c(n.int, Ngrid, ncol(defaultprob)))
    for(t in 1:ncol(defaultprob)){
        S <- 1 - Rstoch[,,t]
        L[,,t] <- t(lossdistCZ(p[,t], issuerweights, S, Ngrid, 0, rho, Z))
        R[,,t] <- t(lossdistCZ(p[,t], issuerweights, 1-S, Ngrid, 0, rho, Z))
    }
    for(i in 1:n.int){
        result[,i] <- tranche.pvvec(Kmodified, L[i,,], R[i,,], cs)
    }
    ## 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[i], Z), program$weight) - defaultprob[i,j])
        }
    }
    errvec <- c(errvec, err)

    ## update the new probabilities
    p <- MFupdate.probC(Z, program$weight, rho, defaultprob)

    errvec <- c(errvec, err)
    w.mod <- program$weight
    cat(err,"\n")
}

Rstoch <- array(0, dim=c(n.credit, n.int, ncol(SurvProb)))
for(t in 1:ncol(SurvProb)){
    for(i in 1:n.credit){
        Rstoch[i,,t] <- stochasticrecovC(recov[i], 0, Z, w.mod, rho[i], defaultprob[i,t], p[i,t])
    }
}

Lw <- matrix(0, Ngrid, n.int)
Rw <- matrix(0, Ngrid, n.int)
L <- matrix(0, Ngrid, ncol(defaultprob))
R <- matrix(0, Ngrid, ncol(defaultprob))
for(t in 1:ncol(defaultprob)){
    S <- 1 - Rstoch[,,t]
    Lw <- lossdistCZ(p[,t], issuerweights, S, Ngrid, 0, rho, Z)
    Rw <- lossdistCZ(p[,t], issuerweights, 1-S, Ngrid, 0, rho, Z)
    L[,t] <- Lw%*%w.mod
    R[,t] <- Rw%*%w.mod
}

dist <- list(L=L, R=R)

write.table(data.frame(Z=Z, w=w.mod),
            file=file.path(root.dir, "Scenarios", "Calibration",
            paste0("calibration-", workdate, ".csv")),
            col.names=T, row.names=F, sep=",")

save(singlenames.data, hy19, tranche.upf, dist,
     file = file.path(root.dir, "Scenarios", "Calibration",
     paste0("marketdata-", workdate, ".RData")), compress="xz")