source("cds_utils.R")

setClass("defaultcurve", representation(dates="Date", hazardrates="numeric"))
setClass("defaultprepaycurve", representation(prepayrates="numeric"), contains="defaultcurve")
setClass("creditcurve", representation(issuer="character", startdate="Date",
                                       recovery="numeric", curve="defaultprepaycurve"))

survivalProbability1 <- function(startdate, date, survival.curve) {
    #based on a flat hazard rate curve
    T <- yearFrac(startdate, survival.curve$dates)
    Tmat <- yearFrac(startdate, date)
    for ( i in 1:length(survival.curve$dates) ){
        if ( date >= survival.curve$dates[i] ) {
            next
        }else{
            if( i > 1 ) {
                w <- ( Tmat - T[i-1] ) / (T[i] - T[i-1])
                logprob <- - (1-w) * survival.curve$hazardrates[i-1] * T[i-1] -
                    w * survival.curve$hazardrates[i] * T[i]
            }else{
                logprob <- - Tmat * survival.curve$hazardrates[1]
                return( exp(as.numeric(logprob)) )
            }
        }
    }
    ## if date is greater than last survival.curve date, keep the hazard rate flate
    logprob <- - yearFrac(startdate, date) * survival.curve$hazardrates[i]
    return( exp(as.numeric(logprob)) )
}

survivalProbability.exact <- function(credit.curve, date) {
    #based on a forward hazard rate curve
    curve <- credit.curve@curve
    T <- c(0, yearFrac(credit.curve@startdate, curve@dates))
    dT <- diff(T)
    Tmat <- yearFrac(credit.curve@startdate, date)
    logprob <- 0
    for ( i in 1:length(dT) ){
        if ( date > curve@dates[i] ) {
            logprob <- logprob - curve@hazardrates[i] * dT[i]
        }else{
            if( i > 1 ){
                logprob <- logprob - curve@hazardrates[i] * (Tmat - T[i])
            }else{
                logprob <- logprob - curve@hazardrates[1] * Tmat
            }
            break
        }
    }
    return( exp(as.numeric(logprob)) )
}

PD <- function(sc){
    ## computes the default probability associated with the survival curve
    if(class(sc)=="defaultprepaycurve"){
        T <- c(0, yearFrac(today(), sc@dates))
        dT <- diff(T)
        return( 1-cumprod(exp(-sc@hazardrates * dT)) )
    }else{
        stop("not a default curve")
    }
}

couponleg <- function(cs, sc, accruedondefault=TRUE){
    ## computes the pv of the risky coupon leg based on a given coupon schedule
    ## and a survival curve. Also called premium leg or fixed leg.
    dT <- diff(c(0, yearFrac(today(), cs$dates)))
    PD <- cumprod(exp(-sc@hazardrates[1:length(dT)] * dT))
    if(accruedondefault){
        PDadj <- 0.5 * (c(1, PD[-length(PD)]) + PD)
    }else{
        PDadj <- PD
    }
    return( crossprod(PDadj, cs$coupons * cs$df) )
}

couponleg.flat <- function(cs, h, accruedondefault=TRUE){
    T <- yearFrac(today(), cs$dates)
    PD <- exp(- h * T )
    if(accruedondefault){
        PDadj <- 0.5 * (c(1, PD[-length(PD)]) + PD)
    }else{
        PDadj <- PD
    }
    return( crossprod(PDadj, cs$coupons * cs$df) )
}

couponleg.prepay <- function(cs, dpc, accruedondefault=TRUE){
    ## computes the pv of the risky coupon leg from the coupon schedule,
    ## a hazard rate curve, and a prepay curve. We assume the poisson process driving
    ## default and prepayment are independent, so the intensity of either event
    ## happenning is the sum of the two.

    dT <- diff(c(0, yearFrac(today(), cs$dates)))
    SP <- cumprod(exp( - (dpc@hazardrates[1:length(dT)] + dpc@prepayrates[1:length(dT)] ) * dT))
    if(accruedondefault){
        SPadj <- 0.5 * (c(1, SP[-length(SP)]) + SP)
    }else{
        SPadj <- SP
    }
    return( crossprod(SPadj, cs$coupons * cs$df) )
}

