source("intex_deals_functions.R")
library("RQuantLib")
library("parallel")
## duration <- function(creditcurve){
##     # computes the duration for a hazard rate curve
##     T <- yearFrac(creditcurve@startdate, creditcurve@curve@dates)
##     r <- DiscountCurve(L3m$params, L3m$tsQuotes, T)$forwards
##     h <- creditcurve@curve@hazardrates
##     if(class(creditcurve)=="creditcurve"){
##         T.ext <- c(0, T)
##         s <- (exp(- (h + r) * T.ext[-1]) - exp(- (h + r)* T.ext[-length(T.ext)]))/(h+r)
##     }else{
##         stop("not of class credit curve")
##     }
##     return( - sum(s) )
## }

## dealnames <- c("babs072", "symph4", "flags5", "cent11", "wasatl", "oceant2", "acacl071", "limes")

## for(dealname in dealnames){
##     cat(dealname,"\n")
##     collatdata <- dbGetQuery(dbCon,
##                              paste("select * from et_aggdealinfo_historical(", dealname, ",",
##                                    Sys.Date(),")",
##                                    sep="'"))
##     dealdata <- dbGetQuery(dbCon,
##                            paste("select maturity, \"Reinv End Date\" from clo_universe where dealname='", dealname, "'",sep=""))
##     portfolio <- buildSC.portfolio(today(), collatdata, dealdata)
## }

## dealname <- "symph4"
## maturity18 <- as.Date("2017-06-20")
## probs <- as.numeric(lapply(portfolio$SC, survivalProbability2, maturity18))
## weights <- portfolio$notional/sum(portfolio$notional)
## S <- weights * (1-as.numeric(lapply(portfolio$SC, attr, "recovery")))
## L.ind <- lossdistrib3(1-probs, S, 0.01)
## EL <- crossprod(seq(0,1,0.01), L.ind)

## library(statmod)
## n.int <- 500
## Z <- gauss.quad.prob(n.int, "normal")$nodes
## w <- gauss.quad.prob(n.int, "normal")$weights
## rho <- 0.35

## g <- 1 - probs
## R <- as.numeric(lapply(portfolio$SC, attr, "recovery"))
## r <- matrix(0, 101, 11)
## rho <- seq(0,1,0.1)
## rho[11] <- 0.999
## for(j in 1:11){
##     Lrho <- matrix(0,101,500)
##     for(i in 1:length(Z)){
##         g.shocked <- shockprob(g, rho[j], Z[i])
##         R.shocked <- stochasticrecov.simple(R, 0, Z[i], rho[j], g)
##         Lrho[,i] <- lossdistrib3(g.shocked, (1-R.shocked) * weights, 0.01)
##     }
##     r[,j] <- Lrho%*%w
## }

## R.shocked <- matrix(0,267,n.int)
## for(i in 1:length(Z)){
##     R.shocked [,i] <- stochasticrecov.simple(R, 0, Z[i], rho[j], g)
## }
## extendcurves <- function(portfolio, dealdata){
##     for(i in 1:length(portfolio)){
##         temp <- portfolio$SC[[i]]@curve
##         temp@prepayrates <- c(temp@prepayrates, temp@prepayrates[length(temp@prepayrates)])
##         temp@dates <- c(temp@dates, dealdata$maturity)
##         temp@hazardrates <- c(temp@hazardrates, temp@hazardrates[length(temp@hazardrates)])
##         portfolio$SC[[i]]@curve <- temp
##     }
##     return( portfolio )
## }

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))


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))

stonln1.portfolio <- buildSC.portfolio("stonln1", global.params)
stonln1.data <- getdealdata("stonln1")
A <- SPmatrix2(stonln1.portfolio$SC, stonln1.data, freq="3 months")
S <- 1 - sapply(stonln1.portfolio$SC, attr, "recov")
issuerweights <- stonln1.portfolio$notional/sum(stonln1.portfolio$notional)


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

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

dpmod <- MFupdate.prob(Z, w.mod, rho, dp)
ppmod <- MFupdate.prob(-Z, w.mod, rho, pp)
cl <- makeCluster(6)
clusterExport(cl, list("shockprob", "rho", "Z", "issuerweights", "lossdistC.prepay.joint",
                       "lossrecovdist.joint.term", "lossdistC.joint", "Ngrid"))
dist <- MFlossrecovdist.prepay(w.mod, Z, rho, dp, dpmod, pp, ppmod, issuerweights, 1-S, Ngrid=201, TRUE)
dist.joint <- MFlossdist.prepay.joint(cl, w.mod, 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.mod, 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, dealdates)
intexrecov <- recoveryfromscenarios(scenariosd, scenariosr)

cdrmonthly <- matrix(0, n.scenarios, ncol(dp)*3)
for(i in 1:n.scenarios){
    cdrmonthly[i,] <- cdr[i, as.numeric(matrix(1:ncol(dp), 3, ncol(dp), byrow=T))]
}
recoverymonthly <- matrix(0, n.scenarios, ncol(dp)*3)
for(i in 1:n.scenarios){
    recoverymonthly[i,] <- intexrecov[i, as.numeric(matrix(1:ncol(dp), 3, ncol(dp), byrow=T))]
}
recoverymonthly
write.table(cdrmonthly, file="stonln1cdr.csv", row.names=F, col.names=F, sep=",")
write.table(recoverymonthly * 100, file="stonln1recovery.csv", row.names=F, col.names=F, sep=",")