library("parallel")
setwd("//WDSENTINEL/share/CorpCDOs/R")
source("cds_utils.R")
source("cds_functions_generic.R")
source("index_definitions.R")
source("tranche_functions.R")
source("yieldcurve.R")
source("optimization.R")

cl <- makeCluster(6)

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 HY17
## calibrate the single names curves
singlenames.data <- read.table(file="clipboard", sep="\t", header=T)
nondefaulted <- singlenames.data[!singlenames.data$ticker %in% hy17$defaulted,]
bps <- 1e-4

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

## clusterEvalQ(cl, {setClass("abstractcurve")
##                   setClass("defaultcurve", contains="abstractcurve",
##                            representation(dates="Date", hazardrates="numeric"))
##                   setClass("creditcurve", representation(issuer="character", startdate="Date",
##                                                          recovery="numeric", curve="defaultcurve"))})
## clusterExport(cl, list("nondefaulted", "cdsdates", "cdshazardrate", "today",
##                        "bps", "couponSchedule", "nextIMMDate", "DiscountCurve", "L3m"))
## test <- parSapply(cl, 1:nrow(nondefaulted), parf)

## parf <- function(i){
##     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,5:9])*0.01,
##                          running = rep(nondefaulted$running[i]*bps,5))
##     return( cdshazardrate(quotes, nondefaulted$recovery[i]/100))
## }

hy17portfolio <- 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,5:9])*0.01,
                         running=rep(nondefaulted$running[i]*bps, 5))
    SC@curve <- cdshazardrate(quotes, nondefaulted$recovery[i]/100)
    hy17portfolio <- c(hy17portfolio, SC)
}

issuerweights <- rep(1/length(hy17portfolio), length(hy17portfolio))
hy17$indexref <- 1.02
hy17portfolio.tweaked <- tweakcurves(hy17portfolio, hy17)

SurvProb <- SPmatrix(hy17portfolio.tweaked, hy17)
## load common parameters
K <- c(0, 0.15, 0.25, 0.35, 1)
Kmodified <- adjust.attachments(K, hy17$loss, hy17$factor)
tranche.upf <- c(44.875, 91.75, 104.8125, 114.3125)
tranche.running <- c(0.05, 0.05, 0.05, 0.05)

Ngrid <- 2*nrow(nondefaulted)+1
recov <- sapply(hy17portfolio.tweaked, attr, "recovery")
cs <- couponSchedule(nextIMMDate(today()), hy17$maturity,"Q", "FIXED", 0.05, 0)

## calibrate the tranches using base correlation
rhovec <- c()
f <- function(rho, ...){
    temp <- BClossdistC(SurvProb, issuerweights, recov, rho, Ngrid)
    bp <- 100*(1+1/(Kmodified[i]-Kmodified[i-1]) *
               (tranche.pv(temp$L, temp$R, cs, 0, Kmodified[i], Ngrid) -
                tranche.pv(oldtemp$L, oldtemp$R, cs, 0, Kmodified[i-1], Ngrid)))
    return( abs(tranche.upf[i-1]-bp))
}

for(i in 2:length(Kmodified)){
    rho <- optimize(f, interval=c(0,1),
                    SurvProb, issuerweights, recov, lu, tranche.upf, Kmodified, cs, oldtemp)$minimum
    oldtemp <- BClossdistC(SurvProb, issuerweights, recov, rho, Ngrid)
    rhovec <- c(rhovec, rho)
}

deltas <- c()
for(i in 2:5){
    deltas <- c(deltas, BCtranche.delta(hy17portfolio.tweaked, hy17, 0.05, K[i-1], K[i], rhovec[i-1], rhovec[i], Ngrid))
}

##calibrate by modifying the factor distribution
bottomup <- 1:3
topdown <- 2:4
n.int <- 100
n.credit <- 96
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("shockprob", "issuerweights", "rho", "Z", "lossrecovdist.term",
                       "lossrecovdist", "lossdistribC", "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, 0, 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")
}
clusterExport(cl, list("shockprob", "issuerweights", "rho", "Z", "lossrecovdist.term",
                       "lossrecovdist", "lossdistribC", "Ngrid",
                       "tranche.pvvec", "tranche.pv", "tranche.pl", "tranche.cl",
                       "trancheloss", "trancherecov", "pos", "Kmodified", "cs"))

MFtranche.pv <- function(cl, cs, w, rho, defaultprob, p, issuerweights,
                         Ngrid=length(issuerweights)+1, Kmodified, n.int=100){
    ## computes the tranches pv using the modified factor distribution
    ## p is the modified probability so that
    n.credit <- length(issuerweights)
    Rstoch <- array(0, dim=c(n.int, n.credit, ncol(defaultprob)))
    for(t in 1:ncol(defaultprob)){
        for(i in 1:n.credit){
            Rstoch[,i,t] <- stochasticrecov(recov[i], 0, Z, w, rho, defaultprob[i,t], p[i,t])
        }
    }
    parf <- function(i){
        pshocked <- apply(p, 2, shockprob, rho=rho, Z=Z[i])
        S <- 1 - Rstoch[i,,]
        dist <- lossrecovdist.term(pshocked, 0, issuerweights, S, Ngrid)
        return( tranche.pvvec(Kmodified, dist$L, dist$R, cs))
    }
    clusterExport(cl, list("Rstoch", "p"))
    result <- parSapply(cl, 1:n.int, parf)
    return( 100*(1+result%*%w.mod) )
}

## computes deltas
newportf <- hy17portfolio.tweaked
eps <- 1e-4
for(i in 1:length(newportf)){
    newportf[[i]]@curve@hazardrates <- hy17portfolio.tweaked[[i]]@curve@hazardrates * (1 + eps)
}
SurvProb2 <- SPmatrix(newportf, hy17)
p2 <- MFupdate.prob(Z, w.mod, rho, 1-SurvProb2)
dPVtranches <- MFtranche.pv(cl, cs, w.mod, rho, 1-SurvProb2, p2, issuerweights) - MFtranche.pv(cl, cs, w.mod, rho, defaultprob, p, issuerweights)
dPVindex <- indexpv(newportf, hy17)-indexpv(hy17portfolio.tweaked, hy17)
MFdeltas <- dPVtranches/dPVindex

#global deltas
PVtranches <- MFtranche.pv(cl, cs, w.mod, rho, defaultprob, p, issuerweights)
PVindex <-