couponleg.prepay.flat <- function(cs, h){
     dpc <- new("defaultprepaycurve", dates = cs$dates, hazardrates=rep(h, length(cs$dates)),
                prepayrates=rep(k(h), length(cs$dates)))
     couponleg.prepay(cs, dpc)
}

dcouponleg <- function(cs, sc, index, accruedondefault=TRUE){
    ## derivative of couponleg with respect to hazardrate
    dT <- diff(c(0, yearFrac(today(), cs$dates)))
    PD <- cumprod(exp(-sc@hazardrates[1:length(dT)] * dT))
    dPD <- cumsum(index * dT) * PD
    if(accruedondefault){
        dPDadj <- 0.5 * (c(0, dPD[-length(PD)]) + dPD)
    }else{
        dPDadj <- dPD
    }
    return( - crossprod( dPDadj, cs$coupons * cs$df) )
}

k <- function(h, alpha=0.25, beta=15){
    ## prepay rate as a function of the hazardrate
    ## this is a decreasing function
    ## alpha is the maximum prepay rate
    return ( alpha * exp(- beta * h) )
}

dcouponleg.prepay <- function(cs, dpc, index, beta, accruedondefault=TRUE){
    ## derivative of couponleg.prepay with respect to hazardrate
    ## If k is the prepay rate, it assumes dk/dh = - beta * k,
    ## which is the case if k(h) = alpha * exp(-beta *h)
    dT <- diff(c(0, yearFrac(today(), cs$dates)))
    SP <- cumprod(exp( - (dpc@hazardrates[1:length(dT)] + dpc@prepayrates[1:length(dT)] ) * dT))
    dSP <- -cumsum(index * dT * (1 - beta * dpc@prepayrates[1:length(dT)])) * SP
    if(accruedondefault){
        dSPadj <- 0.5 * (c(0, dSP[-length(SP)]) + dSP)
    }else{
        SPadj <- dSP
    }
    return( crossprod(dSPadj, cs$coupons * cs$df) )
}

dcouponleg.prepay.flat <- function(cs, h, index, beta, accruedondefault=TRUE){
    hvec <- rep(h, length(cs$dates))
    kvec <- k(h)
    dpc <- new("defaultprepaycurve", dates=cs$dates, hazardrates=hvec
               prepayrates=kvec )
    return( dcouponleg.prepay(cs, dpc, index, beta, accruedondefault))
}

cdsduration <- function(sc, maturity,  accruedondefault=TRUE){
    # computes the risky PV01, also called risky annuity of a cds
    cs <- couponSchedule(nextIMMDate(today()), maturity, "Q", "FIXED", 1)
    couponleg(cs, sc, accruedondefault)
}

defaultleg.flat <- function(cs, h, recovery){
    ## Computes the pv of the default leg of a cds based on a given
    ## coupon schedule, flat hazard rate, and recovery.
    T <- yearFrac(today(), cs$dates)
    Q <- exp(-h * T) * cs$df
    Qmid <- 1/2 * (c(1, Q[-length(Q)]) + Q)
    r <- (1 - recovery) * crossprod(h * Qmid, dT)
    return( r )
}

ddefaultleg.flat <- function(cs, h, recovery){
    ## derivative of defaultleg.flat with respect to hazardrate
    T <- yearFrac(today(), cs$dates)
    dT <- diff(c(0, T))
    dQ <- - T * exp(-h * T) * cs$df
    Q <- exp(-h * T) * cs$df
    Qmid <- 1/2 * (c(1, Q[-length(Q)]) + Q)
    dQmid <- 1/2 * (c(0, dQ[-length(dQ)]) + dQ)
    dr <- (1-recovery) * (crossprod(Qmid, dT) + h * crossprod(dQmid, dT))
    return( dr )
}

