library(RQuantLib)
library(statmod)
library(hash)

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 <- 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 cusip_universe a LEFT JOIN latest_clo_universe b ON a.dealname = b.dealname"
    data <- dbGetQuery(dbCon, sqlstring)
    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){
    ## cat(line.item$issuername, "\n")
    if(!is.na(line.item$iscdo) && line.item$iscdo && is.na(line.item$price)){
        ##we have prices for some cdos e.g. 210795PS3
        if(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@startdate <- startdate + global.params$defaultedlag
        if(!is.na(reinvdate) && SC@startdate<reinvdate && global.params$reinvflag){
            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
    }
    return( list(SC=SC, notional=line.item$currentbalance, price = line.item$price) )
}

buildSC.portfolio <- function(dealname, dealdata, cusipdata, global.params, startdate = Sys.Date()) {
    collatdata <- getcollateral(dealname, startdate)
    index <- hash(cusipdata$cusip, 1:length(cusipdata$cusip))
    notionalvec <- c()
    SCvec <- c()
    betavec <- c()
    pricevec <- c()
    fields <- c("maturity", "fixedorfloat", "spread", "grosscoupon", "orig_moody")
    cdonotional <- 0
    missingpricenotional <- 0
    totalnotional <- 0
    for(i in 1:nrow(collatdata)){
        line.item <- collatdata[i,]
        if( is.na(line.item$maturity) ){
            stop("empty maturity")
        }
        ##most likely equity, doesn't impact the risk anyway
        if(line.item$currentbalance < 1 || line.item$assettype=="Equity"){
            next
        }
        ## overwrite the data with cusipdata info if not NA,
        ## because it should be higher quality
        ## e.g all cdos in acacl061 had wrong maturity
        if(has.key(line.item$cusip, index)){
            line.item$iscdo <- TRUE
            for(field in fields){
                if(!is.na(cusipdata[index[[line.item$cusip]],field])){
                    line.item[field] <- cusipdata[index[[line.item$cusip]],field]
                }
            }
        }
        totalnotional <- totalnotional + line.item$currentbalance
        if(!is.na(line.item$iscdo) && line.item$iscdo){
            cdonotional <- cdonotional + line.item$currentbalance
        }else{
            if(is.na(line.item$price) && line.item$maturity>startdate){
                missingpricenotional <- missingpricenotional + line.item$currentbalance
            }
        }
        temp <- buildSC(line.item, dealdata$"Reinv End Date", dealdata$maturity, global.params, startdate)
        notionalvec <- c(notionalvec, temp$notional)
        SCvec <- c(SCvec, temp$SC)
        pricevec <- c(pricevec, temp$price)
        betavec <- c(betavec, if(!is.na(line.item$iscdo) && line.item$iscdo) 1 else
                     global.params$defaultcorr)
    }
    return( list(notional=notionalvec, SC=SCvec, beta=betavec, price = pricevec,
                 cdopercentage = cdonotional/totalnotional,
                 stale = missingpricenotional/totalnotional, collatbalance = totalnotional) )
}

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