aboutsummaryrefslogtreecommitdiffstats
path: root/R/cds_functions_generic.R
diff options
context:
space:
mode:
Diffstat (limited to 'R/cds_functions_generic.R')
-rw-r--r--R/cds_functions_generic.R1800
1 files changed, 900 insertions, 900 deletions
diff --git a/R/cds_functions_generic.R b/R/cds_functions_generic.R
index 293eb4f5..31c0bd3c 100644
--- a/R/cds_functions_generic.R
+++ b/R/cds_functions_generic.R
@@ -1,900 +1,900 @@
-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=FALSE){
- stopifnot(class(startdate)=="Date")
- stopifnot(is.logical(accruedondefault))
- cs <- cs[cs$dates>=startdate,]
- T <- yearFrac(startdate, cs$dates)
- PD <- exp(- sc@h * T )
- PDadj <- if(accruedondefault){
- 0.5 * (c(1, PD[-length(PD)]) + PD)
- }else{
- 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))
- cs <- cs[cs$dates>=startdate,]
- 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, pay.at.end=FALSE){
- cs <- cs[cs$dates>=startdate,]
- T <- yearFrac(startdate, cs$dates)
- x1T <- yearFrac(startdate, sc@dates)
- hfun <- approxfun(x1T, sc@hazardrates, method="constant", rule=2)
- dT <- diff(c(0, T))
- if(!pay.at.end) {
- Q <- cumprod(exp(-hfun(T) * dT)) * cs$df
- } else {
- Q <- cumprod(exp(-hfun(T) * dT))
- }
- Qmid <- 1/2 * (c(1, Q[-length(Q)]) + Q)
- r <- (1 - recovery) * crossprod(hfun(T) * Qmid, dT)
- if(!pay.at.end) {
- return( as.numeric(r) )
- } else {
- return( as.numeric(r) * cs$df[length(cs$df)] )
- }
- })
-
-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=365/360*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
- previous.hvec <- c()
- eps <- 1e-8
- previous.cs <- data.frame()
- sc <- NULL
- for(i in 1:nrow(quotes)){
- 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 <- if(is.null(previous.hvec)) 0.05 else previous.hvec[length(previous.hvec)]
- flength <- nrow(cs) - nrow(previous.cs)
- hvec <- c(previous.hvec, rep(new.h, flength))
- sc <- new("defaultcurve", dates=cs$dates, hazardrates=hvec)
- if(is.na(quotes$upfront[i])){
- next
- }
- 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<3){
- 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, start.tweak=1L, end.tweak=NA,
- multiplicative=TRUE){
- ## 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) {
- if(is.na(end.tweak)) {
- range <- start.tweak:length(x@curve@hazardrates)
- } else {
- range <- start.tweak:end.tweak
- }
- x@curve@hazardrates[range] <- x@curve@hazardrates[range] * (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, start.tweak=1L, end.tweak=NA) {
- ## 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, start.tweak, end.tweak)
- }else {
- portfolio <- index$portfolio
- }
- startdate <- tradedate + 1
- cs <- index$cs[index$cs$unadj.dates <= maturity,]
- cl.list <- vapply(portfolio, function(x) {
- cl <- couponleg(cs, x@curve, startdate)
- if(is.na(cl)) {
- logerror(paste("couldn't compute single name coupon leg for", x@issuer))
- return( NA )
- }
- return( cl )
- }, numeric(1))
- pl.list <- vapply(portfolio, function(x) {
- pl <- defaultleg(cs, x@curve, x@recovery, startdate)
- if(is.na(pl)) {
- logerror(paste("couldn't compute single name protection leg for", x@issuer))
- return( NA )
- }
- return( pl )
- }, numeric(1))
- spread <- index$quotes$spread[index$quotes$maturity == index$maturity]
- r <- list(cl = spread * crossprod(cl.list, index$issuerweights),
- pl = crossprod(pl.list, index$issuerweights),
- bp = 1 + crossprod(index$issuerweights, 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=Sys.Date(), maturity=index$maturity) {
- ## compute the duration of a portfolio of survival curves
- cs <- index$cs[index$cs$unadj.dates <= maturity,]
- startdate <- tradedate + 1
- cl.list <- unlist(lapply(index$portfolio, function(x) couponleg(cs, x@curve, startdate)))
- return( crossprod(index$issuerweights, 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(), maturity=index$maturity) {
- if(nrow(index$cs) <= 4) {
- return( NA )
- }
- current.pv <- indexpv(index, tradedate=tradedate, maturity=maturity)$bp
- newmaturity <- maturity - lubridate::years(1)
- rolled.pv <- indexpv(index, tradedate=tradedate, maturity=newmaturity)$bp
- theta <- rolled.pv - current.pv + index$quotes$spread[index$quotes$maturity == index$maturity]
- 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){
- start.tweak <- 1
- range <- which(index$cs$unadj.dates<=index$quotes$maturity[i])
- end.tweak <- range[length(range)]
- }else{
- range <- which((index$cs$unadj.dates > index$quotes$maturity[i-1]) &
- (index$cs$unadj.dates <= index$quotes$maturity[i]))
- start.tweak <- range[1]
- end.tweak <- range[length(range)]
- }
- f <- function(epsilon, index, tradedate, start.tweak, end.tweak, i){
- abs(indexpv(index, epsilon, tradedate, maturity=index$quotes$maturity[i],
- TRUE, start.tweak, end.tweak)$bp - index$quotes$price[i])
- }
- epsilon[i] <- optimize(f, c(-0.25, 0.30),
- index, tradedate, start.tweak, end.tweak, i, tol=1e-6)$minimum
- index$portfolio <- tweakportfolio(index$portfolio, epsilon[i], start.tweak, end.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 flat
- 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{
- logprob <- logprob - curve@hazardrates[i] * (Tmat - T[i])
- break
- }
- }
- return( exp(as.numeric(logprob)) )
-}
-
-FEP <- function(index, exerciseDate) {
- ## computes the front end protection for a credit index at a given date
- DP <- 1-unlist(lapply(index$portfolio, survivalProbability.exact, date=exerciseDate))
- R <- unlist(lapply(index$portfolio, function(x)x@recovery))
- fep <- crossprod(index$issuerweights, (1-R)*DP)
- return( fep )
-}
-
-SP <- function(sc, dates, dT){
- ## computes the survival probability associated with the survival curve
- h <- approx(sc@dates, sc@hazardrates, dates, "constant", rule=2)$y
- return( cumprod(exp(-h * dT)) )
-}
-
-SPmatrix <- function(portfolio, 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), length(dates))
- for(i in 1:length(portfolio)){
- T <- c(0, yearFrac(portfolio[[i]]@startdate, dates))
- dT <- diff(T)
- SP[i,] <- SP(portfolio[[i]]@curve, dates, dT)
- }
- colnames(SP) <- as.character(dates)
- rownames(SP) <- unlist(lapply(portfolio, function(x)x@issuer))
- 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=c("Quarterly", "Monthly"),
- startdate=Sys.Date()){
- ## computes the default and prepay probability matrix of a portfolio
- ## at the dates specified from dealdata
- freq <- match.arg(freq)
- 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, startdate)
- 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")
- }
-}
-
-forward.cl <- function(index, exerciseDate){
- aux <- function(creditCurve){
- couponleg(index$cs, creditCurve@curve, startdate=exerciseDate)
- }
- return( crossprod(index$issuerweights, unlist(lapply(index$portfolio, aux))) )
-}
-
-forward.prot <- function(index, exerciseDate){
- prot <- unlist(lapply(index$portfolio, function(cc){
- defaultleg(index$cs, cc@curve, cc@recovery, startdate=exerciseDate)
- }))
- return( crossprod(prot, index$issuerweights) )
-}
-
-defaultAdjustedForwardIndexPrice <- function(index, exerciseDate, fixedRate=0.05){
- tes <- addBusDay(exerciseDate)
- df <- YC$discount(tes)
- price <- 1 - FEP(index, exerciseDate) +
- 1/df * (forward.cl(index, exerciseDate) * fixedRate -
- forward.prot(index, exerciseDate) - cdsAccrued(exerciseDate, fixedRate))
- return( price )
-}
-
-forwardflatcds <- function(h, cs, tradeDate, exerciseDate, fixedRate=0.05, R=0.4){
- tes <- addBusDay(exerciseDate)
- fep <- (1-R)*(1-exp(-h*yearFrac(tradeDate, exerciseDate)))
- df <- YC$discount(tes)
- sc <- new("flatcurve", h=h)
- cl <- couponleg(cs, sc, startdate=exerciseDate)
- pl <- defaultleg(cs, sc, recovery=R, startdate=exerciseDate)
- price <- 1 - fep +
- 1/df*(cl*fixedRate-pl-cdsAccrued(exerciseDate, fixedRate))
- return( price )
-}
+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=FALSE){
+ stopifnot(class(startdate)=="Date")
+ stopifnot(is.logical(accruedondefault))
+ cs <- cs[cs$dates>=startdate,]
+ T <- yearFrac(startdate, cs$dates)
+ PD <- exp(- sc@h * T )
+ PDadj <- if(accruedondefault){
+ 0.5 * (c(1, PD[-length(PD)]) + PD)
+ }else{
+ 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))
+ cs <- cs[cs$dates>=startdate,]
+ 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, pay.at.end=FALSE){
+ cs <- cs[cs$dates>=startdate,]
+ T <- yearFrac(startdate, cs$dates)
+ x1T <- yearFrac(startdate, sc@dates)
+ hfun <- approxfun(x1T, sc@hazardrates, method="constant", rule=2)
+ dT <- diff(c(0, T))
+ if(!pay.at.end) {
+ Q <- cumprod(exp(-hfun(T) * dT)) * cs$df
+ } else {
+ Q <- cumprod(exp(-hfun(T) * dT))
+ }
+ Qmid <- 1/2 * (c(1, Q[-length(Q)]) + Q)
+ r <- (1 - recovery) * crossprod(hfun(T) * Qmid, dT)
+ if(!pay.at.end) {
+ return( as.numeric(r) )
+ } else {
+ return( as.numeric(r) * cs$df[length(cs$df)] )
+ }
+ })
+
+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=365/360*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
+ previous.hvec <- c()
+ eps <- 1e-8
+ previous.cs <- data.frame()
+ sc <- NULL
+ for(i in 1:nrow(quotes)){
+ 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 <- if(is.null(previous.hvec)) 0.05 else previous.hvec[length(previous.hvec)]
+ flength <- nrow(cs) - nrow(previous.cs)
+ hvec <- c(previous.hvec, rep(new.h, flength))
+ sc <- new("defaultcurve", dates=cs$dates, hazardrates=hvec)
+ if(is.na(quotes$upfront[i])){
+ next
+ }
+ 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<3){
+ 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, start.tweak=1L, end.tweak=NA,
+ multiplicative=TRUE){
+ ## 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) {
+ if(is.na(end.tweak)) {
+ range <- start.tweak:length(x@curve@hazardrates)
+ } else {
+ range <- start.tweak:end.tweak
+ }
+ x@curve@hazardrates[range] <- x@curve@hazardrates[range] * (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, start.tweak=1L, end.tweak=NA) {
+ ## 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, start.tweak, end.tweak)
+ }else {
+ portfolio <- index$portfolio
+ }
+ startdate <- tradedate + 1
+ cs <- index$cs[index$cs$unadj.dates <= maturity,]
+ cl.list <- vapply(portfolio, function(x) {
+ cl <- couponleg(cs, x@curve, startdate)
+ if(is.na(cl)) {
+ logerror(paste("couldn't compute single name coupon leg for", x@issuer))
+ return( NA )
+ }
+ return( cl )
+ }, numeric(1))
+ pl.list <- vapply(portfolio, function(x) {
+ pl <- defaultleg(cs, x@curve, x@recovery, startdate)
+ if(is.na(pl)) {
+ logerror(paste("couldn't compute single name protection leg for", x@issuer))
+ return( NA )
+ }
+ return( pl )
+ }, numeric(1))
+ spread <- index$quotes$spread[index$quotes$maturity == index$maturity]
+ r <- list(cl = spread * crossprod(cl.list, index$issuerweights),
+ pl = crossprod(pl.list, index$issuerweights),
+ bp = 1 + crossprod(index$issuerweights, 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=Sys.Date(), maturity=index$maturity) {
+ ## compute the duration of a portfolio of survival curves
+ cs <- index$cs[index$cs$unadj.dates <= maturity,]
+ startdate <- tradedate + 1
+ cl.list <- unlist(lapply(index$portfolio, function(x) couponleg(cs, x@curve, startdate)))
+ return( crossprod(index$issuerweights, 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(), maturity=index$maturity) {
+ if(nrow(index$cs) <= 4) {
+ return( NA )
+ }
+ current.pv <- indexpv(index, tradedate=tradedate, maturity=maturity)$bp
+ newmaturity <- maturity - lubridate::years(1)
+ rolled.pv <- indexpv(index, tradedate=tradedate, maturity=newmaturity)$bp
+ theta <- rolled.pv - current.pv + index$quotes$spread[index$quotes$maturity == index$maturity]
+ 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){
+ start.tweak <- 1
+ range <- which(index$cs$unadj.dates<=index$quotes$maturity[i])
+ end.tweak <- range[length(range)]
+ }else{
+ range <- which((index$cs$unadj.dates > index$quotes$maturity[i-1]) &
+ (index$cs$unadj.dates <= index$quotes$maturity[i]))
+ start.tweak <- range[1]
+ end.tweak <- range[length(range)]
+ }
+ f <- function(epsilon, index, tradedate, start.tweak, end.tweak, i){
+ abs(indexpv(index, epsilon, tradedate, maturity=index$quotes$maturity[i],
+ TRUE, start.tweak, end.tweak)$bp - index$quotes$price[i])
+ }
+ epsilon[i] <- optimize(f, c(-0.25, 0.30),
+ index, tradedate, start.tweak, end.tweak, i, tol=1e-6)$minimum
+ index$portfolio <- tweakportfolio(index$portfolio, epsilon[i], start.tweak, end.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 flat
+ 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{
+ logprob <- logprob - curve@hazardrates[i] * (Tmat - T[i])
+ break
+ }
+ }
+ return( exp(as.numeric(logprob)) )
+}
+
+FEP <- function(index, exerciseDate) {
+ ## computes the front end protection for a credit index at a given date
+ DP <- 1-unlist(lapply(index$portfolio, survivalProbability.exact, date=exerciseDate))
+ R <- unlist(lapply(index$portfolio, function(x)x@recovery))
+ fep <- crossprod(index$issuerweights, (1-R)*DP)
+ return( fep )
+}
+
+SP <- function(sc, dates, dT){
+ ## computes the survival probability associated with the survival curve
+ h <- approx(sc@dates, sc@hazardrates, dates, "constant", rule=2)$y
+ return( cumprod(exp(-h * dT)) )
+}
+
+SPmatrix <- function(portfolio, 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), length(dates))
+ for(i in 1:length(portfolio)){
+ T <- c(0, yearFrac(portfolio[[i]]@startdate, dates))
+ dT <- diff(T)
+ SP[i,] <- SP(portfolio[[i]]@curve, dates, dT)
+ }
+ colnames(SP) <- as.character(dates)
+ rownames(SP) <- unlist(lapply(portfolio, function(x)x@issuer))
+ 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=c("Quarterly", "Monthly"),
+ startdate=Sys.Date()){
+ ## computes the default and prepay probability matrix of a portfolio
+ ## at the dates specified from dealdata
+ freq <- match.arg(freq)
+ 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, startdate)
+ 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")
+ }
+}
+
+forward.cl <- function(index, exerciseDate){
+ aux <- function(creditCurve){
+ couponleg(index$cs, creditCurve@curve, startdate=exerciseDate)
+ }
+ return( crossprod(index$issuerweights, unlist(lapply(index$portfolio, aux))) )
+}
+
+forward.prot <- function(index, exerciseDate){
+ prot <- unlist(lapply(index$portfolio, function(cc){
+ defaultleg(index$cs, cc@curve, cc@recovery, startdate=exerciseDate)
+ }))
+ return( crossprod(prot, index$issuerweights) )
+}
+
+defaultAdjustedForwardIndexPrice <- function(index, exerciseDate, fixedRate=0.05){
+ tes <- addBusDay(exerciseDate)
+ df <- YC$discount(tes)
+ price <- 1 - FEP(index, exerciseDate) +
+ 1/df * (forward.cl(index, exerciseDate) * fixedRate -
+ forward.prot(index, exerciseDate) - cdsAccrued(exerciseDate, fixedRate))
+ return( price )
+}
+
+forwardflatcds <- function(h, cs, tradeDate, exerciseDate, fixedRate=0.05, R=0.4){
+ tes <- addBusDay(exerciseDate)
+ fep <- (1-R)*(1-exp(-h*yearFrac(tradeDate, exerciseDate)))
+ df <- YC$discount(tes)
+ sc <- new("flatcurve", h=h)
+ cl <- couponleg(cs, sc, startdate=exerciseDate)
+ pl <- defaultleg(cs, sc, recovery=R, startdate=exerciseDate)
+ price <- 1 - fep +
+ 1/df*(cl*fixedRate-pl-cdsAccrued(exerciseDate, fixedRate))
+ return( price )
+}
generated by cgit v1.2.3-70-g09d2 (git 2.52.0) at 2026-01-03 08:20:51 +0000