defaultleg <- function(cs, sc, recovery){
    ## Computes the pv of the default leg of a cds based on a given
    ## coupon schedule, hazard rate curve, and recovery.
    T <- yearFrac(today(), cs$dates)
    dT <- diff(c(0, T))
    Q <- cumprod(exp(-sc@hazardrates[1:length(dT)] * dT)) * cs$df
    Qmid <- 1/2 * (c(1, Q[-length(Q)]) + Q)
    r <- (1 - recovery) * crossprod(sc@hazardrates[1:length(dT)] * Qmid, dT)
    return( r )
}

ddefaultleg <- function(cs, sc, recovery, index){
    ## derivative of defaultleg with respect to hazardrate
    T <- yearFrac(today(), cs$dates)
    dT <- diff(c(0,T))
    Q <- cumprod(exp(-sc@hazardrates[1:length(dT)] * dT)) * cs$df
    dQ <- - cumsum(index * dT) * Q
    Qmid <- 1/2 * (c(1, Q[-length(Q)]) + Q)
    dQmid <- 1/2 *(c(0, dQ[-length(dQ)]) + dQ)
    dr <- (1-recovery) * (crossprod(index * Qmid, dT) + crossprod(sc@hazardrates * dQmid, dT))
    return( dr )
}

defaultleg.prepay <- function(cs, sc, recovery){
    ## Computes the pv of the default leg of a cds based on a given
    ## coupon schedule, hazard rates curve, prepay curves, and recovery.
    T <- yearFrac(today(), cs$dates)
    dT <- diff(c(0, T))
    Q <- cumprod(exp(- (dpc@hazardrates[1:length(dT)]+dpc@prepayrates[1:length(dT)]) * dT)) * cs$df
    Qmid <- 1/2 * (c(1, Q[-length(Q)]) + Q)
    r <- (1 - recovery) * crossprod(dpc@hazardrates[1:length(dT)] * Qmid, dT)
    +crossprod(dpc@prepayrates[1:length(dT)] * Qmid, dT)
    return( r )
}

cdspv <- function(cs, sc, R){
    return ( couponleg(cs, sc) - defaultleg(cs, sc, R) )
}

cdsspread <- function(sc, maturity, recovery){
    ## computes exact cds running spread for a cds
    ## should be very close to hazardrate * (1-recovery)
    cs <- couponSchedule(nextIMMDate(today()), maturity, "Q", "FIXED", 1)
    defaultleg(cs, sc, recovery)/couponleg(cs, sc)
}

dcdspv <- function(cs, sc, recovery, index){
    ## derivative of cdspv with respect to hazardrate
    return ( dcouponleg(cs, sc, index) - ddefaultleg(cs, sc, recovery, index) )
}

cdshazardrate.flat <- function(upfront, running, maturity, R=0.4){
    ## computes the implied hazard rate of the cds based on the upfront
    ## and running quotes, as well as maturity and recovery
    cs <- couponSchedule(nextIMMDate(today()), maturity, "Q", "FIXED", running)
    lambda <- 0.05
    eps <- 1e-8
    while(abs(cdspv(lambda, cs, R) + upfront) > eps){
        lambda <- lambda - (upfront + cdspv(lambda, cs, R))/dcdspv(lambda, cs, R)
    }
    return(lambda)
}

