library("methods")
library("doParallel")
library("itertools")

source("cds_utils.R")

## TODO:
## - switch hazard rates and prepay curves to functions
## - allow to compute forward quantities (works for couponleg and defauleg now,
##   if provide startdate

setClass("abstractcurve")
## flat hazard rate curve
setClass("flatcurve", contains="abstractcurve", representation(h="numeric"))
## shaped curve: the hazard rate curve is given by scaling a base shape by a factor
## k(h)=alpha*exp(-beta*h)
setClass("shapedcurve", contains="abstractcurve", representation(h="numeric",
                        shape="function", alpha="numeric", beta="numeric"))
setClass("defaultcurve", contains="abstractcurve", representation(dates="Date", hazardrates="numeric"))
setClass("defaultprepaycurve", representation(prepayrates="numeric"), contains="defaultcurve")
setClass("creditcurve", representation(issuer="character", startdate="Date",
                                       recovery="numeric", curve="defaultcurve"))

shapedtodpc <- function(cs, sc, startdate=Sys.Date()){
    ## convert a shaped curve to a standard defaultprepaycuve
    T <- yearFrac(startdate, cs$dates)
    hvec <- sc@shape(T) * sc@h
    kvec <- sc@alpha * exp(-sc@beta * hvec)
    dpc <- new("defaultprepaycurve", hazardrates=hvec, prepayrates=kvec, dates=cs$dates)
    return (dpc)
}

## define couponleg generic
setGeneric("couponleg", function(cs, sc, ...) {
    standardGeneric("couponleg")
})

## write couponleg methods for the four types of curves
setMethod("couponleg", signature("data.frame", "flatcurve"),
          function(cs, sc, startdate=Sys.Date() + 1, accruedondefault=TRUE){
              stopifnot(class(startdate)=="Date")
              stopifnot(is.logical(accruedondefault))
              T <- yearFrac(startdate, cs$dates)
              PD <- exp(- sc@h * T )
              if(accruedondefault){
                  PDadj <- 0.5 * (c(1, PD[-length(PD)]) + PD)
              }else{
                  PDadj <- PD
              }
              return( as.numeric(crossprod(PDadj, cs$coupons * cs$df)) )
          })

setMethod("couponleg", signature("data.frame", "defaultcurve"),
          ## 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.
          function(cs, sc, startdate=Sys.Date() + 1, accruedondefault=TRUE){
              stopifnot(class(startdate)=="Date")
              stopifnot(is.logical(accruedondefault))
              x1T <- yearFrac(startdate, sc@dates)
              x2T <- yearFrac(startdate, cs$dates)
              dT <- diff(c(0, yearFrac(startdate, cs$dates)))
              hfun <- approxfun(x1T, sc@hazardrates, method="constant", rule=2)
              PD <- cumprod(exp(- hfun(x2T) * dT))
              if(accruedondefault){
                  PDadj <- 0.5 * (c(1, PD[-length(PD)]) + PD)
              }else{
                  PDadj <- PD
              }
              return( as.numeric(crossprod(PDadj, cs$coupons * cs$df)) )
          })


