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=today()){
    ## 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, accruedondefault=TRUE){
               T <- yearFrac(today(), 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, accruedondefault=TRUE){
              x1T <- yearFrac(today(), sc@dates)
              x2T <- yearFrac(today(), cs$dates)
              dT <- diff(c(0, x2T))
              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, accruedondefault=TRUE, startdate=today()){
              x1T <- yearFrac(today(), sc@dates)
              x2T <- yearFrac(startdate, cs$dates)
              dT <- diff(c(0, x2T))
              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, accruedondefault=TRUE, startdate=today()){
              return( couponleg(cs, shapedtodpc(cs, sc, startdate), accruedondefault, startdate) )
          })

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

setMethod("dcouponleg", signature("data.frame", "flatcurve", "missing"),
          function(cs, sc, accruedondefault=TRUE){
              T <- yearFrac(today(), 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", "numeric"),
          ## derivative of couponleg with respect to hazardrate
          function(cs, sc, index, accruedondefault=TRUE) {
              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( as.numeric(crossprod( dPDadj, cs$coupons * cs$df)) )
          })

setMethod("dcouponleg", signature("data.frame", "shapedcurve", "missing"),
          function(cs, sc, accruedondefault = TRUE){
              T <- yearFrac(today(), 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{
                  SPadj <- 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(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( as.numeric(crossprod(dSPadj, cs$coupons * cs$df)) )
##       })

## test.shape <- splinefun(0:5, rep(1,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)

## (couponleg(cs, test.scplus)-couponleg(cs, test.scminus))/(2 * eps)


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

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

## 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){
              T <- yearFrac(today(), 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){
              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( 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, hazard rates curve, prepay curves, and recovery.
          function(cs, sc, recovery, startdate=today()){
              x2T <- yearFrac(startdate, cs$dates)
              x1T <- yearFrac(today(), 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){
              return( defaultleg(cs, shapedtodpc(cs, sc), recovery) )
          })

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

setMethod("ddefaultleg", signature("data.frame", "flatcurve", "numeric", "missing"),
          ## derivative of defaultleg with respect to flat hazardrate
          function(cs, sc, recovery){
              T <- yearFrac(today(), 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", "numeric"),
          ## derivative of defaultleg with respect to hazardrate
          function(cs, sc, recovery, index){
              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( as.numeric(dr) )
})

setMethod("ddefaultleg", signature("data.frame", "shapedcurve", "numeric", "missing"),
          function(cs, sc, recovery) {
              T <- yearFrac(today(), 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){
              T <- yearFrac(today(), 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){
              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 <- 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=today()) {
              x1T <- yearFrac(today(), sc@dates)
              x2T <- yearFrac(startdate, cs$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 <- 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=today()){
             return( contingentleg(cs, shapedtodpc(cs, sc), recovery, startdate) )
         })

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

setMethod("dcontingentleg", signature("data.frame", "defaultcurve", "numeric", "numeric"),
          ## derivative of contingentleg with respect to hazardrate
          function(cs, sc, recovery, index){
              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 <- 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", "missing"),
          ## derivative of contingentleg with respect to hazardrate
          function(cs, sc, recovery){
              ## derivative of contingentleg with respect to hazardrate
              T <- yearFrac(today(), 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", "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){
              T <- yearFrac(today(), 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){
    return ( couponleg(cs, sc) - defaultleg(cs, sc, recovery))
}

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=NULL){
    if(is.null(index)){
        return(dcouponleg(cs, sc)-ddefaultleg(cs, sc, recovery))
    }else{
        return ( dcouponleg(cs, sc, index) - ddefaultleg(cs, sc, recovery, index) )
    }
}

bondpv <- function(cs, sc, recovery){
    return( contingentleg(cs, sc, recovery)+couponleg(cs, sc) )
}

dbondpv <- function(cs, sc, recovery, index=NULL){
    if(is.null(index)){
        return( dcontingentleg(cs, sc, recovery) + dcouponleg(cs, sc))
    }else{
        return( dcontingentleg(cs, sc, recovery, index)+dcouponleg(cs, sc, 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)
    sc <- new("flatcurve", h = 0.05)
    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@h)
}

cdshazardrate.shaped <- function(upfront, running, maturity, shape, R=0.4) {
    cs <- couponSchedule(nextIMMDate(today()), 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){
    ## 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)){
        if(is.na(quotes$upfront[i])){
            next
        }
        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("defaultcurve", 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)
}

bondhazardrate.shaped <- function(collateral, shape, R=0.4, alpha=0.25, beta=15){
    ## 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)
    sc <- new("shapedcurve", h=0.05, shape=shape, alpha=alpha, beta=beta)
    eps <- 1e-8
    counter <- 0
    while(abs(bondpv(cs, sc, R) - collateral$price/100) > eps){
        dh <- (collateral$price/100 - bondpv(cs, sc, R))/dbondpv(cs, sc, R)
        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) )
}

indexpv <- function(portfolio, index, epsilon=0){
    ## computes the intrinsic index pv of a portfolio of cds
    pl <- rep(0, length(portfolio))
    cl <- 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)
            cl[i] <- couponleg(cs, tweakedcurve, portfolio[[i]]@recovery)
            pl[i] <- defaultleg(cs, tweakedcurve, portfolio[[i]]@recovery)
        }else{
            cl[i] <- couponleg(cs, portfolio[[i]]@curve, portfolio[[i]]@recovery)
            pl[i] <- defaultleg(cs, portfolio[[i]]@curve, portfolio[[i]]@recovery)
        }
    }
    return( list(cl=mean(cl), pl=mean(pl), bp=1+mean(cl-pl)))
}

indexduration <- function(portfolio, index){
    ## compute the duration of a portfolio of survival curves
    durations <- sapply(sapply(portfolio, attr, "curve"), cdsduration, index$maturity)
    return( mean(durations) )
}

indexspread <- function(portfolio, index){
    ## 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)
    ## }
    S <- hy17$coupon-(indexpv(portfolio, index)-1)/indexduration(portfolio, index)
    return(S)
}

tweakcurves <- function(portfolio, index){
    ## computes the tweaking factor
    epsilon <- 0
    f <- function(epsilon, ...){
        abs(indexpv(portfolio, index, epsilon)$bp-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 )
}

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

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


SPmatrix <- function(portfolio, index){
    ## 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
    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,] <- SP(portfolio[[i]]@curve)[1:length(cs$dates)]
    }
    return( SP )
}

DP2 <- function(sc, dates){
    ## computes the default probability and prepay probability associated
    ## with the survival curve at the dates specified by dates
    x2T <- yearFrac(today(), dates)
    dT <- diff(c(0, x2T))
    x1T <- yearFrac(today(), 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))
    list(defaultprob = cumsum(hfun(x2T) * Qmid * dT),
         prepayprob = cumsum(pfun(x2T) * Qmid * dT))
}

getdealschedule <- function(dealdata, freq="3 months"){
    dates <- seq(dealdata$"Deal Next Pay Date", dealdata$maturity, by=freq)
    dates <- dates[dates>today()]
    return( dates )
}

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

forwardportfolioprice <- function(portfolio, startdate, rollingmaturity, coupontype, margin, recovery){
    forwardcs <- couponSchedule(nextpaydate=startdate+45, maturity=startdate+rollingmaturity,
                                      frequency="Q", "FLOAT", margin, margin, startdate=startdate)
    r <- rep(0, length(portfolio$SC))
    for(i in 1:length(portfolio$SC)){
        cl <- couponleg(forwardcs, portfolio$SC[[i]]@curve, startdate=startdate)
        pl <- contingentleg(forwardcs, portfolio$SC[[1]]@curve, portfolio$SC[[i]]@recovery, startdate=startdate)
        r[i] <- pl+cl
    }
    return(mean(r))
}