library(RQuantLib)
library(data.table)
library(doParallel)
library(lossdistrib)

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

source("cds_functions_generic.R")
source("db.R")
etdb <- dbConn("ET")

getdealdata <- function(dealname, workdate){
    sqlstring <- paste0("select marketvalue from et_deal_model_numbers where dealname=$1 and ",
                        "updatedate in (select max(updatedate) from et_deal_model_numbers where ",
                        "dealname = $2 and updatedate<=$3)")
    mv <- dbGetQuery(etdb, sqlstring,
                     params = list(dealname, dealname, workdate))$marketvalue

    sqlstring <- paste0("select \"Curr Collat Bal\", reinv_end_date, ",
                        "first_pay_date , maturity, \"Principal Bal\" , pay_day from ",
                        "historical_clo_universe($1, $2)")
    dealdata <- dbGetQuery(etdb, sqlstring, params=list(dealname, workdate))
    if(!length(mv)){
        dealdata$mv <- NA
    }else{
        dealdata$mv <- mv
    }
    return(dealdata)
}

getcollateral <- function(dealname, date){
    if(missing(date)){
        collatdata <- dbGetQuery(etdb, "select * from et_aggdealinfo($1)",
                                 params=list(dealname))
    }else{
        collatdata <-  dbGetQuery(etdb,
                                  "select * from et_aggdealinfo_historical($1, $2)",
                                  params=list(dealname, date))
    }
    return(collatdata)
}

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

cusip.data <- function(workdate = Sys.Date()){
    sqlstring <- "SELECT DISTINCT ON (cusip) cusip, maturity, coupon AS grosscoupon,
spread, CASE WHEN floater_index like 'LIBOR%' THEN 'FLOAT' ELSE 'FIXED' END
AS fixedorfloat, orig_moody FROM cusip_universe JOIN deal_indicative USING (dealname)
WHERE updatedate<=$1 ORDER BY cusip, updatedate DESC"
    data <- dbGetQuery(etdb, sqlstring, workdate)
    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(etdb, sqlstr)
    return( r$p_dealname )
}

cusipsfromdealnames <- function(dealnames){
    sqlstring <-
        sprintf("select unnest(\"Deal Cusip List\") from deal_indicative where dealname in ('%s')",
                paste(dealnames, collapse="','"))
    return( dbGetQuery(etdb, 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{
            SC <- NULL
        }

    }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, startdate)
        line.item$price <- bondpv(cs, SC@curve, recovery(line.item), startdate) * 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, startdate)
        if(is.null(try)){
            SC <- stackcurve(SC, line.item, global.params, SC@startdate)
        }else{
            SC@curve <- try
        }
    }
    if(!is.na(reinvdate) && !is.null(SC) && 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,
               `:=`(maturity=i.maturity, fixedorfloat=i.fixedorfloat,
                    spread=i.spread, grosscoupon=i.grosscoupon, orig_moody=i.orig_moody, iscdo=TRUE),
               allow.cartesian=TRUE]

    portfolio <- foreach(line.item = iter(collatdata, by='row')) %:% {
        when( !is.na(line.item$maturity) && line.item$currentbalance > 1
             && !is.na(line.item$assettype) && line.item$assettype!="Equity") } %dopar% {
                 buildSC(line.item, dealdata$reinv_end_date, dealdata$maturity, global.params, startdate)
             }
    ## non-parallel version for debugging
    ## portfolio <- c()
    ## for(i in 1:nrow(collatdata)){
    ##     line.item <- collatdata[i,]
    ##     if(is.na(line.item$maturity) || line.item$currentbalance <= 1){
    ##         next
    ##     }
    ##     portfolio <- c(portfolio, 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[!is.na(iscdo)&(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, tradedate){
    ## 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(tradedate, 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)
}

get.reinvassets <- function(dealname, tradedate){
  r <- list()
  sqlstr <- "select * from et_historicaldealinfo($1, $2) where ReinvFlag Is true"
  data <- dbGetQuery(etdb, sqlstr, params=list(dealname, tradedate))
  if(nrow(data)>0){
      for(i in 1:nrow(data)){
          r[[data$issuername[i]]] <- list(coupontype=data$fixedorfloat[i], liborfloor=data$liborfloor[i])
      }
  }
  return( r )
}