setMethod("couponleg", signature("data.frame", "defaultprepaycurve"),
          ## computes the pv of the risky coupon leg from the coupon schedule,
          ## a hazard rate curve, and a prepay curve. We assume the poisson processes driving
          ## default and prepayment are independent, so the intensity of either event
          ## happenning is the sum of the two.
          function(cs, sc, startdate=Sys.Date() + 1, accruedondefault=TRUE){
              stopifnot(class(startdate)=="Date")
              stopifnot(is.logical(accruedondefault))
              x1T <- yearFrac(startdate, sc@dates)
              x2T <- yearFrac(startdate, cs$dates)
              dT <- diff(c(0, yearFrac(startdate, cs$dates)))
              hfun <- approxfun(x1T, sc@hazardrates, method="constant", rule=2)
              pfun <- approxfun(x1T, sc@prepayrates, method="constant", rule=2)
              SP <- cumprod(exp( - (hfun(x2T) + pfun(x2T)) * dT))
              if(accruedondefault){
                  SPadj <- 0.5 * (c(1, SP[-length(SP)]) + SP)
              }else{
                  SPadj <- SP
              }
              return( as.numeric(crossprod(SPadj, 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) )
}

setMethod("couponleg", signature("data.frame", "shapedcurve"),
          ## computes the pv of the risky coupon leg from the coupon schedule,
          ## a hazard rate curve, and a prepay curve. We assume the poisson processes driving
          ## default and prepayment are independent, so the intensity of either event
          ## happenning is the sum of the two.
          function(cs, sc, startdate=Sys.Date() + 1, accruedondefault=TRUE){
              stopifnot(class(startdate)=="Date")
              stopifnot(is.logical(accruedondefault))
              return( couponleg(cs, shapedtodpc(cs, sc, startdate), startdate, accruedondefault) )
          })

## define dcouponleg generic
setGeneric("dcouponleg", function(cs, sc, startdate, ...) {
    standardGeneric("dcouponleg")
})

setMethod("dcouponleg", signature("data.frame", "flatcurve", "Date"),
          function(cs, sc, startdate=Sys.Date() + 1, accruedondefault=TRUE){
              stopifnot(is.logical(accruedondefault))
              T <- yearFrac(startdate, cs$dates)
              dPD <- -T * exp(- sc@h * T )
              if(accruedondefault){
                  dPDadj <- 0.5 * (c(0, dPD[-length(dPD)]) + dPD)
              }else{
                  dPDadj <- dPD
              }
              return( as.numeric(crossprod(dPDadj, cs$coupons * cs$df)) )
          })

setMethod("dcouponleg", signature("data.frame", "defaultcurve", "Date"),
          ## derivative of couponleg with respect to hazardrate
          function(cs, sc, startdate=Sys.Date() + 1, index, accruedondefault=TRUE) {
              stopifnot(is.logical(accruedondefault))
              dT <- diff(c(0, yearFrac(startdate, 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( as.numeric(crossprod( dPDadj, cs$coupons * cs$df)) )
          })

setMethod("dcouponleg", signature("data.frame", "shapedcurve", "Date"),
          function(cs, sc, startdate=Sys.Date() + 1, accruedondefault = TRUE){
              T <- yearFrac(startdate, cs$dates)
              dT <- diff(c(0, T))
              hvec <- sc@h * sc@shape(T)
              kvec <- sc@alpha*exp(-sc@beta*hvec)
              SP <- cumprod(exp( - (hvec + kvec ) * dT))
              dSP <- -cumsum( dT * sc@shape(T) * (1 - sc@beta * kvec)) * SP
              if(accruedondefault) {
                  dSPadj <- 0.5 *(c(0, dSP[-length(SP)]) + dSP)
              }else{
                  dSPadj <- dSP
              }
              return( as.numeric(crossprod(dSPadj, cs$coupons * cs$df)) )
          })

## setMethod("dcouponleg", signature("data.frame", "defaultprepaycurve", "numeric"),
##           ## 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)
##           function(cs, dpc, index, beta, accruedondefault=TRUE) {
##           dT <- diff(c(0, yearFrac(Sys.Date(), 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( as.numeric(crossprod(dSPadj, cs$coupons * cs$df)) )
##       })

## define cdsduration generic
setGeneric("cdsduration", function(sc, maturity, ...) {
    standardGeneric("cdsduration")
})

## duration is based on the standard IMM schedule
setMethod("cdsduration", signature("abstractcurve", "Date"),
          ## computes the risky PV01, also called risky annuity of a cds
          function(sc, maturity, tradedate=Sys.Date(), accruedondefault=TRUE){
              stopifnot(is.logical(accruedondefault))
              cs <- couponSchedule(IMMDate(tradedate, noadj=TRUE), maturity,"Q", "FIXED", 1,
                                   0, tradedate, IMMDate(tradedate, "prev"))
              startdate <- tradedate+1
              acc <- cdsAccrued(tradedate, 1)
              return( couponleg(cs, sc, startdate, accruedondefault=accruedondefault) - acc)
          })

## define cdsduration generic
setGeneric("cdsspread", function(sc, recovery,  maturity, ...) {
    standardGeneric("cdsspread")
})

## duration is based on the standard IMM schedule
setMethod("cdsspread", signature("abstractcurve", "numeric", "Date"),
          ## computes the risky PV01, also called risky annuity of a cds
          function(sc, recovery, maturity, tradedate=Sys.Date(), accruedondefault=TRUE){
              stopifnot(is.logical(accruedondefault))
              cs <- couponSchedule(IMMDate(tradedate, noadj=TRUE), maturity,"Q", "FIXED", 1,
                                   0, tradedate, IMMDate(tradedate, "prev"))
              startdate <- tradedate+1
              acc <- cdsAccrued(tradedate, 1)

              return( defaultleg(cs, sc, recovery, startdate)/
                     (couponleg(cs, sc, startdate, accruedondefault=accruedondefault) - acc))
          })

## define defaultleg generic
setGeneric("defaultleg", function(cs, sc, recovery, ...) {
    standardGeneric("defaultleg")
})

## write defaultleg methods for the four types of curves
setMethod("defaultleg", signature("data.frame", "flatcurve", "numeric"),
          ## Computes the pv of the default leg of a cds based on a given
          ## coupon schedule, flat hazard rate, and recovery.
          function(cs, sc, recovery, startdate=Sys.Date()+1){
              T <- yearFrac(startdate, cs$dates)
              dT <- diff(c(0, T))
              Q <- exp(-sc@h * T) * cs$df
              Qmid <- 1/2 * (c(1, Q[-length(Q)]) + Q)
              r <- (1 - recovery) * crossprod(sc@h * Qmid, dT)
              return( as.numeric(r) )
          })

setMethod("defaultleg", signature("data.frame", "defaultcurve", "numeric"),
          ## Computes the pv of the default leg of a cds based on a given
          ## coupon schedule, hazard rate curve, and recovery.
          function(cs, sc, recovery, startdate=Sys.Date()+1){
              T <- yearFrac(startdate, cs$dates)
              x1T <- yearFrac(startdate, sc@dates)
              hfun <- approxfun(x1T, sc@hazardrates, method="constant", rule=2)
              dT <- diff(c(0, T))
              Q <- cumprod(exp(-hfun(T) * dT)) * cs$df
              Qmid <- 1/2 * (c(1, Q[-length(Q)]) + Q)
              r <- (1 - recovery) * crossprod(hfun(T) * Qmid, dT)
              return( as.numeric(r) )
          })

setMethod("defaultleg", signature("data.frame", "defaultprepaycurve", "numeric"),
          ## Computes the pv of the default leg of a cds based on a given
          ## coupon schedule, defaultprepaycurve, and recovery.
          function(cs, sc, recovery, startdate=Sys.Date()+1){
              stopifnot(class(startdate)=="Date")
              x2T <- yearFrac(startdate, cs$dates)
              x1T <- yearFrac(startdate, sc@dates)
              dT <- diff(c(0, x2T))
              hfun <- approxfun(x1T, sc@hazardrates, method = "constant", rule=2)
              pfun <- approxfun(x1T, sc@prepayrates, method = "constant", rule=2)
              Q <- cumprod(exp(- (hfun(x2T)+pfun(x2T)) * dT)) * cs$df
              Qmid <- 1/2 * (c(1, Q[-length(Q)]) + Q)
              r <- (1 - recovery) * crossprod(hfun(x2T) * Qmid, dT)
              return( as.numeric(r) )
})

setMethod("defaultleg", signature("data.frame", "shapedcurve", "numeric"),
          ## Computes the pv of the default leg of a cds based on a shaped curve.
          function(cs, sc, recovery, startdate=Sys.Date()+1){
              return( defaultleg(cs, shapedtodpc(cs, sc), recovery) )
          })

## define ddefaultleg generic
setGeneric("ddefaultleg", function(cs, sc, recovery, startdate, ...) {
    standardGeneric("ddefaultleg")
})

setMethod("ddefaultleg", signature("data.frame", "flatcurve", "numeric", "Date"),
          ## derivative of defaultleg with respect to flat hazardrate
          function(cs, sc, recovery, startdate=Sys.Date()+1){
              T <- yearFrac(startdate, cs$dates)
              dT <- diff(c(0, T))
              dQ <- - T * exp(-sc@h * T) * cs$df
              Q <- exp(-sc@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) +sc@h * crossprod(dQmid, dT))
              return( as.numeric(dr) )
          })

setMethod("ddefaultleg", signature("data.frame", "defaultcurve", "numeric", "Date"),
          ## derivative of defaultleg with respect to hazardrate
          function(cs, sc, recovery, startdate=Sys.Date()+1, index){
              T <- yearFrac(startdate, 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( as.numeric(dr) )
})

setMethod("ddefaultleg", signature("data.frame", "shapedcurve", "numeric", "Date"),
          function(cs, sc, recovery, startdate=Sys.Date()+1) {
              T <- yearFrac(startdate, cs$dates)
              dT <- diff(c(0, T))
              hvec <- sc@shape(T) * sc@h
              kvec <- sc@alpha * exp(-sc@beta * hvec)
              Q <- cumprod(exp( -(hvec + kvec) * dT)) * cs$df
              dQ <- -cumsum(sc@shape(T) * dT * (1 - sc@beta * kvec)) * Q
              Qmid <- 1/2 * (c(1, Q[-length(Q)]) + Q)
              dQmid <- 1/2 * (c(0, dQ[-length(dQ)]) + dQ)
              dr <- (1-recovery)* (crossprod(sc@shape(T) * Qmid, dT) + crossprod(hvec * dQmid, dT))
              return( as.numeric(dr) )
          })

## test.shape <- splinefun(0:5, seq(0.5,1,length=6))
## eps <- 1e-8
## test.sc.flat <- new("flatcurve", h=0.07)
## test.sc <- new("shapedcurve", h=0.07, shape=test.shape, alpha=.25, beta=15)
## test.scplus <- new("shapedcurve", h=0.07+eps, shape=test.shape, alpha=.25, beta=15)
## test.scminus <- new("shapedcurve", h=0.07-eps, shape=test.shape, alpha=.25, beta=15)
## (contingentleg(cs, test.scplus, 0.4) - contingentleg(cs, test.scminus, 0.4))/(2*eps)
## dcontingentleg(cs, test.sc)

## define contingentleg generic
setGeneric("contingentleg", function(cs, sc, recovery, ...) {
    standardGeneric("contingentleg")
})

## write contingentleg methods for the four types of curves
setMethod("contingentleg", signature("data.frame", "flatcurve", "numeric"),
          ## Computes the pv of the contingent leg of a cds based on a given
          ## coupon schedule, flat hazard rate, and recovery.
          function(cs, sc, recovery, startdate=Sys.Date()){
              stopifnot(class(startdate)=="Date")
              T <- yearFrac(startdate, cs$dates)
              dT <- diff(c(0, T))
              Q <- exp(-sc@h * T) * cs$df
              Qmid <- 1/2 * (c(1, Q[-length(Q)]) + Q)
              r <- Q[length(cs$dates)] + recovery  * sc@h * crossprod(Qmid, dT)
              return( as.numeric(r))
          })

setMethod("contingentleg", signature("data.frame", "defaultcurve", "numeric"),
          ## Computes the pv of the contingent leg of a cds based on a given
          ## coupon schedule, flat hazard rate, and recovery.
          function(cs, sc, recovery, startdate=Sys.Date()){
              T <- yearFrac(startdate, 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 <- Q[length(cs$dates)] + recovery * crossprod(sc@hazardrates[1:length(dT)] * Qmid, dT)
              return( as.numeric(r))
          })

setMethod("contingentleg", signature("data.frame", "defaultprepaycurve", "numeric"),
          ## Computes the pv of the contingent leg of a cds based on a given
          ## coupon schedule, hazard rates curve, prepay curve, and recovery.
          function(cs, sc, recovery, startdate=Sys.Date()) {
              stopifnot(class(startdate)=="Date")
              x1T <- yearFrac(startdate, sc@dates)
              x2T <- yearFrac(startdate, cs$dates)
              dT <- diff(c(0, yearFrac(startdate, cs$dates)))
              hfun <- approxfun(x1T, sc@hazardrates, method="constant", rule=2)
              pfun <- approxfun(x1T, sc@prepayrates, method="constant", rule=2)
              Q <- cumprod(exp( -(hfun(x2T)+pfun(x2T)) * dT)) * cs$df
              Qmid <- 1/2 * (c(1, Q[-length(Q)]) + Q)
              r <- recovery * crossprod(hfun(x2T) * Qmid, dT) +
                  crossprod(pfun(x2T) * Qmid, dT) + Q[length(cs$dates)]
              return( as.numeric(r) )
          })

setMethod("contingentleg", signature("data.frame", "shapedcurve", "numeric"),
          function(cs, sc, recovery, startdate=Sys.Date()+1){
              stopifnot(class(startdate)=="Date")
             return( contingentleg(cs, shapedtodpc(cs, sc, startdate), recovery, startdate) )
         })

## define dcontingentleg generic
setGeneric("dcontingentleg", function(cs, sc, recovery, startdate, index, ...) {
    standardGeneric("dcontingentleg")
})

setMethod("dcontingentleg", signature("data.frame", "defaultcurve", "numeric", "Date", "numeric"),
          ## derivative of contingentleg with respect to hazardrate
          function(cs, sc, recovery, startdate=Sys.Date(), index){
              T <- yearFrac(startdate, 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 <- dQ[length(cs$dates)] + recovery * crossprod(index * Qmid, dT) +
                  recovery * crossprod(sc@hazardrates[1:length(dT)] * dQmid, dT)
              return( as.numeric(dr) )
          })

setMethod("dcontingentleg", signature("data.frame", "defaultcurve", "numeric", "Date", "missing"),
          ## derivative of contingentleg with respect to hazardrate
          function(cs, sc, recovery, startdate=Sys.Date()){
              ## derivative of contingentleg with respect to hazardrate
              T <- yearFrac(startdate, cs$dates)
              dT <- diff(c(0, T))
              Q <- exp(-sc@h * T) * cs$df
              dQ <- -T * exp(- sc@h * T) * cs$df
              Qmid <-  1/2 * (c(1, Q[-length(Q)]) + Q)
              dr <- dQ[length(dQ)] + recovery * crossprod(Qmid, dT) +
                  recovery * h * crossprod(1/2 * (c(0, dQ[-length(dQ)]) + dQ), dT)
              return( as.numeric(dr) )
          })

setMethod("dcontingentleg", signature("data.frame", "shapedcurve", "numeric", "Date", "missing"),
          ## Computes the pv of the contingent leg of a cds based on a given
          ## coupon schedule, hazard rates curve, prepay curve, and recovery.
          function(cs, sc, recovery, startdate=Sys.Date()){
              T <- yearFrac(startdate, cs$dates)
              dT <- diff(c(0, T))
              hvec <- sc@shape(T) * sc@h
              kvec <- sc@alpha * exp( - sc@beta *hvec)
              Q <- cumprod(exp( -(hvec + kvec) * dT)) * cs$df
              dQ <- -cumsum(sc@shape(T) * dT * (1 - sc@beta * kvec)) * Q
              Qmid <- 1/2 * (c(1, Q[-length(Q)]) + Q)
              dQmid <- 1/2 * (c(0, dQ[-length(dQ)]) + dQ)
              dr <- recovery * (crossprod(sc@shape(T) * Qmid, dT) + crossprod(hvec * dQmid, dT)) +
                  crossprod(-sc@beta * sc@shape(T) * kvec * Qmid, dT) +
                      crossprod(kvec * dQmid, dT) + dQ[length(cs$dates)]
              return( as.numeric(dr) )
          })

cdspv <- function(cs, sc, recovery, tradedate){
    startdate <- tradedate + 1
    return ( couponleg(cs, sc, startdate) - defaultleg(cs, sc, recovery, startdate))
}

cds.parspread <- function(cs, sc, recovery, tradedate){
    ## computes exact cds running spread for a cds
    ## should be very close to hazardrate * (1-recovery)
    startdate <- tradedate + 1
    return (defaultleg(cs, sc, recovery, startdate)/
            (couponleg(cs, sc, startdate)-cdsAccrued(tradedate, 1)))
}

dcds.parspread <- function(cs, sc, recovery, tradedate){
    startdate <- tradedate + 1
    u <- defaultleg(cs, sc, recovery, startdate)
    v <- couponleg(cs, sc, startdate)-cdsAccrued(tradedate, 1)
    u_prime <- ddefaultleg(cs, sc, recovery, startdate)
    v_prime <- dcouponleg(cs, sc, startdate)
    return( (u_prime*v-u*v_prime)/v^2 )
}

snacpv <- function(cs, parspread, coupon, recovery, tradedate){
    ## compute the equivalent pv of a snac quote
    startdate <- tradedate + 1
    sc <- new("flatcurve", h=parspread/(1-recovery))
    epsilon <- 1e-8
    while(abs(cds.parspread(cs, sc, recovery, tradedate)-parspread)>epsilon){
        sc@h <- sc@h - (cds.parspread(cs, sc, recovery, tradedate)-parspread)/
            dcds.parspread(cs, sc, recovery, tradedate)
    }
    acc <- cdsAccrued(tradedate, coupon)
    return( 1 + couponleg(cs, sc, startdate)*coupon-defaultleg(cs, sc, recovery, startdate)-acc )
}

dcdspv <- function(cs, sc, recovery, tradedate, index){
    startdate <- tradedate + 1
    if(missing(index)){
        return(dcouponleg(cs, sc, startdate) - ddefaultleg(cs, sc, recovery, startdate))
    }else{
        return ( dcouponleg(cs, sc, startdate, index) - ddefaultleg(cs, sc, recovery, startdate, index) )
    }
}

bondpv <- function(cs, sc, recovery, tradedate=Sys.Date()){
    return( contingentleg(cs, sc, recovery, tradedate) + couponleg(cs, sc, tradedate) )
}

dbondpv <- function(cs, sc, recovery, startdate, index){
    if(missing(index)){
        return( dcontingentleg(cs, sc, recovery, startdate) + dcouponleg(cs, sc, startdate))
    }else{
        return( dcontingentleg(cs, sc, recovery, startdate, index) + dcouponleg(cs, sc, startdate, index) )
    }
}

cdshazardrate.flat <- function(upfront, running, maturity, R=0.4, tradedate){
    ## computes the implied hazard rate of the cds based on the upfront
    ## and running quotes, as well as maturity and recovery
    cs <- couponSchedule(IMMDate(tradedate, noadj=TRUE), maturity,"Q", "FIXED", running,
                           0, tradedate, IMMDate(tradedate, "prev"))
    sc <- new("flatcurve", h = 0.05)
    eps <- 1e-8
    acc <- cdsAccrued(tradedate, running)
    upfront <- upfront - acc
    while(abs(cdspv(cs, sc, R, tradedate) + upfront) > eps){
        sc@h <- sc@h - (upfront + cdspv(cs, sc, R, tradedate))/dcdspv(cs, sc, R, tradedate)
    }
    return(sc@h)
}

cdshazardrate.shaped <- function(upfront, running, maturity, shape, R=0.4) {
    cs <- couponSchedule(IMMDate(Sys.Date()), maturity, "Q", "FIXED", running)
    sc <- new("shapedcurve", shape=shape, h=0.05, alpha=0.25, beta=15)
    eps <- 1e-8
    while(abs(cdspv(cs, sc, R) + upfront) > eps){
        sc@h <- sc@h - (upfront + cdspv(cs, sc, R))/dcdspv(cs, sc, R)
    }
    return(sc)
}

cdshazardrate <- function(quotes, R=0.4, tradedate=Sys.Date(), cs.all){
    ## bootstrap the implied hazard rate curve of the cds based on the upfront
    ## and running quotes, as well as maturity and recovery
    previous.maturity <- tradedate
    hvec <- c()
    previous.hvec <- c()
    eps <- 1e-8
    previous.cs <- data.frame()
    sc <- NULL
    for(i in 1:nrow(quotes)){
        if(is.na(quotes$upfront[i])){
            next
        }
        maturity <- quotes$maturity[i]
        if(tradedate+1>=maturity){
            next
        }
        cs <- cs.all[cs.all$unadj.dates<=maturity,]
        cs$coupons <- cs$coupons * quotes$running[i]
        new.h <- 0.05
        flength <- nrow(cs) - nrow(previous.cs)
        hvec <- c(previous.hvec, rep(new.h, flength))
        sc <- new("defaultcurve", dates=cs$dates, hazardrates=hvec)
        index <- c(rep(0, length(previous.hvec)), rep(1, flength))
        acc <- cdsAccrued(tradedate, quotes$running[i])
        while(abs(cdspv(cs, sc, R, tradedate) + quotes$upfront[i] - acc) > eps){
            new.h <- new.h - (quotes$upfront[i] + cdspv(cs, sc, R, tradedate) - acc)/
                dcdspv(cs, sc, R, tradedate, index)
            hvec <- c(previous.hvec, rep(new.h, flength))
            sc@hazardrates <- hvec
        }
        previous.hvec <- hvec
        previous.maturity <- maturity
        previous.cs <- cs
    }
    return(sc)
}

bondhazardrate.shaped <- function(collateral, shape, R=0.4, alpha=0.25, beta=15, startdate=Sys.Date()){
    ## calibrate a default prepay curve to the collateral information
    cs <- couponSchedule(collateral$nextpaydate, collateral$maturity,
                         collateral$frequency, collateral$fixedorfloat,
                         collateral$grosscoupon*0.01, collateral$spread*0.01, startdate)
    sc <- new("shapedcurve", h=0.05, shape=shape, alpha=alpha, beta=beta)
    eps <- 1e-8
    counter <- 0
    if(collateral$price<2){
        cat("price is too low\n")
        sc <- new("shapedcurve", h=1e6, shape=shape, alpha=alpha, beta=beta)
        return( shapedtodpc(cs, sc, startdate) )
    }

    while(abs(bondpv(cs, sc, R, startdate) - collateral$price/100) > eps){
        dh <- (collateral$price/100 - bondpv(cs, sc, R, startdate))/dbondpv(cs, sc, R, startdate)
        while(sc@h+dh<0){
            dh <- 0.5 * dh
        }
        sc@h <- sc@h+dh
        counter <- counter+1
        if(counter >= 100){
            return( NULL )
            stop("didn't reach convergence")
        }
    }
    return( shapedtodpc(cs, sc, startdate) )
}

tweakportfolio <- function(portfolio, epsilon, multiplicative=TRUE, forward.tweak=1){
    ## tweak a portfolio of creditcurves
    ## if multiplicative is TRUE apply a multiplicative tweak
    ## otherwise apply an additive one
    if(multiplicative){
        r <- lapply(portfolio, function(x) {
            x@curve@hazardrates[forward.tweak:length(x@curve@hazardrates)] <-
                x@curve@hazardrates[forward.tweak:length(x@curve@hazardrates)] * (1+epsilon)
            x
        })
    }else{
        ## we do a tweak to the spread
        r <- lapply(portfolio, function(x) {
            x@curve@hazardrates <- x@curve@hazardrates + epsilon/(1-x@recovery)
            x
        })
    }
    return( r )
}

indexpv <- function(index, epsilon=0, tradedate=Sys.Date(), clean=TRUE,
                    maturity=index$maturity, forward.tweak=1, check=FALSE){
    ## computes the intrinsic price of a portfolio of cds
    ## If maturity is provided, only computes the pv up to that point
    ## (Say we compute the 3 year pv based on 5 year curves
    ## forward.tweak only makes sense if epsilon is non zero
    ## and will teak the curves starting from forward.index
    if(epsilon!=0){
        portfolio <- tweakportfolio(index$portfolio, epsilon, forward.tweak=forward.tweak)
    }else{
        portfolio <- index$portfolio
    }
    startdate <- tradedate + 1
    #ugly hack
    cs <- index$cs[index$cs$dates-maturity< 3,]
    cl.list <- unlist(lapply(portfolio, function(x){couponleg(cs, x@curve, startdate)}))
    pl.list <- unlist(lapply(portfolio, function(x){defaultleg(cs, x@curve, x@recovery, startdate)}))
    if(check){
        if(any(is.na(cl.list))){
            logerror(paste("couldn't compute single name protection leg for",
                           index$portfolio[[which(is.na(pl.list))]]@issuer))
            return(NA)
        }
        if(any(is.na(pl.list))){
            logerror(paste("couldn't compute single name protection leg for",
                           index$portfolio[[which(is.na(pl.list))]]@issuer))
            return(NA)
        }
    }
    spread <- index$quotes$spread[index$quotes$maturity==maturity]
    r <- list(cl = spread * mean(cl.list), pl = mean(pl.list), bp = 1+mean(spread*cl.list-pl.list))
    if(clean){
        accrued <- cdsAccrued(tradedate, spread)
        r$bp <- r$bp-accrued
        r$cl <- r$cl-accrued
    }
    return(r)
}

indexduration <- function(index, tradedate, maturity=index$maturity){
    ## compute the duration of a portfolio of survival curves
    cs <- index$cs[index$cs$dates-maturity< 3,]
    startdate <- tradedate+1
    cl.list <- unlist(lapply(index$portfolio, function(x){couponleg(cs, x@curve, startdate)}))
    return( mean(cl.list) )
}

indexspread <- function(index, tradedate=Sys.Date()){
    ## computes the spread of a portfolio of survival curves
    ## S <- 0
    ## d <- rep(0, length(portfolio))
    ## for(i in 1:length(portfolio)){
    ##     d[i] <- cdsduration(portfolio[[i]]@curve, index$maturity)
    ##     S <- S + d[i] * cdsspread(portfolio[[i]]@curve, index$maturity, portfolio[[i]]@recovery)
    ## }
    temp <- indexpv(index, tradedate=tradedate)
    S <- index$quotes$spread*(1-(temp$bp-1)/temp$cl)
    return(S)
}

indextheta <- function(index, tradedate=Sys.Date()){
    current.pv <- indexpv(index, tradedate=tradedate)$bp
    newmaturity <- index$cs$unadj.dates[nrow(index$cs)-4]
    index$quotes$maturity <- newmaturity
    forward.pv <- indexpv(index, tradedate=tradedate, maturity=newmaturity)$bp
    theta <- forward.pv-current.pv+index$quotes$spread
    return( theta )
}


portfoliospread <- function(portfolio, maturity, tradedate=Sys.Date()){
    ## computes the spread of a portfolio defined by notionals and survivalcurves
    ## for a given maturity.
    ## if maturity is missing, we use the intrinsic maturity for each curve
    if(missing(maturity)){
        maturityvec <- as.Date(sapply(portfolio$SC, creditcurve.maturity), origin="1970-01-01")
    }else{
        maturityvec <- rep(maturity, length(portfolio$SC))
    }

    data <- foreach(it = izip(creditcurve = portfolio$SC, notional = portfolio$notional,
                        maturity = maturityvec, count=icount(length(portfolio$notional))), .combine=rbind) %:% {
            when(length(it$creditcurve@curve@hazardrates)!=0 && it$maturity > it$creditcurve@startdate
                 && it$creditcurve@curve@hazardrates[1]<=10)
        } %dopar% {
            d <- cdsduration(it$creditcurve@curve, it$maturity, tradedate)
            S <- cdsspread(it$creditcurve@curve, it$creditcurve@recovery, it$maturity, tradedate)
            c(d * S * it$notional, d * it$notional)
        }
    if(is.null(dim(data))){
        return(data[1]/data[2])
    }else{
        return(sum(data[,1])/sum(data[,2]))
    }
}

portfolioduration <- function(portfolio, maturity){
    durations <- unlist(lapply(portfolio$SC, function(x){cdsduration(x@curve, maturity)}))
    return( crossprod(durations, portfolio$notional)/sum(portfolio$notional) )
}

tweakcurves <- function(index, tradedate=Sys.Date()){
    ## computes the tweaking factor
    epsilon <- rep(0, nrow(index$quotes))
    for(i in 1:nrow(index$quotes)){
        if(i==1){
            forward.tweak <- 1
        }else{
            forward.tweak <- which.min(index$cs$dates>index$quotes$maturity[i-1])
        }
        f <- function(epsilon, ...){
            abs(indexpv(index, epsilon, tradedate, maturity=index$quotes$maturity[i],
                        forward.tweak=forward.tweak)$bp - index$quotes$price[i])
        }
        epsilon[i] <- optimize(f, c(-0.15, 0.30), index, tol=1e-6)$minimum
        index$portfolio <- tweakportfolio(index$portfolio, epsilon[i], forward.tweak=forward.tweak)
        loginfo(paste("tweak =", epsilon[i]))
    }
    return( list(portfolio=index$portfolio, basis=epsilon) )
}

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

SP <- function(sc, startdate=Sys.Date()){
    ## computes the survival probability associated with the survival curve
    T <- c(0, yearFrac(startdate, sc@dates))
    dT <- diff(T)
    return( cumprod(exp(-sc@hazardrates * dT)) )
}

SPmatrix <- function(portfolio, n.dates){
    ## computes matrix of survival probability
    ## inputs:
    ##   portfolio: portfolio of survival curves
    ##   index: index representation
    ## ouput:
    ##   matrix of survival probabilities of dimensions dim1 x dim2
    ##   with dim1 number of issuers and dim2 number of dates in the
    ##   coupon schedule of index
    SP <- matrix(0, length(portfolio), n.dates)
    for(i in 1:length(portfolio)){
        SP[i,] <- SP(portfolio[[i]]@curve, portfolio[[i]]@startdate)[1:n.dates]
    }
    return( SP )
}

DP2 <- function(sc, dates, startdate=Sys.Date()){
    ## computes the default probability and prepay probability associated
    ## with the survival curve at the dates specified by dates
    if(length(sc@hazardrates)==0){
        return(list(defaultprob=rep(0,length(dates)),
                    prepayprob=rep(0, length(dates))))
    }
    x2T <- yearFrac(startdate, dates)
    dT <- diff(c(0, x2T))
    x1T <- yearFrac(startdate, sc@dates)
    hfun <- approxfun(x1T, sc@hazardrates, method="constant", rule=2)
    pfun <- approxfun(x1T, sc@prepayrates, method="constant", rule=2)
    Qmid <- exp(-cumsum((hfun(x2T)+pfun(x2T)) * dT))
    return(list(defaultprob = cumsum(hfun(x2T) * Qmid * dT),
                prepayprob = cumsum(pfun(x2T) * Qmid * dT)))
}

SPmatrix2 <- function(portfolio, dealdata, freq="3 months", startdate=Sys.Date()){
    ## computes the default and prepay probability matrix of a portfolio
    ## at the dates specified from dealdata
    dates <- getdealschedule(dealdata, freq)
    dates <- dates[dates>=addBusDay(startdate, 3)]
    DP <- matrix(0, length(portfolio), length(dates))
    PP <- matrix(0, length(portfolio), length(dates))
    for(i in seq_along(portfolio)){
        temp <- DP2(portfolio[[i]]@curve, dates)
        DP[i,] <- temp$defaultprob
        PP[i,] <- temp$prepayprob
        colnames(DP) <- as.character(dates)
        colnames(PP) <- as.character(dates)
    }
    return(list(DP=DP, PP=PP))
}

creditcurve.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")
    }
}