cdshazardrate <- function(quotes, R=0.4){
    ## bootstrap the implied hazard rate curve of the cds based on the upfront
    ## and running quotes, as well as maturity and recovery
    previous.maturity <- today()
    hvec <- c()
    previous.hvec <- c()
    eps <- 1e-8
    previous.cs <- data.frame()
    for(i in 1:nrow(quotes)){
        maturity <- quotes$maturity[i]
        cs <- couponSchedule(nextIMMDate(today()), maturity, "Q", "FIXED", quotes$running[i])
        new.h <- 0.05
        flength <- nrow(cs) - nrow(previous.cs)
        hvec <- c(previous.hvec, rep(new.h, flength))
        sc <- new("defaultprepaycurve", dates=cs$dates, hazardrates=hvec)
        index <- c(rep(0, length(previous.hvec)), rep(1, flength))

        while(abs(cdspv(cs, sc, R) + quotes$upfront[i]) > eps){
            new.h <- new.h - (quotes$upfront[i] + cdspv(cs, sc, R))/dcdspv(cs, sc, R, index)
            hvec <- c(previous.hvec, rep(new.h, flength))
            sc@hazardrates <- hvec
        }
        previous.hvec <- hvec
        previous.maturity <- maturity
        previous.cs <- cs
    }
    return(sc)
}


bonddp <- function(collateral, R=0.7){
    ## bootstrap the implied hazard rate curve of a bond based on the upfront
    ## and running quotes, as well as maturity and recovery
    previous.maturity <- today()
    hvec <- c()
    previous.hvec <- c()
    eps <- 1e-8
    previous.cs <- data.frame()
    for(i in 1:nrow(quotes)){
        maturity <- quotes$maturity[i]
        cs <- couponSchedule(collateral$nextpaydate, collateral$maturity,
                             collateral$frequency, collateral$fixedorfloat,
                             collateral$grosscoupon*0.01, collateral$spread*0.01)
        new.h <- 0.05
        flength <- nrow(cs)-nrow(previous.cs)
        hvec <- c(previous.hvec, rep(new.h, flength))
        sc <- new("survivalcurve", dates=cs$dates, hazardrates=hvec)
        index <- c(rep(0, length(previous.hvec)), rep(1, flength))

        while(abs(cdspv(cs, sc, R) + quotes$upfront[i]) > eps){
            new.h <- new.h - (quotes$upfront[i] + cdspv(cs, sc, R))/dcdspv(cs, sc, R, index)
            hvec <- c(previous.hvec, rep(new.h, flength))
            sc@hazardrates <- hvec
        }
        previous.hvec <- hvec
        previous.maturity <- maturity
        previous.cs <- cs
    }
    return(sc)
}

indexpv <- function(portfolio, index, epsilon=0){
    ## computes the intrinsic index pv of a portfolio of cds
    r <- rep(0, length(portfolio))
    cs <- couponSchedule(nextIMMDate(today()), index$maturity, "Q", "FIXED", index$coupon)
    for(i in 1:length(portfolio)){
        if(epsilon!=0){
            tweakedcurve <- portfolio[[i]]@curve
            tweakedcurve@hazardrates <- tweakedcurve@hazardrates * (1 + epsilon)
            r[i] <- cdspv(cs, tweakedcurve, index$recovery)
        }else{
            r[i] <- cdspv(cs, portfolio[[i]]@curve, index$recovery)
        }
    }
    return( 1+mean(r) )
}

tweakcurves <- function(portfolio, index){
    ## computes the tweaking factor
    epsilon <- 0
    f <- function(epsilon, ...){
        abs(indexpv(portfolio, index, epsilon)-index$indexref)
    }
    epsilon <- optimize(f, c(-0.5, 0.5), portfolio, index, tol=1e-6)$minimum
    portfolio.new <- portfolio
    for(i in 1:length(portfolio)){
        portfolio.new[[i]]@curve@hazardrates <- portfolio[[i]]@curve@hazardrates * (1 + epsilon)
    }
    return( portfolio.new )
}

SPmatrix <- function(portfolio, index){
    cs <- couponSchedule(nextIMMDate(today()), index$maturity, "Q", "FIXED", index$coupon)
    SP <- matrix(0, length(portfolio), length(cs$dates))
    for(i in 1:length(portfolio)){
        SP[i,] <- 1 - PD(portfolio[[i]]@curve)[1:length(cs$dates)]
    }
    return( SP )
}