#!/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", "tranche_functions.R"))
source(file.path(root.dir, "code", "R", "yieldcurve.R"))
source(file.path(root.dir, "code", "R", "optimization.R"))
load.index("hy21")

##figure out the tradedate
if(length(args) >= 1){
  tradedate <- as.Date(args[1])
}else{
  tradedate <- Sys.Date()
}

exportYC(tradedate)
## calibrate HY21
## calibrate the single names curves
singlenames.data <- read.csv(file.path(root.dir, "Scenarios", "Calibration",
                                       paste0("hy21_singlenames_", tradedate, ".csv")))
nondefaulted <- singlenames.data[!singlenames.data$ticker %in% hy21$defaulted,]
bps <- 1e-4

cdsdates <- as.Date(character(0))
for(tenor in paste0(1:5, "y")){
    cdsdates <- c(cdsdates, cdsMaturity(tenor, date=tradedate))
}

hy21portfolio <- c()
cs <- couponSchedule(IMMDate(tradedate), cdsdates[length(cdsdates)], "Q", "FIXED",
                     1, tradedate, IMMDate(tradedate, "prev"))
for(i in 1:nrow(nondefaulted)){
    SC <- new("creditcurve",
              recovery=nondefaulted$recovery[i]/100,
              startdate=tradedate,
              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, tradedate, cs)
    hy21portfolio <- c(hy21portfolio, SC)
}

issuerweights <- rep(1/length(hy21portfolio), length(hy21portfolio))

## load tranche data
K <- c(0, 0.15, 0.25, 0.35, 1)
Kmodified <- adjust.attachments(K, hy21$loss, hy21$factor)
market.data <- read.csv(file.path(root.dir, "Scenarios", "Calibration",
                                  paste0("hy21_tranches_", tradedate, ".csv")))

tranche.upf <- market.data$Mid
tranche.running <- c(0.05, 0.05, 0.05, 0.05)
# get the index ref
hy21$indexref <- market.data$bidRefPrice[1]/100
hy21portfolio.tweaked <- tweakcurves(hy21portfolio, hy21, tradedate)$portfolio
cs <- couponSchedule(IMMDate(tradedate), hy21$maturity,"Q", "FIXED", 0.05, 0, tradedate, IMMDate(tradedate, "prev"))
SurvProb <- SPmatrix(hy21portfolio.tweaked, length(cs$dates))

Ngrid <- 2 * nrow(nondefaulted) + 1
recov <- sapply(hy21portfolio.tweaked, attr, "recovery")
acc <- cdsAccrued(tradedate, 0.05)
##calibrate by modifying the factor distribution
bottomup <- 1:3
topdown <- 2:4
n.int <- 500
n.credit <- length(hy21portfolio)
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)
err <- Inf
while(err >0.01){
    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) - acc
    }
    ## 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-", tradedate, ".csv")),
            col.names=T, row.names=F, sep=",")

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