library(RQuantLib)
library(statmod)
root = "//WDSENTINEL/share/CorpCDOs/"
source(file.path(root, "R", "yieldCurve.R"))
source(file.path(root, "R", "cds_functions_generic.R"))
source(file.path(root, "R", "etdb.R"))
source(file.path(root, "R", "tranche_functions.R"))
load(file.path(root, "R", "bloomberg_data.RData"))

cdorating <- function(cusip){
    return( sub("[0-9]","", dataMtge[dataMtge$CUSIP %in% cusip,]$RTG_MDY_INITIAL ))
}

getcollateral <- function(dealname, date=Sys.Date()){
    sqlstring <- sprintf("select * from et_aggdealinfo_historical('%s', '%s')", dealname, date)
    collatdata <- dbGetQuery(dbCon, sqlstring)
    return(collatdata)
}

getdealdata <- function(dealnames){
    sqlstring <- sprintf("select * from latest_clo_universe where dealname in ('%s')",
                         paste(dealnames, collapse="','"))
    return( dbGetQuery(dbCon, sqlstring) )
}

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 procedure, not the fastest probably
    r <- NULL
    for(i in 1:length(cusips)){
        sqlstr <- sprintf("select * from dealname_from_cusip('%s')", cusips[i])
        r <- c(r, as.character(dbGetQuery(dbCon, sqlstr)))
    }
    return( r )
}

cusipsfromdealnames <- function(dealnames){
    unlist(strsplit(getdealdata(dealnames)$"Deal Cusip List", ","))
}

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


maturity <- function(creditcurve){
    if(class(creditcurve)=="creditcurve"){
        dates <- creditcurve@curve@dates
        if(length(dates)){
            return( dates[length(dates)] )
        }else{
            return( creditcurve@startdate )
        }
    }else{
        stop("not of class creditcurve")
    }
}

stackcurve <- function(SC, line.item, global.params, startdate){
    newdates <- seq(startdate, line.item$maturity, by="3 months")
    if(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$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(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@curveprepayrates <- 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
        orig.moody <- cdorating(line.item$cusip)
        if(length(orig.moody)==0){
            orig.moody <- "NR"
        }
        line.item$price <- as.numeric(global.params$cdoprices[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 * line.item$price/100
        }else{
            line.item$currentbalance <- line.item$currentbalance * recovery(line.item)
        }
        SC@startdate <- startdate + global.params$defaultedlag
        line.item$maturity <- min(dealmaturity, SC@startdate + global.params$rollingmaturity)
        ## 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)
    }else if(is.na(line.item$price)){ #missing price
        SC <- stackcurve(SC, line.item, global.params, SC@startdate)
    }else{ #normal case
        if(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)
        if(is.null(try)){
            SC <- stackcurve(SC, line.item, global.params, SC@startdate)
        }else{
            SC@curve <- try
        }
    }
    if(maturity(SC) <= reinvdate){ #we reinvest
        newstartdate <- line.item$maturity
        line.item$maturity <- min(dealmaturity, newstartdate + global.params$rollingmaturity)
        SC <- stackcurve(SC, line.item, global.params, newstartdate)
    }
    return( list(SC=SC, notional=line.item$currentbalance) )
}

buildSC.portfolio <- function(dealname, global.params, startdate=today()) {
    dealdata <- getdealdata(dealname)
    collatdata <- getcollateral(dealname)
    notionalvec <- c()
    SCvec <- c()
    betavec <- c()
    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){
            next
        }
        temp <- buildSC(line.item, dealdata$"Reinv End Date", dealdata$maturity, global.params, startdate)
        notionalvec <- c(notionalvec, temp$notional)
        SCvec <- c(SCvec, temp$SC)
        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) )
}

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 it this happens
    intexrecov <- matrix(0, n.scenarios, ncol(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
            }
        }
    }
    return(intexrecov)
}