library(RQuantLib)
library(statmod)
library(data.table)
options(datatable.nomatch=0)
library(doParallel)

hostname <- system("hostname", intern=TRUE)
if(hostname=="debian"){
    registerDoParallel(8)
}else{
    registerDoParallel(4)
}

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

source(file.path(root.dir, "code", "R", "yieldCurve.R"))
source(file.path(root.dir, "code", "R", "cds_functions_generic.R"))
source(file.path(root.dir, "code", "R", "etdb.R"))
source(file.path(root.dir, "code", "R", "tranche_functions.R"))

getdealdata <- function(dealnames, date){
    if(missing(date)){
        sqlstring <- sprintf("select * from latest_clo_universe where dealname in ('%s')",
                             paste(dealnames, collapse="','"))
    }else{
        sqlstring <- paste("SELECT b.* FROM (SELECT MAX(\"Latest Update\") AS latestdate,",
                           "dealname FROM clo_universe WHERE \"Latest Update\"<='%s' AND dealname in ('%s')",
                           "GROUP BY dealname) a JOIN clo_universe b ON a.dealname = b.dealname",
                           "AND a.latestdate= b.\"Latest Update\" ORDER by dealname asc")
        sqlstring <- sprintf(sqlstring, date, paste(dealnames, collapse="','"))
    }
    data <- dbGetQuery(dbCon, sqlstring)
    data <- data[match(dealnames, data$dealname),]
    rownames(data) <- c()
    return( data )
}

getcollateral <- function(dealname, date){
    if(missing(date)){
        sqlstring <- sprintf("select * from et_aggdealinfo('%s')", dealname)
    }else{
        sqlstring <- sprintf("select * from et_aggdealinfo_historical('%s', '%s')", dealname, date)
    }
    collatdata <- dbGetQuery(dbCon, sqlstring)
    return(collatdata)
}

listdealnames <- function(){
    sqlstring <- "select distinct dealname from clo_universe order by dealname"
    return( dbGetQuery(dbCon, sqlstring))
}

cusip.data <- function(){
    ## TODO: make it date dependent
    sqlstring <- "SELECT a.cusip, b.maturity, a.coupon AS grosscoupon, a.spread,
CASE WHEN a.floater_index like 'LIBOR%' THEN 'FLOAT' ELSE 'FIXED' END
AS fixedorfloat, a.orig_moody from latest_cusip_universe a
LEFT JOIN latest_clo_universe b ON a.dealname = b.dealname"
    data <- dbGetQuery(dbCon, sqlstring)
    data <- data.table(data)
    setkey(data, "cusip")
    return( data )
}

recovery <- function(collateral) {
    ## return assumed recovery based on assumptions from recovery.assumptions
    if(!is.na(collateral$secondlien) && collateral$secondlien){
        collateral$assettype <- "SecondLien"
    }
    recovery <- with(global.params,
                     as.numeric(recovery.assumptions[collateral$assettype]))
    if( !is.na(collateral$covlite) && collateral$covlite) {
        recovery <- recovery - global.params$recovery.assumptions$Adj_Covlite
    }
    if( !is.na(collateral$iscdo) && collateral$iscdo ){
        recovery <- 0
    }
    ## price is too low need to lower the assumed recovery
    if(!is.na(collateral$price) && recovery > collateral$price/100 - 0.1){
        recovery <- max(collateral$price/100-0.2, 0)
    }
    return(recovery)
}

dealnamefromcusip <- function(cusips){
    ## wrapper around the sql function
    sqlstr <- sprintf("select * from dealname_from_cusip('%s')",
                      paste(cusips, collapse="','"))
    r <- dbGetQuery(dbCon, sqlstr)
    return( r$p_dealname )
}

cusipsfromdealnames <- function(dealnames){
    sqlstring <-
        sprintf("select unnest(\"Deal Cusip List\") from latest_clo_universe where dealname in ('%s')",
                paste(dealnames, collapse="','"))
    return( dbGetQuery(dbCon, sqlstring)$unnest )
}