getpayday <- function(dealdata, tradedate){
    ## try to compute the previous pay date of a deal
    ## it relies on thwo things to be accurate: the pay_day
    ## as well as the first_pay_date (that's how we get the month)
    m <- as.numeric(format(dealdata$first_pay_date, "%m"))
    m <- m %%3+9
    y <- as.numeric(format(tradedate, "%Y"))
    y <- y - 1
    payday <- as.Date(sprintf("%s-%s-%s", y, m, dealdata$pay_day))
    i <- 1
    cal <- Calendar$new("UnitedStates")
    nextdate <- cal$advance(dates=payday, period = "3m", bdc = "Unadjusted")
    while(nextdate < tradedate){
        payday <- nextdate
        i <- i+1
        nextdate <- cal$advance(payday, 3*i, 2, bdc = "Unadjusted")
    }
    return(payday)
}

getdealschedule <- function(dealdata, freq = c("Monthly", "Quarterly"), tradedate=Sys.Date(),
                            bdc = c("Unadjusted", "Following", "ModifiedFollowing")) {
    payday <- getpayday(dealdata, tradedate)
    freq <- match.arg(freq)
    bdc <- match.arg(bdc)
    params <- list(effectiveDate = getpayday(dealdata, tradedate),
                   maturityDate = dealdata$maturity,
                   period = freq,
                   businessDayConvention = bdc,
                   terminationDateConvention = "Unadjusted",
                   dateGeneration = "Forward")
    return( Schedule(params) )
}

intexportfolio.forwardprice <- function(cdrmonthly, recoverymonthly, startdate, maturity,
                                   coupontype, margin, liborfloor){
    if(missing(liborfloor)||is.na(liborfloor)){
        currentcoupon <- margin
    }else{
        currentcoupon <- margin + liborfloor
    }
    forwardcs <- data.table(couponSchedule(nextpaydate=startdate+45, maturity,
                                           frequency="Q", coupontype, margin,
                                           currentcoupon, tradedate=startdate), key="dates")
    notionals <-  cdrmonthly[date>=startdate, lapply(.SD,function(x)cumprod(1-x/100*1/12)),
                             .SDcols=paste0("V",1:100)]
    recovery <- as.matrix(recoverymonthly[date>=startdate, .SD, .SDcols=paste0("V",1:100)])*
        -apply(rbind(1,as.matrix(notionals)), 2, diff)
    if(nrow(recovery)==1){
        recovery <- recovery*last(forwardcs[,df])
    }else{
        recovery <- data.table(dates=cdrmonthly[date>=startdate,date],apply(recovery, 2, cumsum),key="dates")
        recovery <- recovery[forwardcs, roll=TRUE]
        df <- recovery[,df]
        recovery <- t(df)%*%as.matrix(recovery[,lapply(.SD,function(x)diff(c(0,x))),.SDcols=paste0("V",1:100)])
    }
    notionals <- data.table(dates=cdrmonthly[date>=startdate,date], notionals, key="dates")
    outstanding <- notionals[forwardcs, roll=TRUE]
    mat.outstanding <- as.matrix(outstanding[,.SD,.SDcols=paste0("V",1:100)])
    po <- mat.outstanding[nrow(mat.outstanding),]*last(outstanding)[,df]
    io <- outstanding[, df*coupons]%*%mat.outstanding
    mean(recovery+po+io)
}

compute.reinvprices <- function(dealname, cdrmonthly, recoverymonthly, params, rollingmaturity, tradedate){
    reinvassets <- get.reinvassets(dealname, tradedate)
    reinvprices <- list()
    if(length(reinvassets)>0){
        maturity <- cdrmonthly$date[nrow(cdrmonthly)]
        for(assetname in names(reinvassets)){
            asset <- reinvassets[[assetname]]
            coupon <- if(asset$coupontype=="FLOAT") params$reinvfloat else params$reinvfixed
            reinvprices[[assetname]] <- foreach(date = iter(cdrmonthly$date), .combine=c) %dopar% {
                100 * intexportfolio.forwardprice(cdrmonthly.dt, recoverymonthly.dt, date,
                                                  min(date+rollingmaturity, maturity),
                                                  asset$coupontype, coupon, asset$liborfloor/100)
            }
        }
    }
    return(reinvprices)
}