library("RQuantLib")
library("parallel")
root = "//WDSENTINEL/share/CorpCDOs"
source(file.path(root, "R", "intex_deals_functions.R"))

MarkitData <- getMarkitIRData()
L1m <- buildMarkitYC(MarkitData, dt = 1/12)
L2m <- buildMarkitYC(MarkitData, dt = 1/6)
L3m <-  buildMarkitYC(MarkitData)
L6m <- buildMarkitYC(MarkitData, dt = 1/2)
L12m <- buildMarkitYC(MarkitData, dt = 1)
setEvaluationDate(as.Date(MarkitData$effectiveasof))
today <- function(){
    return(as.Date(MarkitData$effectiveasof))
}

bps <- 1e-4
global.params <- list()
global.params$recovery.assumptions <- list("Loan"=0.7,
                                           "SecondLien"=0.3,
                                           "Bond"=0.4,
                                           "Mezzanine"=0.15,
                                           "Adj_Covlite"=0.1)

global.params$cdoprices <- list("Aaa"=90,
                                "Aa"=80,
                                "A"=70,
                                "Baa"=60,
                                "Ba"=50,
                                "B"=40,
                                "NR"=40)
#reinvest in 7 years assets
global.params$rollingmaturity <- 7 * 365
global.params$defaultedlag <- 90
global.params$defaultcorr <- 0.4
global.params$defaultbondhazardrate <- 500 * bps
global.params$defaultloanhazardrate <- 400 * bps
global.params$alpha <- 0.25
global.params$beta <- 15
global.params$shape <- function(T)0.25+(1-exp(-T/5))

## dealnames <- c("stonln1", "babs072", "symph4", "flags5", "cent11", "wasatl", "oceant2", "acacl071", "limes")
dealnames <- c("acacl071", "limes")
calibration.date <- "2012-11-08"
calibration <- read.table(file.path(root, "Scenarios", paste0("calibration-",calibration.date,".csv")),
                          sep=",", header=T)
Z <- calibration$Z
w <- calibration$w

rho <- 0.45
Ngrid <- 201
cl <- makeCluster(6)
clusterExport(cl, list("shockprob", "lossdistC.prepay.joint", "lossrecovdist.joint.term", "lossdistC.joint"))
support <- seq(0, 1, length=Ngrid)
for(deal.name in dealnames){
    deal.portfolio <- buildSC.portfolio(deal.name, global.params)
    deal.data <- getdealdata(deal.name)
    A <- SPmatrix2(deal.portfolio$SC, deal.data, freq="3 months")
    S <- 1 - sapply(deal.portfolio$SC, attr, "recov")
    deal.weights <- deal.portfolio$notional/sum(deal.portfolio$notional)
    clusterExport(cl, list("deal.weights"))

    deal.dates <- getdealschedule(deal.data)
    ## compute reinvestment price
    reinvloanprice <- rep(0, length(deal.dates))
    reinvmezzprice <- rep(0, length(deal.dates))
    for(i in 1:length(deal.dates)){
        reinvloanprice[i] <- forwardportfolioprice(deal.portfolio, deal.dates[i], global.params$rollingmaturity, "FLOAT", 0.0235)
        reinvmezzprice[i] <- forwardportfolioprice(deal.portfolio, deal.dates[i], global.params$rollingmaturity, "FLOAT", 0.0385)
    }

    dp <- A$DP
    pp <- A$PP
    Smat <- matrix(S, length(deal.weights), ncol(dp))
    ##no correlation setup
    ##test <- lossrecovdist.joint.term(dp, pp, issuerweights, Smat, Ngrid)

    dpmod <- MFupdate.prob(Z, w, rho, dp)
    ppmod <- MFupdate.prob(-Z, w, rho, pp)
    ## dist <- MFlossrecovdist.prepay(w, Z, rho, dp, dpmod, pp, ppmod, deal.weights, 1-S, Ngrid, TRUE)
    dist.joint <- MFlossdist.prepay.joint(cl, w, Z, rho, dp, dpmod,
                                          pp, ppmod, deal.weights, 1-S, Ngrid, FALSE, n.chunks=4)
    clusterCall(cl, gc)
    gc()
    ## if don't want to use cluster (use less memory)
    ## dist.joint <- MFlossdist.prepay.joint(NULL, w, Z, rho, dp, dpmod,
    ##                                       pp, ppmod, issuerweights, 1-S, Ngrid=201, FALSE)

    ## two ways to compute the joint (D, R) distribution
    ## first using the actual function (numerically instable)
    ## dist.joint2 <- MFlossdist.prepay.joint(cl, w, Z, rho, dp, dpmod,
    ##                                        pp, ppmod, issuerweights, 1-S, Ngrid=201, TRUE)

    ## second, by doing a change of variable seems to work better for now
    distDR <- dist.transform(dist.joint)

    ## compute E(R|D)
    R <- matrix(0, Ngrid, ncol(dp))
    for(t in 1:ncol(dp)){
        R[,t] <- (sweep(distDR[t,,], 1, rowSums(distDR[t,,]), "/") %*% support)/support
    }
    R[1,] <- 0
    n.scenarios <- 100
    percentiles <- (seq(0, 1, length=n.scenarios+1)[-1]+
                    seq(0, 1, length=n.scenarios+1)[-(n.scenarios+1)])/2

    ## compute scenariosd
    scenariosd <- matrix(0, n.scenarios, ncol(dp))
    for(t in 1:ncol(dp)){
        D <- rowSums(distDR[t,,])
        D <- D/sum(D)
        Dfun <- splinefun(c(0, cumsum(D)), c(0, support), "monoH.FC")
        ## dvallow <- floor(Dfun(percentiles)*(Ngrid-1))
        ## dvalup <- ceil(Dfun(percentiles)*(Ngrid-1))
        scenariosd[,t] <- Dfun(percentiles)
    }

    ## compute scenariosr
    scenariosr <- matrix(0, n.scenarios, ncol(dp))
    for(t in 1:ncol(dp)){
        Rfun <- approxfun(support, R[,t], rule=2)
        scenariosr[,t] <- Rfun(scenariosd[,t])
    }

    cdr <- cdrfromscenarios(scenariosd, deal.dates)
    intexrecov <- recoveryfromscenarios(scenariosd, scenariosr)

    ## linear approximation for monthly scenarios
    deal.datesmonthly <- getdealschedule(deal.data, "1 month")
    cdrmonthly <- matrix(0, n.scenarios, length(deal.datesmonthly))
    recoverymonthly <- matrix(0, n.scenarios, length(deal.datesmonthly))
    for(i in 1:n.scenarios){
        cdrmonthly[i,] <- approx(deal.dates, cdr[i,], deal.datesmonthly, rule=2)$y
        recoverymonthly[i,] <- approx(deal.dates, intexrecov[i,], deal.datesmonthly, rule=2)$y
    }

    write.table(cdrmonthly, file=file.path(root, "Scenarios", paste0(deal.name,"-cdr.csv")), row.names=F, col.names=F, sep=",")
    write.table(recoverymonthly * 100, file=file.path(root, "Scenarios", paste0(deal.name,"-recovery.csv")), row.names=F, col.names=F, sep=",")
    cat("generated scenarios for:", deal.name, "\n")
}