fithazardrate.fast <- function(collateral, eps=1e-6){
    lambda <- 0.05
    cs <- couponSchedule(collateral$nextpaydate, collateral$maturity,
                         collateral$frequency, collateral$fixedorfloat,
                         collateral$grosscoupon * 0.01, collateral$spread*0.01)
    R <- recovery(collateral)
    while(abs(bondprice(lambda, cs, R) * 100 - collateral$price) > eps){
        lambda <- lambda - (bondprice(lambda, cs, R) - 0.01*collateral$price)/dbondprice(lambda, cs, R)
    }
    return( lambda )
}

vanillabondprice <- function(h, collateral, prepay=TRUE) {
    R <- recovery(collateral)
    cs <- couponSchedule(collateral$nextpaydate, collateral$maturity,
                         collateral$frequency, collateral$fixedorfloat,
                         collateral$grosscoupon*0.01, collateral$spread*0.01)
    if(prepay){
        dpc <- new("defaultprepaycurve", dates=cs$dates, hazardrates=rep(h,length(cs$dates)),
                   prepayrates=rep(k(h), length(cs$dates)))
    }else{
        dpc <- new("defaultprepaycurve", dates=cs$dates, hazardrates=rep(h,length(cs$dates)),
                   prepayrates=numeric(0))
    }
    return( bondprice(cs, dpc, R) )
}

dvanillabondprice <- function(hazardrate, collateral) {
    R <- recovery(collateral)
    cs <- couponSchedule(collateral$nextpaydate, collateral$maturity,
                         collateral$frequency, collateral$fixedorfloat,
                         collateral$grosscoupon*0.01, collateral$spread*0.01)
    return( bondprice(hazardrate, cs, R) )
}

fithazardrate <- function(collateral){
    R <- recovery(collateral)
    cs <- couponSchedule(collateral$nextpaydate, collateral$maturity,
                         collateral$frequency, collateral$fixedorfloat,
                         collateral$grosscoupon*0.01, collateral$spread*0.01)
    f <- function(lambda){
        u <- bondprice(lambda, cs, R )
        return( (u * 100-collateral$price)^2 )
    }
    return( optimize(f, c(0,1), tol=1e-6)$minimum )
}

stackcurve <- function(SC, line.item, global.params, startdate){
    if(line.item$nextpaydate> line.item$maturity){
        SC@curve@hazardrates <- 0
        SC@curve@prepayrates <- 0
        SC@curve@dates <- line.item$maturity
        return( SC )
    }
    newdates <- seq(line.item$nextpaydate, line.item$maturity, by="3 months")
    if(newdates[length(newdates)]<line.item$maturity){
        newdates <- c(newdates, line.item$maturity)
    }
    newdates <- c(startdate, newdates[newdates>startdate])
    if(is.na(line.item$assettype) || line.item$assettype=="Loan"){
        hvec <- global.params$shape(yearFrac(today(), newdates[-1])) * global.params$defaultloanhazardrate
        kvec <- global.params$alpha * exp(-global.params$beta * hvec)
    }else if(line.item$assettype=="Bond" || line.item$assettype=="Credit Default Swap"
             || (!is.na(line.item$iscdo) && line.item$iscdo)){
        hvec <- global.params$shape(yearFrac(today(), newdates[-1])) * global.params$defaultbondhazardrate
        kvec <- rep(0, length(hvec))
    }
    SC@curve@hazardrates <- c(SC@curve@hazardrates, hvec)
    SC@curve@prepayrates <- c(SC@curve@prepayrates, kvec)
    SC@curve@dates <- c(SC@curve@dates, newdates[-1])
    return(SC)
}

buildSC.matured <- function(SC, line.item, reinvdate, dealmaturity, global.params, startdate){
    if(!is.na(reinvdate) && startdate<=reinvdate){ #reinvest
        line.item$maturity <- min(dealmaturity, startdate + global.params$rollingmaturity)
        SC <- stackcurve(SC, line.item, global.params, startdate)
    }else{ #no reinvestment
        SC@curve@dates <- startdate
        SC@curve@hazardrates <- 0
        SC@curve@prepayrates <- 0
    }
    return( SC )
}

