library(RQuantLib)
library(statmod)
root = "//WDSENTINEL/share/CorpCDOs/R"
setwd(root)
source(file.path(root, "yieldCurve.R"))
source(file.path(root, "cds_functions_generic.R"))
source(file.path(root, "etdb.R"))
load(file.path(root, "bloomberg_data.RData"))
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.3,
                                           "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))

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(dealname, date=Sys.Date()){
    sqlstring <- sprintf("select * from clo_universe where dealname='%s' order by \"Latest Update\"", dealname)
    return( dbGetQuery(dbCon, sqlstring)[1,] )
}

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

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, startdate, global.params){
    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 <- function(line.item, global.params, startdate=today()){
     ## cat(i, "\n")
    line.item <- collatdata[i,]
    ## cat(line.item$issuername, "\n")
    if( is.na(line.item$maturity) ){
        stop("empty maturity")
    }
    ##most likely equity, doesn't impact the risk anyway
    if(line.item$currentbalance < 1){
        next
    }
    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){
            line.item$currentbalance <- line.item$currentbalance * line.item$price/100
            SC@startdate <- startdate + global.params$defaultedlag
            line.item$maturity <- min(dealdata$maturity, SC@startdate + global.params$rollingmaturity)
            ## automatic reinvest
            SC<- stackcurve(SC, line.item, SC@startdate, global.params)
        }else if( is.na(line.item$price) ){ #missing price
            if(line.item$maturity <= startdate){
                if(startdate<=dealdata$"Reinv End Date"){ #reinvest
                    line.item$maturity <- min(dealdata$maturity, startdate + global.params$rollingmaturity)
                    SC <- stackcurve(SC, line.item, SC@startdate, global.params)
                }else{ #no reinvestment
                    SC@dates <- startdate
                    SC@hazardrates <- 0
                    SC@prepayrates <- 0
                }
                SC <- stackcurve(SC, line.item, SC@startdate, global.params)
            }else{
                SC <- stackcurve(SC, line.item, SC@startdate, global.params)
            }
        }else{
            ## normal case
            if(line.item$maturity > startdate){
                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@curve <- try
                }
            }
        }
    ## if(length(maturity(SC))==0){
    ##     browser()
    ## }
    if(maturity(SC) <= dealdata$"Reinv End Date"){ #we reinvest
        newstartdate <- line.item$maturity
        line.item$maturity <- min(dealdata$maturity, newstartdate + global.params$rollingmaturity)
        SC <- stackcurve(SC, line.item, newstartdate, global.params)
    }
}

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)){
        ## cat(i, "\n")
        line.item <- collatdata[i,]
        ## cat(line.item$issuername, "\n")
        if( is.na(line.item$maturity) ){
            stop("empty maturity")
        }
        #most likely equity, doesn't impact the risk anyway
        if(line.item$currentbalance < 1){
            next
        }
        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)))
        ##expired asset
        ## if(i ==185){
        ##     browser()
        ## }
        ## defaulted asset
        if(!is.na(line.item$defaultedflag) && line.item$defaultedflag){
            line.item$currentbalance <- line.item$currentbalance * line.item$price/100
            SC@startdate <- startdate + global.params$defaultedlag
            line.item$maturity <- min(dealdata$maturity, SC@startdate + global.params$rollingmaturity)
            ## automatic reinvest
            SC<- stackcurve(SC, line.item, SC@startdate, global.params)
        }else if( is.na(line.item$price) ){ #missing price
            if(line.item$maturity <= startdate){
                if(startdate<=dealdata$"Reinv End Date"){ #reinvest
                    line.item$maturity <- min(dealdata$maturity, startdate + global.params$rollingmaturity)
                    SC <- stackcurve(SC, line.item, SC@startdate, global.params)
                }else{ #no reinvestment
                    SC@dates <- startdate
                    SC@hazardrates <- 0
                    SC@prepayrates <- 0
                }
                SC <- stackcurve(SC, line.item, SC@startdate, global.params)
            }else{
                SC <- stackcurve(SC, line.item, SC@startdate, global.params)
            }
        }else{
            ## normal case
            if(line.item$maturity > startdate){
                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@curve <- try
                }
            }
        }
        ## if(length(maturity(SC))==0){
        ##     browser()
        ## }
        if(maturity(SC) <= dealdata$"Reinv End Date"){ #we reinvest
            newstartdate <- line.item$maturity
            line.item$maturity <- min(dealdata$maturity, newstartdate + global.params$rollingmaturity)
            SC <- stackcurve(SC, line.item, newstartdate, global.params)
        }

        notionalvec <- c(notionalvec, line.item$currentbalance)
        SCvec <- c(SCvec, 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) )
}

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


dpc <- new("defaultprepaycurve", dates=cs$dates, hazardrates=rep(0.05,length(cs$dates)),
           prepayrates=rep(0.01, length(cs$dates)))

dc <- new("defaultcurve", dates=cs$dates, hazardrates=rep(0.05,length(cs$dates)))

octagon8.collateral <- getcollateral("octagon8")
octagon8.dealdata <- getdealdata("octagon8")
pomme <- buildSC.portfolio(today(), octagon8.collateral, octagon8.dealdata)

collateral <- octagon8.collateral[1,]
cs <- couponSchedule(collateral$nextpaydate, collateral$maturity,
                     collateral$frequency, collateral$fixedorfloat,
                     collateral$grosscoupon*0.01, collateral$spread*0.01)

k <- function(h, gamma = 15){
    0.25*exp(-gamma * h)
}

test <- buildSC.portfolio("stonln1", global.params)