buildSC <- function(line.item, reinvdate, dealmaturity, global.params, startdate){
    if(!is.na(line.item$iscdo) && line.item$iscdo && is.na(line.item$price)){
        ##we have prices for some cdos e.g. 210795PS3
        if(is.na(line.item$orig) || line.item$orig_moody == "NA" || length(line.item$orig_moody)==0){
            line.item$orig_moody <- "NR"
        }
        line.item$price <- as.numeric(global.params$cdoprices[gsub("\\d", "", line.item$orig_moody)])
    }
    ##build survival curve
    SC <- new("creditcurve", recovery=recovery(line.item), startdate=startdate,
              issuer=line.item$issuername)
    SC@curve <- new("defaultprepaycurve", dates=as.Date(character(0)))
    ## defaulted asset
    if(!is.na(line.item$defaultedflag) && line.item$defaultedflag){
        if(is.na(line.item$price)){
             line.item$currentbalance <- line.item$currentbalance * recovery(line.item)
        }else{
            line.item$currentbalance <- line.item$currentbalance * line.item$price/100
        }
        line.item$price <- 100
        SC@recovery <- 0.7
        SC@startdate <- startdate + global.params$defaultedlag
        if(global.params$reinvflag){#we reinvest recovery assets
            line.item$maturity <- min(dealmaturity, SC@startdate + global.params$rollingmaturity)
            line.item$nextpaydate <- SC@startdate
            ## automatic reinvest
            SC<- stackcurve(SC, line.item, global.params, SC@startdate)
        }
    }else if(line.item$maturity <= startdate){#matured asset
        SC <- buildSC.matured(SC, line.item, reinvdate, dealmaturity, global.params, startdate)
        if(is.na(line.item$price))line.item$price <- 100
    }else if(is.na(line.item$price)){ #missing price
        SC <- stackcurve(SC, line.item, global.params, SC@startdate)
        cs <- couponSchedule(line.item$nextpaydate, line.item$maturity,
                             line.item$frequency, line.item$fixedorfloat,
                             line.item$grosscoupon*0.01, line.item$spread*0.01)
        line.item$price <- bondpv(cs, SC@curve, recovery(line.item)) * 100
    }else{ #normal case
        if(!is.na(line.item$assettype) && line.item$assettype=="Bond"){ #no prepay rate
            alpha <- 0
        }else{
            alpha <- global.params$alpha
        }
        try <- bondhazardrate.shaped(line.item, global.params$shape,
                                     recovery(line.item), alpha, global.params$beta)
        if(is.null(try)){
            SC <- stackcurve(SC, line.item, global.params, SC@startdate)
        }else{
            SC@curve <- try
        }
    }
    if(!is.na(reinvdate) && creditcurve.maturity(SC) <= reinvdate){
        ## if reinvdate is missing, assume no reinvestment
        ## otherwise reinvest
        newstartdate <- line.item$maturity
        line.item$maturity <- min(dealmaturity, newstartdate + global.params$rollingmaturity)
        SC <- stackcurve(SC, line.item, global.params, newstartdate)
    }
    if(is.na(line.item$price)){
        ## TODO
    }
    beta <- if(!is.na(line.item$iscdo) && line.item$iscdo) 1 else global.params$defaultcorr
    return( list(SC=SC, notional=line.item$currentbalance, price = line.item$price, beta = beta) )
}

buildSC.portfolio <- function(dealname, dealdata, cusipdata, global.params, startdate = Sys.Date()) {
    collatdata <- data.table(getcollateral(dealname, startdate))
    setkey(collatdata, "cusip")

    ## replace the cdo fields by bloomberg data
    collatdata[cusipdata[collatdata$cusip],
               `:=`(maturity=i.maturity, fixedorfloat=i.fixedorfloat,
                    spread=i.spread, grosscoupon=i.grosscoupon, orig_moody=i.orig_moody, iscdo=TRUE)]

    portfolio <- foreach(line.item = iter(collatdata, by='row')) %:% {
        when( !is.na(line.item$maturity) && line.item$currentbalance > 1
             && line.item$assettype!="Equity") } %dopar% {
                 buildSC(line.item, dealdata$"Reinv End Date", dealdata$maturity, global.params, startdate)
             }
    missingpricenotional <- sum(collatdata[is.na(price) & maturity>startdate &
                                           (is.na(iscdo)|!iscdo), currentbalance])
    cdonotional <- sum(collatdata[iscdo==TRUE,currentbalance])
    collatbalance <- sum(collatdata[,currentbalance])

    return( list(notional=vapply(portfolio, function(x)x$notional, numeric(1)),
                 beta = vapply(portfolio, function(x)x$beta, numeric(1)),
                 price = vapply(portfolio, function(x)x$price, numeric(1)),
                 SC = lapply(portfolio, function(x)x$SC),
                 stale = missingpricenotional/collatbalance,
                 cdopercentage = cdonotional/collatbalance,
                 collatbalance = collatbalance
                 ) )
}

cdrfromscenarios <- function(scenarios, dates){
    ## compute the forward cdr rates to pass to intex
    ## so that we match the default curves in scenarios
    cdr <- matrix(0, nrow(scenarios), ncol(scenarios))
    for(i in 1:nrow(scenarios)){
        cdr[i,] <- 100*(1-exp(diff(c(0, log(1-scenarios[i,])))))/diff(c(0, yearFrac(today(), dates)))
    }
    return( cdr )
}

recoveryfromscenarios <- function(scenariosd, scenariosr){
    ## compute the forward recovery rate based on the term
    ## structure of recovery scenarios
    ## we run into trouble for very stressed scenarios
    ## this code should cap the scenarios at 0 if this happens
    intexrecov <- matrix(0, n.scenarios, ncol(scenariosr))
    scenariosr <- scenariosr
    for(i in 1:n.scenarios){
        current <- 1
        intexrecov[i,1] <- scenariosr[i,1]
        for(t in 2:ncol(scenariosr)){
            w <- scenariosd[i,current]/scenariosd[i,t]
            ## if(scenariosr[i,t]-w*scenariosr[i,current]>=0){
            ##     intexrecov[i,t] <- (scenariosr[i,t]-w*scenariosr[i,current])/(1-w)
            ##     current <- current+1
            ## }else{
            ##     intexrecov[i,t] <- 0
            ## }
            intexrecov[i,t] <- (scenariosr[i,t]-w*scenariosr[i,current])/(1-w)
            current <- current + 1
        }
    }
    return(intexrecov)
}

recoveryfromscenarios.fast <- function(scenariosr, scenariosd){
    r <- rbind(scenariosr[,1]/scenariosd[,1],
               apply(scenariosr, 1, diff)/apply(scenariosd, 1, diff))
    return( t(r) )
}

severityfromscenarios <- function(scenariosd, scenariosr){
    ## compute the forward recovery rate based on the term
    ## structure of recovery scenarios
    ## we run into trouble for very stressed scenarios
    ## this code should cap the scenarios at 0 if this happens
    intexseverity <- matrix(0, n.scenarios, ncol(scenariosr))
    for(i in 1:n.scenarios){
        current <- 1
        intexseverity[i,1] <- 1-scenariosr[i,1]
        for(t in 2:ncol(scenariosr)){
            w <- scenariosd[i,current]/scenariosd[i,t]
            intexseverity[i,t] <- 1 - (scenariosr[i,t]-w*scenariosr[i,current])/(1-w)
            current <- current+1
        }
    }
    return(intexseverity)
}

load.portfolio <- function(dealname){
    load(file.path(root.dir, "Scenarios", paste("Portfolios", calibration.date, sep="_"),
                   paste0(dealname, ".RData")), .GlobalEnv)
}

dealupdatedates <- function(dealnames){
    sqlstring <- sprintf("select dealname, \"Latest Update\" from latest_clo_universe where dealname in ('%s')",
                         paste(dealnames, collapse="','"))
    data <- dbGetQuery(dbCon, sqlstring)
    data <- data[match(dealnames, data$dealname),]
    rownames(data) <- c()
    return ( data )
}