diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/bandits.R | 112 | ||||
| -rw-r--r-- | R/build_SC.R | 140 | ||||
| -rw-r--r-- | R/calibrate_tranches.R | 274 | ||||
| -rw-r--r-- | R/cds_functions.R | 343 | ||||
| -rw-r--r-- | R/cds_functions_generic.R | 664 | ||||
| -rw-r--r-- | R/cds_utils.R | 146 | ||||
| -rw-r--r-- | R/clopricer.R | 87 | ||||
| -rw-r--r-- | R/deal_pricer.R | 22 | ||||
| -rw-r--r-- | R/dumpdata.R | 4 | ||||
| -rw-r--r-- | R/ema.R | 3 | ||||
| -rw-r--r-- | R/etdb.R | 3 | ||||
| -rw-r--r-- | R/filterruns.R | 25 | ||||
| -rw-r--r-- | R/index_definitions.R | 24 | ||||
| -rw-r--r-- | R/interpweights.R | 94 | ||||
| -rw-r--r-- | R/interpweights_2.R | 56 | ||||
| -rw-r--r-- | R/intex_deals_functions.R | 253 | ||||
| -rw-r--r-- | R/l1tf.R | 113 | ||||
| -rw-r--r-- | R/latestprices.R | 44 | ||||
| -rw-r--r-- | R/load_bloomberg_data.R | 49 | ||||
| -rw-r--r-- | R/loadcashflows.R | 38 | ||||
| -rw-r--r-- | R/loan_universe.R | 70 | ||||
| -rw-r--r-- | R/lossdistrib.c | 374 | ||||
| -rw-r--r-- | R/mapping.R | 140 | ||||
| -rw-r--r-- | R/mapping_fast.R | 62 | ||||
| -rw-r--r-- | R/optimization.R | 440 | ||||
| -rw-r--r-- | R/parse_intex.R | 18 | ||||
| -rw-r--r-- | R/plot_distributions.R | 35 | ||||
| -rw-r--r-- | R/reinvestingdistribution.R | 56 | ||||
| -rw-r--r-- | R/test.cds_functions.R | 46 | ||||
| -rw-r--r-- | R/time_of_default.R | 99 | ||||
| -rw-r--r-- | R/tranche_functions.R | 798 | ||||
| -rw-r--r-- | R/transactions.R | 26 | ||||
| -rw-r--r-- | R/yieldcurve.R | 63 |
33 files changed, 4721 insertions, 0 deletions
diff --git a/R/bandits.R b/R/bandits.R new file mode 100644 index 00000000..ee547e3a --- /dev/null +++ b/R/bandits.R @@ -0,0 +1,112 @@ +require(RBloomberg)
+conn <- blpConnect()
+
+data <- c()
+for(ticker in sp500.tickers.extended){
+ data <- cbind(data,blpGetData(conn,paste(ticker,"US Equity"),"PX_LAST",start=as.chron("01/01/2000")))
+}
+
+spy <- blpGetData(conn,"spy US Equity","PX_LAST",start=as.chron("01/01/2000"))
+
+spdr <- c("XLY","XLP","XLE","XLF","XLV","XLI","XLB","XLK","XLU")
+data <- c()
+for(ticker in spdr){
+ data <- cbind(data,blpGetData(conn,paste(ticker,"US Equity"),"PX_LAST",start=as.chron("01/01/2000")))
+}
+colnames(data) <- spdr
+
+n.stocks <- ncol(P)
+N <- nrow(P)
+current.wealth <- 1
+w <- rep(1/n.stocks,n.stocks)
+dP <- apply(P,2,diff)
+L <- rep(0,n.stocks)
+V <- 0
+W <- w
+pnl <- 0
+for(i in 1:(N-1)){
+ r <- dP[i,]/as.numeric(P[i,])
+ #r <- c(r,-r)
+ pnl <- cbind(pnl,current.wealth*crossprod(w,r))
+ current.wealth <- 1+sum(pnl)
+ L <- rbind(L,-r)
+ V <- V+crossprod(w,r^2)
+ T <- 2/3*sqrt(log(N)/V)
+ #w <- exp(-T*(apply(1+L,2,prod)-1))
+ w <- exp(-T*colSums(L))
+ w <- w/sum(w)
+ W <- rbind(W,w)
+ if(i%%10==0){
+ cat(current.wealth,sep="\n")
+ }
+}
+price2return <- function(x){
+ diff(x)/x[-length(x)]
+}
+#number of shares implementations
+tc <- 0.005+0.02
+days <- nrow(data.bus)
+init.capital <- 1000000
+tickers <- memb(sp500.tickers,add,as.Date(time(data.bus)[1]))
+tickers.index <- which(sp500.tickers.extended%in%tickers)
+n.stocks <- length(tickers)
+w <- rep(0,length(sp500.tickers.extended))
+w[tickers.index] <- rep(1/n.stocks,n.stocks)
+N <- round((capital*w)/as.numeric(data.bus[1,tickers.index]))
+dP <- apply(data.bus,2,diff)
+L <- rep(0,n.stocks)
+V <- 0
+W <- w
+pnl <- 0
+tcvec <- sum(N)*tc
+for(i in 1:days){
+ tickers <- memb(sp500.tickers,add,as.Date(time(data.bus)[1]))
+ tickers.index <- which(sp500.tickers.extended%in%tickers)
+ r <- dP[i,]/as.numeric(data.bus[i,])
+ pnl <- cbind(pnl,crossprod(N,dP[i,]))
+ capital <- init.capital+sum(pnl)-sum(tcvec)
+ L <- rbind(L,-r)
+ V <- V+crossprod(w,r^2)
+ T <- 2/3*sqrt(log(days)/V)
+ #w <- exp(-T*(apply(1+L,2,prod)-1))
+ w <- exp(-T*colSums(L))
+ w <- w/sum(w)
+ newN <- round((capital*w)/as.numeric(data.bus[(i+1),tickers.index]))
+ tcvec <- c(tcvec,sum(abs(newN-N)*tc))
+ N <- newN
+ if(i%%10==0){
+ cat(capital,sep="\n")
+ }
+}
+
+fixed.rebal <- function(P,delta){
+ init.capital <- 1
+ dP <- apply(P,2,diff)
+ capital <- init.capital
+ pnl <- c()
+ for(i in 1:(nrow(P)-1)){
+ r <- dP[i,]/P[i,]
+ pnl <- c(pnl,capital*crossprod(delta,r))
+ capital <- init.capital+sum(pnl)
+ }
+ return( pnl )
+}
+memb <- function(index,add,date){
+ #return the list of index constituents
+ toreverse <- add[add$Date>=date,]
+ current.index <- index
+ for(i in 1:nrow(toreverse)){
+ if(!is.na(toreverse$Add[i])){
+ current.index <- current.index[-which(current.index==toreverse$Add[i])]
+ }
+ if(!is.na(toreverse$Del[i])){
+ current.index <- sort(c(current.index,toreverse$Del[i]))
+ }
+ }
+ current.index
+}
+
+tickerslist <- c()
+for(i in 101:1000){
+ tickerslist <- rbind(tickerslist,memb(sp500.tickers,add,as.Date(time(data)[i])))
+}
diff --git a/R/build_SC.R b/R/build_SC.R new file mode 100644 index 00000000..11e2ec7b --- /dev/null +++ b/R/build_SC.R @@ -0,0 +1,140 @@ +library("RQuantLib")
+library("parallel")
+root = "//WDSENTINEL/share/CorpCDOs"
+source(file.path(root, "R", "intex_deals_functions.R"))
+
+MarkitData <- getMarkitIRData()
+L1m <- buildMarkitYC(MarkitData, dt = 1/12)
+L2m <- buildMarkitYC(MarkitData, dt = 1/6)
+L3m <- buildMarkitYC(MarkitData)
+L6m <- buildMarkitYC(MarkitData, dt = 1/2)
+L12m <- buildMarkitYC(MarkitData, dt = 1)
+setEvaluationDate(as.Date(MarkitData$effectiveasof))
+today <- function(){
+ return(as.Date(MarkitData$effectiveasof))
+}
+
+bps <- 1e-4
+global.params <- list()
+global.params$recovery.assumptions <- list("Loan"=0.7,
+ "SecondLien"=0.3,
+ "Bond"=0.4,
+ "Mezzanine"=0.15,
+ "Adj_Covlite"=0.1)
+
+global.params$cdoprices <- list("Aaa"=90,
+ "Aa"=80,
+ "A"=70,
+ "Baa"=60,
+ "Ba"=50,
+ "B"=40,
+ "NR"=40)
+#reinvest in 7 years assets
+global.params$rollingmaturity <- 7 * 365
+global.params$defaultedlag <- 90
+global.params$defaultcorr <- 0.4
+global.params$defaultbondhazardrate <- 500 * bps
+global.params$defaultloanhazardrate <- 400 * bps
+global.params$alpha <- 0.25
+global.params$beta <- 15
+global.params$shape <- function(T)0.25+(1-exp(-T/5))
+
+## dealnames <- c("stonln1", "babs072", "symph4", "flags5", "cent11", "wasatl", "oceant2", "acacl071", "limes")
+dealnames <- c("acacl071", "limes")
+calibration.date <- "2012-11-08"
+calibration <- read.table(file.path(root, "Scenarios", paste0("calibration-",calibration.date,".csv")),
+ sep=",", header=T)
+Z <- calibration$Z
+w <- calibration$w
+
+rho <- 0.45
+Ngrid <- 201
+cl <- makeCluster(6)
+clusterExport(cl, list("shockprob", "lossdistC.prepay.joint", "lossrecovdist.joint.term", "lossdistC.joint"))
+support <- seq(0, 1, length=Ngrid)
+for(deal.name in dealnames){
+ deal.portfolio <- buildSC.portfolio(deal.name, global.params)
+ deal.data <- getdealdata(deal.name)
+ A <- SPmatrix2(deal.portfolio$SC, deal.data, freq="3 months")
+ S <- 1 - sapply(deal.portfolio$SC, attr, "recov")
+ deal.weights <- deal.portfolio$notional/sum(deal.portfolio$notional)
+ clusterExport(cl, list("deal.weights"))
+
+ deal.dates <- getdealschedule(deal.data)
+ ## compute reinvestment price
+ reinvloanprice <- rep(0, length(deal.dates))
+ reinvmezzprice <- rep(0, length(deal.dates))
+ for(i in 1:length(deal.dates)){
+ reinvloanprice[i] <- forwardportfolioprice(deal.portfolio, deal.dates[i], global.params$rollingmaturity, "FLOAT", 0.0235)
+ reinvmezzprice[i] <- forwardportfolioprice(deal.portfolio, deal.dates[i], global.params$rollingmaturity, "FLOAT", 0.0385)
+ }
+
+ dp <- A$DP
+ pp <- A$PP
+ Smat <- matrix(S, length(deal.weights), ncol(dp))
+ ##no correlation setup
+ ##test <- lossrecovdist.joint.term(dp, pp, issuerweights, Smat, Ngrid)
+
+ dpmod <- MFupdate.prob(Z, w, rho, dp)
+ ppmod <- MFupdate.prob(-Z, w, rho, pp)
+ ## dist <- MFlossrecovdist.prepay(w, Z, rho, dp, dpmod, pp, ppmod, deal.weights, 1-S, Ngrid, TRUE)
+ dist.joint <- MFlossdist.prepay.joint(cl, w, Z, rho, dp, dpmod,
+ pp, ppmod, deal.weights, 1-S, Ngrid, FALSE, n.chunks=4)
+ clusterCall(cl, gc)
+ gc()
+ ## if don't want to use cluster (use less memory)
+ ## dist.joint <- MFlossdist.prepay.joint(NULL, w, Z, rho, dp, dpmod,
+ ## pp, ppmod, issuerweights, 1-S, Ngrid=201, FALSE)
+
+ ## two ways to compute the joint (D, R) distribution
+ ## first using the actual function (numerically instable)
+ ## dist.joint2 <- MFlossdist.prepay.joint(cl, w, Z, rho, dp, dpmod,
+ ## pp, ppmod, issuerweights, 1-S, Ngrid=201, TRUE)
+
+ ## second, by doing a change of variable seems to work better for now
+ distDR <- dist.transform(dist.joint)
+
+ ## compute E(R|D)
+ R <- matrix(0, Ngrid, ncol(dp))
+ for(t in 1:ncol(dp)){
+ R[,t] <- (sweep(distDR[t,,], 1, rowSums(distDR[t,,]), "/") %*% support)/support
+ }
+ R[1,] <- 0
+ n.scenarios <- 100
+ percentiles <- (seq(0, 1, length=n.scenarios+1)[-1]+
+ seq(0, 1, length=n.scenarios+1)[-(n.scenarios+1)])/2
+
+ ## compute scenariosd
+ scenariosd <- matrix(0, n.scenarios, ncol(dp))
+ for(t in 1:ncol(dp)){
+ D <- rowSums(distDR[t,,])
+ D <- D/sum(D)
+ Dfun <- splinefun(c(0, cumsum(D)), c(0, support), "monoH.FC")
+ ## dvallow <- floor(Dfun(percentiles)*(Ngrid-1))
+ ## dvalup <- ceil(Dfun(percentiles)*(Ngrid-1))
+ scenariosd[,t] <- Dfun(percentiles)
+ }
+
+ ## compute scenariosr
+ scenariosr <- matrix(0, n.scenarios, ncol(dp))
+ for(t in 1:ncol(dp)){
+ Rfun <- approxfun(support, R[,t], rule=2)
+ scenariosr[,t] <- Rfun(scenariosd[,t])
+ }
+
+ cdr <- cdrfromscenarios(scenariosd, deal.dates)
+ intexrecov <- recoveryfromscenarios(scenariosd, scenariosr)
+
+ ## linear approximation for monthly scenarios
+ deal.datesmonthly <- getdealschedule(deal.data, "1 month")
+ cdrmonthly <- matrix(0, n.scenarios, length(deal.datesmonthly))
+ recoverymonthly <- matrix(0, n.scenarios, length(deal.datesmonthly))
+ for(i in 1:n.scenarios){
+ cdrmonthly[i,] <- approx(deal.dates, cdr[i,], deal.datesmonthly, rule=2)$y
+ recoverymonthly[i,] <- approx(deal.dates, intexrecov[i,], deal.datesmonthly, rule=2)$y
+ }
+
+ write.table(cdrmonthly, file=file.path(root, "Scenarios", paste0(deal.name,"-cdr.csv")), row.names=F, col.names=F, sep=",")
+ write.table(recoverymonthly * 100, file=file.path(root, "Scenarios", paste0(deal.name,"-recovery.csv")), row.names=F, col.names=F, sep=",")
+ cat("generated scenarios for:", deal.name, "\n")
+}
diff --git a/R/calibrate_tranches.R b/R/calibrate_tranches.R new file mode 100644 index 00000000..3e2599fc --- /dev/null +++ b/R/calibrate_tranches.R @@ -0,0 +1,274 @@ +library("parallel")
+setwd("//WDSENTINEL/share/CorpCDOs/R")
+source("cds_utils.R")
+source("cds_functions_generic.R")
+source("index_definitions.R")
+source("tranche_functions.R")
+source("yieldcurve.R")
+source("optimization.R")
+
+cl <- makeCluster(6)
+
+MarkitData <- getMarkitIRData()
+L1m <- buildMarkitYC(MarkitData, dt = 1/12)
+L2m <- buildMarkitYC(MarkitData, dt = 1/6)
+L3m <- buildMarkitYC(MarkitData)
+L6m <- buildMarkitYC(MarkitData, dt = 1/2)
+setEvaluationDate(as.Date(MarkitData$effectiveasof))
+
+## calibrate HY19
+## calibrate the single names curves
+singlenames.data <- read.table(file="clipboard", sep="\t", header=T)
+nondefaulted <- singlenames.data[!singlenames.data$ticker %in% hy17$defaulted,]
+bps <- 1e-4
+
+cdsdates <- as.Date(character(0))
+for(tenor in paste0(1:5, "y")){
+ cdsdates <- c(cdsdates, cdsMaturity(tenor))
+}
+
+## clusterEvalQ(cl, {setClass("abstractcurve")
+## setClass("defaultcurve", contains="abstractcurve",
+## representation(dates="Date", hazardrates="numeric"))
+## setClass("creditcurve", representation(issuer="character", startdate="Date",
+## recovery="numeric", curve="defaultcurve"))})
+## clusterExport(cl, list("nondefaulted", "cdsdates", "cdshazardrate", "today",
+## "bps", "couponSchedule", "nextIMMDate", "DiscountCurve", "L3m"))
+## test <- parSapply(cl, 1:nrow(nondefaulted), parf)
+
+## parf <- function(i){
+## SC <- new("creditcurve",
+## recovery=nondefaulted$recovery[i]/100,
+## startdate=today(),
+## issuer=as.character(nondefaulted$ticker[i]))
+## quotes <- data.frame(maturity=cdsdates, upfront = as.numeric(nondefaulted[i,5:9])*0.01,
+## running = rep(nondefaulted$running[i]*bps,5))
+## return( cdshazardrate(quotes, nondefaulted$recovery[i]/100))
+## }
+
+hy19portfolio <- c()
+for(i in 1:nrow(nondefaulted)){
+ SC <- new("creditcurve",
+ recovery=nondefaulted$recovery[i]/100,
+ startdate=today(),
+ issuer=as.character(nondefaulted$ticker[i]))
+ quotes <- data.frame(maturity=cdsdates, upfront = as.numeric(nondefaulted[i,4:8])*0.01,
+ running=rep(nondefaulted$running[i]*bps, 5))
+ SC@curve <- cdshazardrate(quotes, nondefaulted$recovery[i]/100)
+ hy19portfolio <- c(hy19portfolio, SC)
+}
+
+issuerweights <- rep(1/length(hy19portfolio), length(hy19portfolio))
+hy19$indexref <- 0.99
+hy19portfolio.tweaked <- tweakcurves(hy19portfolio, hy19)
+
+SurvProb <- SPmatrix(hy19portfolio.tweaked, hy19)
+## load common parameters
+K <- c(0, 0.15, 0.25, 0.35, 1)
+Kmodified <- adjust.attachments(K, hy19$loss, hy19$factor)
+tranche.upf <- c(37.5625, 87.25, 101.8125, 115.0625)
+tranche.running <- c(0.05, 0.05, 0.05, 0.05)
+
+Ngrid <- 2*nrow(nondefaulted)+1
+recov <- sapply(hy19portfolio.tweaked, attr, "recovery")
+cs <- couponSchedule(nextIMMDate(today()), hy19$maturity,"Q", "FIXED", 0.05, 0)
+
+## calibrate the tranches using base correlation
+rhovec <- c()
+f <- function(rho, ...){
+ temp <- BClossdistC(SurvProb, issuerweights, recov, rho, Ngrid)
+ bp <- 100*(1+1/(Kmodified[i]-Kmodified[i-1]) *
+ (tranche.pv(temp$L, temp$R, cs, 0, Kmodified[i], Ngrid) -
+ tranche.pv(oldtemp$L, oldtemp$R, cs, 0, Kmodified[i-1], Ngrid)))
+ return( abs(tranche.upf[i-1]-bp))
+}
+
+for(i in 2:length(Kmodified)){
+ rho <- optimize(f, interval=c(0,1),
+ SurvProb, issuerweights, recov, Ngrid, tranche.upf, Kmodified, cs, oldtemp)$minimum
+ oldtemp <- BClossdistC(SurvProb, issuerweights, recov, rho, Ngrid)
+ rhovec <- c(rhovec, rho)
+}
+
+rhovec <- c(0, rhovec)
+deltas <- c()
+for(i in 2:5){
+ deltas <- c(deltas, BCtranche.delta(hy19portfolio.tweaked, hy19, 0.05, K[i-1], K[i], rhovec[i-1], rhovec[i], Ngrid))
+}
+
+##calibrate by modifying the factor distribution
+bottomup <- 1:3
+topdown <- 2:4
+n.int <- 500
+n.credit <- length(hy19portfolio)
+errvec <- c()
+quadrature <- gauss.quad.prob(n.int, "normal")
+w <- quadrature$weights
+Z <- quadrature$nodes
+w.mod <- w
+defaultprob <- 1 - SurvProb
+p <- defaultprob
+rho <- 0.45
+
+clusterExport(cl, list("shockprob", "issuerweights", "rho", "Z", "lossrecovdist.term",
+ "lossrecovdist", "lossdistC", "Ngrid",
+ "tranche.pvvec", "tranche.pv", "tranche.pl", "tranche.cl",
+ "trancheloss", "trancherecov", "pos", "Kmodified", "cs"))
+## TODO: investigate if this is the right thing w.r.t recovery
+parf <- function(i){
+ pshocked <- apply(p, 2, shockprob, rho=rho, Z=Z[i])
+ S <- 1 - Rstoch[i,,]
+ dist <- lossrecovdist.term(pshocked,, issuerweights, S, Ngrid)
+ return( tranche.pvvec(Kmodified, dist$L, dist$R, cs))
+}
+
+for(l in 1:100){
+ Rstoch <- array(0, dim=c(n.int, n.credit, ncol(SurvProb)))
+ for(t in 1:ncol(SurvProb)){
+ for(i in 1:n.credit){
+ Rstoch[,i,t] <- stochasticrecov(recov[i], 0, Z, w.mod, rho, defaultprob[i,t], p[i,t])
+ }
+ }
+
+ clusterExport(cl, list("Rstoch", "p"))
+ result <- parSapply(cl, 1:n.int, parf)
+ ## solve the optimization problem
+ program <- KLfit(100*(result[bottomup,]+1), w, tranche.upf[bottomup])
+
+
+ err <- 0
+ for(i in 1:n.credit){
+ for(j in 1:ncol(p)){
+ err <- err + abs(crossprod(shockprob(p[i,j], rho, Z), program$weight) - defaultprob[i,j])
+ }
+ }
+ errvec <- c(errvec, err)
+
+ ## update the new probabilities
+ p <- MFupdate.prob(Z, program$weight, rho, defaultprob)
+
+ errvec <- c(errvec, err)
+ w.mod <- program$weight
+ cat(err,"\n")
+}
+
+write.table(data.frame(Z=Z, w=w.mod), file=paste0("calibration-", Sys.Date(), ".csv"), col.names=T, row.names=F, sep=",")
+
+## computes MFdeltas
+newportf <- hy19portfolio.tweaked
+eps <- 1e-4
+for(i in 1:length(newportf)){
+ newportf[[i]]@curve@hazardrates <- hy19portfolio.tweaked[[i]]@curve@hazardrates * (1 + eps)
+}
+SurvProb2 <- SPmatrix(newportf, hy19)
+p2 <- MFupdate.prob(Z, w.mod, rho, 1-SurvProb2)
+dPVtranches <- MFtranche.pv(cl, cs, w.mod, rho, 1-SurvProb2, p2, issuerweights) - MFtranche.pv(cl, cs, w.mod, rho, defaultprob, p, issuerweights)
+dPVindex <- indexpv(newportf, hy19)-indexpv(hy19portfolio.tweaked, hy19)
+MFdeltas <- dPVtranches/dPVindex
+
+#global deltas
+PVtranches <- MFtranche.pv(cl, cs, w.mod, rho, 1-SurvProb, p, issuerweights, Kmodified)
+pnlindex <- 1+t(PVtranches$pv.w)%*%diff(Kmodified)-hy19$indexref
+plot(1-cumsum(w.mod),PVindex, main="pnl of going long the index", xlab="Market factor cumulative probability distribution", type="l")
+pnltranches <- sweep(1+t(PVtranches$pv.w), 2, tranche.upf/100)
+matplot(1-cumsum(w.mod), pnltranches, 2, tranche.upf/100, main="pnl of going long the tranches", xlab="Market factor cumulative probability distribution", type="l")
+global.deltas <- rep(0,4)
+for(i in 1:4){
+ global.deltas[i] <- lm(pnltranches[,i]~0+pnlindex, weights=w.mod)$coef
+}
+pnlhedged.model <- pnltranches-pnlindex%*%t(global.deltas)
+pnlhedged.market <- pnltranches-pnlindex%*%t(market.deltas)
+## generate a bunch of plots
+postscript("hedged tranches global.eps")
+matplot(1-cumsum(w.mod), pnlhedged.model, type="l", ylab="pnl of the hedged package (global deltas)", xlab="Market factor cumulative probability distribution")
+legend(0.35, y=0.75, c("0-15","15-25","25-35","35-100"), col=1:4, lty=1:4)
+dev.off()
+postscript("hedged tranches market.eps")
+matplot(1-cumsum(w.mod), pnlhedged.market, type="l", ylab="Pnl", main="pnl of the hedged package (market deltas)", xlab="Market factor cumulative probability distribution")
+legend(0.15, y=1.15, c("0-15","15-25","25-35","35-100"), col=1:4, lty=1:4)
+dev.off()
+postscript("0-15 hedged.eps")
+matplot(1-cumsum(w.mod), cbind(pnlhedged.model[,1], pnlhedged.market[,1]), type="l", ylab="Pnl", main="Pnl of the hedged package (market vs global deltas)\n 0-15 tranche", xlab="Market factor cumulative probability distribution")
+dev.off()
+postscript("15-25 hedged.eps")
+matplot(1-cumsum(w.mod), cbind(pnlhedged.model[,2], pnlhedged.market[,2]), type="l", ylab="Pnl", main="Pnl of the hedged package (market vs global deltas)\n 15-25 tranche", xlab="Market factor cumulative probability distribution")
+dev.off()
+postscript("25-35 hedged.eps")
+matplot(1-cumsum(w.mod), cbind(pnlhedged.model[,3], pnlhedged.market[,3]), type="l", ylab="Pnl", main="Pnl of the hedged package (market vs global deltas)\n 25-35 tranche", xlab="Market factor cumulative probability distribution")
+dev.off()
+postscript("35-100 hedged.eps")
+matplot(1-cumsum(w.mod), cbind(pnlhedged.model[,4], pnlhedged.market[,4]), type="l", ylab="Pnl", main="Pnl of the hedged package (market vs global deltas)\n 35-100 tranche", xlab="Market factor cumulative probability distribution")
+dev.off()
+
+## scenario based pricing
+## generate the scenarios by finding the quantiles of the loss and recovery distributions
+n.scenarios <- 100
+percentiles <- (seq(0, 1, length=n.scenarios+1)[-1]+
+ seq(0, 1, length=n.scenarios+1)[-(n.scenarios+1)])/2
+l <- matrix(0, ncol(defaultprob), n.scenarios)
+r <- matrix(0, ncol(defaultprob), n.scenarios)
+
+MFdist <- MFlossdistrib(w.mod, rho, defaultprob, p, issuerweights, Ngrid)
+MFdist.orig <- MFlossdistrib(w, rho, defaultprob, defaultprob, issuerweights, Ngrid)
+lossdist.orig <- BClossdistC(SurvProb, issuerweights, recov, rho, 101, 500)
+for(i in 1:17){
+ Lfun <- splinefun(c(0, cumsum(MFdist$L[,i])),c(0, seq(0, 1, length=Ngrid)), "monoH.FC")
+ Rfun <- splinefun(c(0, cumsum(MFdist$R[,i])),c(0, seq(0, 1, length=Ngrid)), "monoH.FC")
+ for(j in 1:99){
+ l[i, j] <- Lfun(d1[j])
+ r[i, j] <- Rfun(d1[j])
+ }
+}
+
+## joint generation of scenarios for loss and recovery distribution
+clusterExport(cl, list("lossrecovdist.joint.term", "lossdistribC.joint"))
+MFdist2 <- MFlossdistrib2(cl, w.mod, rho, defaultprob, p, issuerweights, Ngrid)
+## memory never gets released by the clusters for some reasons, needs to call gc twice
+clusterCall(cl, gc)
+clusterCall(cl, gc)
+gc()
+#L+R<=1
+test <- apply(MFdist2[21,,], 2, rev)
+test.tri <- lower.tri(test)
+sum(test[test.tri])#close to 1
+
+## D=L+R
+dist <- MFdist2[21,,]
+supportD <- outer(seq(0,1, length=Ngrid), seq(0,1, length=Ngrid), "+")
+## compute the joint density of (L/D, D)
+## u=x/(x+y)
+## v=x+y
+x <- seq(0,1, length=Ngrid) %o% seq(0,1, length=Ngrid)
+y <- sweep(-x, 2, seq(0, 1, length=Ngrid), "+")
+xgrid <- round(x/0.005)
+ygrid <- round(y/0.005)
+distchv <- matrix(0, Ngrid, Ngrid)
+for(i in 1:Ngrid){
+ for(j in 1:Ngrid){
+ distchv[i,j] <- dist[xgrid[i,j]+1, ygrid[i,j]+1]*seq(0,1, length=Ngrid)[j]
+ }
+}
+
+## compute reinvesting distribution
+beta <- 1.1
+f <- function(l, r)(1-l-r)/(1-r)^beta
+test <- outer(seq(0, 1, length=Ngrid), seq(0, 1, length=Ngrid), f)
+plot(density(x=test[-Ngrid*Ngrid], weights=dist[-Ngrid*Ngrid], from=0, to=5))
+
+## two dimensional scenarios
+n.scenarios <- 100
+percentiles <- (seq(0, 1, length=n.scenarios+1)[-1]+
+ seq(0, 1, length=n.scenarios+1)[-(n.scenarios+1)])/2
+l <- matrix(0, ncol(defaultprob), n.scenarios)
+r <- matrix(0, ncol(defaultprob), n.scenarios)
+
+MFdist <- MFlossdistrib(w.mod, rho, defaultprob, p, issuerweights, Ngrid)
+i <- 21
+Lfun <- splinefun(c(0, cumsum(MFdist$L[,i])),c(0, seq(0, 1, length=Ngrid)), "monoH.FC")
+
+Rfun <- splinefun(c(0, cumsum(MFdist$R[,i])),c(0, seq(0, 1, length=Ngrid)), "monoH.FC")
+ for(j in 1:99){
+ l[i, j] <- Lfun(d1[j])
+ r[i, j] <- Rfun(d1[j])
+ }
+}
diff --git a/R/cds_functions.R b/R/cds_functions.R new file mode 100644 index 00000000..7404328d --- /dev/null +++ b/R/cds_functions.R @@ -0,0 +1,343 @@ +source("cds_utils.R")
+
+setClass("defaultcurve", representation(dates="Date", hazardrates="numeric"))
+setClass("defaultprepaycurve", representation(prepayrates="numeric"), contains="defaultcurve")
+setClass("creditcurve", representation(issuer="character", startdate="Date",
+ recovery="numeric", curve="defaultprepaycurve"))
+
+survivalProbability1 <- function(startdate, date, survival.curve) {
+ #based on a flat hazard rate curve
+ T <- yearFrac(startdate, survival.curve$dates)
+ Tmat <- yearFrac(startdate, date)
+ for ( i in 1:length(survival.curve$dates) ){
+ if ( date >= survival.curve$dates[i] ) {
+ next
+ }else{
+ if( i > 1 ) {
+ w <- ( Tmat - T[i-1] ) / (T[i] - T[i-1])
+ logprob <- - (1-w) * survival.curve$hazardrates[i-1] * T[i-1] -
+ w * survival.curve$hazardrates[i] * T[i]
+ }else{
+ logprob <- - Tmat * survival.curve$hazardrates[1]
+ return( exp(as.numeric(logprob)) )
+ }
+ }
+ }
+ ## if date is greater than last survival.curve date, keep the hazard rate flate
+ logprob <- - yearFrac(startdate, date) * survival.curve$hazardrates[i]
+ return( exp(as.numeric(logprob)) )
+}
+
+survivalProbability.exact <- function(credit.curve, date) {
+ #based on a forward hazard rate curve
+ curve <- credit.curve@curve
+ T <- c(0, yearFrac(credit.curve@startdate, curve@dates))
+ dT <- diff(T)
+ Tmat <- yearFrac(credit.curve@startdate, date)
+ logprob <- 0
+ for ( i in 1:length(dT) ){
+ if ( date > curve@dates[i] ) {
+ logprob <- logprob - curve@hazardrates[i] * dT[i]
+ }else{
+ if( i > 1 ){
+ logprob <- logprob - curve@hazardrates[i] * (Tmat - T[i])
+ }else{
+ logprob <- logprob - curve@hazardrates[1] * Tmat
+ }
+ break
+ }
+ }
+ return( exp(as.numeric(logprob)) )
+}
+
+PD <- function(sc){
+ ## computes the default probability associated with the survival curve
+ if(class(sc)=="defaultprepaycurve"){
+ T <- c(0, yearFrac(today(), sc@dates))
+ dT <- diff(T)
+ return( 1-cumprod(exp(-sc@hazardrates * dT)) )
+ }else{
+ stop("not a default curve")
+ }
+}
+
+couponleg <- function(cs, sc, accruedondefault=TRUE){
+ ## computes the pv of the risky coupon leg based on a given coupon schedule
+ ## and a survival curve. Also called premium leg or fixed leg.
+ dT <- diff(c(0, yearFrac(today(), cs$dates)))
+ PD <- cumprod(exp(-sc@hazardrates[1:length(dT)] * dT))
+ if(accruedondefault){
+ PDadj <- 0.5 * (c(1, PD[-length(PD)]) + PD)
+ }else{
+ PDadj <- PD
+ }
+ return( crossprod(PDadj, cs$coupons * cs$df) )
+}
+
+couponleg.flat <- function(cs, h, accruedondefault=TRUE){
+ T <- yearFrac(today(), cs$dates)
+ PD <- exp(- h * T )
+ if(accruedondefault){
+ PDadj <- 0.5 * (c(1, PD[-length(PD)]) + PD)
+ }else{
+ PDadj <- PD
+ }
+ return( crossprod(PDadj, cs$coupons * cs$df) )
+}
+
+couponleg.prepay <- function(cs, dpc, accruedondefault=TRUE){
+ ## computes the pv of the risky coupon leg from the coupon schedule,
+ ## a hazard rate curve, and a prepay curve. We assume the poisson process driving
+ ## default and prepayment are independent, so the intensity of either event
+ ## happenning is the sum of the two.
+
+ dT <- diff(c(0, yearFrac(today(), cs$dates)))
+ SP <- cumprod(exp( - (dpc@hazardrates[1:length(dT)] + dpc@prepayrates[1:length(dT)] ) * dT))
+ if(accruedondefault){
+ SPadj <- 0.5 * (c(1, SP[-length(SP)]) + SP)
+ }else{
+ SPadj <- SP
+ }
+ return( crossprod(SPadj, cs$coupons * cs$df) )
+}
+
+couponleg.prepay.flat <- function(cs, h){
+ dpc <- new("defaultprepaycurve", dates = cs$dates, hazardrates=rep(h, length(cs$dates)),
+ prepayrates=rep(k(h), length(cs$dates)))
+ couponleg.prepay(cs, dpc)
+}
+
+dcouponleg <- function(cs, sc, index, accruedondefault=TRUE){
+ ## derivative of couponleg with respect to hazardrate
+ dT <- diff(c(0, yearFrac(today(), cs$dates)))
+ PD <- cumprod(exp(-sc@hazardrates[1:length(dT)] * dT))
+ dPD <- cumsum(index * dT) * PD
+ if(accruedondefault){
+ dPDadj <- 0.5 * (c(0, dPD[-length(PD)]) + dPD)
+ }else{
+ dPDadj <- dPD
+ }
+ return( - crossprod( dPDadj, cs$coupons * cs$df) )
+}
+
+k <- function(h, alpha=0.25, beta=15){
+ ## prepay rate as a function of the hazardrate
+ ## this is a decreasing function
+ ## alpha is the maximum prepay rate
+ return ( alpha * exp(- beta * h) )
+}
+
+dcouponleg.prepay <- function(cs, dpc, index, beta, accruedondefault=TRUE){
+ ## derivative of couponleg.prepay with respect to hazardrate
+ ## If k is the prepay rate, it assumes dk/dh = - beta * k,
+ ## which is the case if k(h) = alpha * exp(-beta *h)
+ dT <- diff(c(0, yearFrac(today(), cs$dates)))
+ SP <- cumprod(exp( - (dpc@hazardrates[1:length(dT)] + dpc@prepayrates[1:length(dT)] ) * dT))
+ dSP <- -cumsum(index * dT * (1 - beta * dpc@prepayrates[1:length(dT)])) * SP
+ if(accruedondefault){
+ dSPadj <- 0.5 * (c(0, dSP[-length(SP)]) + dSP)
+ }else{
+ SPadj <- dSP
+ }
+ return( crossprod(dSPadj, cs$coupons * cs$df) )
+}
+
+dcouponleg.prepay.flat <- function(cs, h, index, beta, accruedondefault=TRUE){
+ hvec <- rep(h, length(cs$dates))
+ kvec <- k(h)
+ dpc <- new("defaultprepaycurve", dates=cs$dates, hazardrates=hvec
+ prepayrates=kvec )
+ return( dcouponleg.prepay(cs, dpc, index, beta, accruedondefault))
+}
+
+cdsduration <- function(sc, maturity, accruedondefault=TRUE){
+ # computes the risky PV01, also called risky annuity of a cds
+ cs <- couponSchedule(nextIMMDate(today()), maturity, "Q", "FIXED", 1)
+ couponleg(cs, sc, accruedondefault)
+}
+
+defaultleg.flat <- function(cs, h, recovery){
+ ## Computes the pv of the default leg of a cds based on a given
+ ## coupon schedule, flat hazard rate, and recovery.
+ T <- yearFrac(today(), cs$dates)
+ Q <- exp(-h * T) * cs$df
+ Qmid <- 1/2 * (c(1, Q[-length(Q)]) + Q)
+ r <- (1 - recovery) * crossprod(h * Qmid, dT)
+ return( r )
+}
+
+ddefaultleg.flat <- function(cs, h, recovery){
+ ## derivative of defaultleg.flat with respect to hazardrate
+ T <- yearFrac(today(), cs$dates)
+ dT <- diff(c(0, T))
+ dQ <- - T * exp(-h * T) * cs$df
+ Q <- exp(-h * T) * cs$df
+ Qmid <- 1/2 * (c(1, Q[-length(Q)]) + Q)
+ dQmid <- 1/2 * (c(0, dQ[-length(dQ)]) + dQ)
+ dr <- (1-recovery) * (crossprod(Qmid, dT) + h * crossprod(dQmid, dT))
+ return( dr )
+}
+
+defaultleg <- function(cs, sc, recovery){
+ ## Computes the pv of the default leg of a cds based on a given
+ ## coupon schedule, hazard rate curve, and recovery.
+ T <- yearFrac(today(), cs$dates)
+ dT <- diff(c(0, T))
+ Q <- cumprod(exp(-sc@hazardrates[1:length(dT)] * dT)) * cs$df
+ Qmid <- 1/2 * (c(1, Q[-length(Q)]) + Q)
+ r <- (1 - recovery) * crossprod(sc@hazardrates[1:length(dT)] * Qmid, dT)
+ return( r )
+}
+
+ddefaultleg <- function(cs, sc, recovery, index){
+ ## derivative of defaultleg with respect to hazardrate
+ T <- yearFrac(today(), cs$dates)
+ dT <- diff(c(0,T))
+ Q <- cumprod(exp(-sc@hazardrates[1:length(dT)] * dT)) * cs$df
+ dQ <- - cumsum(index * dT) * Q
+ Qmid <- 1/2 * (c(1, Q[-length(Q)]) + Q)
+ dQmid <- 1/2 *(c(0, dQ[-length(dQ)]) + dQ)
+ dr <- (1-recovery) * (crossprod(index * Qmid, dT) + crossprod(sc@hazardrates * dQmid, dT))
+ return( dr )
+}
+
+defaultleg.prepay <- function(cs, sc, recovery){
+ ## Computes the pv of the default leg of a cds based on a given
+ ## coupon schedule, hazard rates curve, prepay curves, and recovery.
+ T <- yearFrac(today(), cs$dates)
+ dT <- diff(c(0, T))
+ Q <- cumprod(exp(- (dpc@hazardrates[1:length(dT)]+dpc@prepayrates[1:length(dT)]) * dT)) * cs$df
+ Qmid <- 1/2 * (c(1, Q[-length(Q)]) + Q)
+ r <- (1 - recovery) * crossprod(dpc@hazardrates[1:length(dT)] * Qmid, dT)
+ +crossprod(dpc@prepayrates[1:length(dT)] * Qmid, dT)
+ return( r )
+}
+
+cdspv <- function(cs, sc, R){
+ return ( couponleg(cs, sc) - defaultleg(cs, sc, R) )
+}
+
+cdsspread <- function(sc, maturity, recovery){
+ ## computes exact cds running spread for a cds
+ ## should be very close to hazardrate * (1-recovery)
+ cs <- couponSchedule(nextIMMDate(today()), maturity, "Q", "FIXED", 1)
+ defaultleg(cs, sc, recovery)/couponleg(cs, sc)
+}
+
+dcdspv <- function(cs, sc, recovery, index){
+ ## derivative of cdspv with respect to hazardrate
+ return ( dcouponleg(cs, sc, index) - ddefaultleg(cs, sc, recovery, index) )
+}
+
+cdshazardrate.flat <- function(upfront, running, maturity, R=0.4){
+ ## computes the implied hazard rate of the cds based on the upfront
+ ## and running quotes, as well as maturity and recovery
+ cs <- couponSchedule(nextIMMDate(today()), maturity, "Q", "FIXED", running)
+ lambda <- 0.05
+ eps <- 1e-8
+ while(abs(cdspv(lambda, cs, R) + upfront) > eps){
+ lambda <- lambda - (upfront + cdspv(lambda, cs, R))/dcdspv(lambda, cs, R)
+ }
+ return(lambda)
+}
+
+cdshazardrate <- function(quotes, R=0.4){
+ ## bootstrap the implied hazard rate curve of the cds based on the upfront
+ ## and running quotes, as well as maturity and recovery
+ previous.maturity <- today()
+ hvec <- c()
+ previous.hvec <- c()
+ eps <- 1e-8
+ previous.cs <- data.frame()
+ for(i in 1:nrow(quotes)){
+ maturity <- quotes$maturity[i]
+ cs <- couponSchedule(nextIMMDate(today()), maturity, "Q", "FIXED", quotes$running[i])
+ new.h <- 0.05
+ flength <- nrow(cs) - nrow(previous.cs)
+ hvec <- c(previous.hvec, rep(new.h, flength))
+ sc <- new("defaultprepaycurve", dates=cs$dates, hazardrates=hvec)
+ index <- c(rep(0, length(previous.hvec)), rep(1, flength))
+
+ while(abs(cdspv(cs, sc, R) + quotes$upfront[i]) > eps){
+ new.h <- new.h - (quotes$upfront[i] + cdspv(cs, sc, R))/dcdspv(cs, sc, R, index)
+ hvec <- c(previous.hvec, rep(new.h, flength))
+ sc@hazardrates <- hvec
+ }
+ previous.hvec <- hvec
+ previous.maturity <- maturity
+ previous.cs <- cs
+ }
+ return(sc)
+}
+
+
+bonddp <- function(collateral, R=0.7){
+ ## bootstrap the implied hazard rate curve of a bond based on the upfront
+ ## and running quotes, as well as maturity and recovery
+ previous.maturity <- today()
+ hvec <- c()
+ previous.hvec <- c()
+ eps <- 1e-8
+ previous.cs <- data.frame()
+ for(i in 1:nrow(quotes)){
+ maturity <- quotes$maturity[i]
+ cs <- couponSchedule(collateral$nextpaydate, collateral$maturity,
+ collateral$frequency, collateral$fixedorfloat,
+ collateral$grosscoupon*0.01, collateral$spread*0.01)
+ new.h <- 0.05
+ flength <- nrow(cs)-nrow(previous.cs)
+ hvec <- c(previous.hvec, rep(new.h, flength))
+ sc <- new("survivalcurve", dates=cs$dates, hazardrates=hvec)
+ index <- c(rep(0, length(previous.hvec)), rep(1, flength))
+
+ while(abs(cdspv(cs, sc, R) + quotes$upfront[i]) > eps){
+ new.h <- new.h - (quotes$upfront[i] + cdspv(cs, sc, R))/dcdspv(cs, sc, R, index)
+ hvec <- c(previous.hvec, rep(new.h, flength))
+ sc@hazardrates <- hvec
+ }
+ previous.hvec <- hvec
+ previous.maturity <- maturity
+ previous.cs <- cs
+ }
+ return(sc)
+}
+
+indexpv <- function(portfolio, index, epsilon=0){
+ ## computes the intrinsic index pv of a portfolio of cds
+ r <- rep(0, length(portfolio))
+ cs <- couponSchedule(nextIMMDate(today()), index$maturity, "Q", "FIXED", index$coupon)
+ for(i in 1:length(portfolio)){
+ if(epsilon!=0){
+ tweakedcurve <- portfolio[[i]]@curve
+ tweakedcurve@hazardrates <- tweakedcurve@hazardrates * (1 + epsilon)
+ r[i] <- cdspv(cs, tweakedcurve, index$recovery)
+ }else{
+ r[i] <- cdspv(cs, portfolio[[i]]@curve, index$recovery)
+ }
+ }
+ return( 1+mean(r) )
+}
+
+tweakcurves <- function(portfolio, index){
+ ## computes the tweaking factor
+ epsilon <- 0
+ f <- function(epsilon, ...){
+ abs(indexpv(portfolio, index, epsilon)-index$indexref)
+ }
+ epsilon <- optimize(f, c(-0.5, 0.5), portfolio, index, tol=1e-6)$minimum
+ portfolio.new <- portfolio
+ for(i in 1:length(portfolio)){
+ portfolio.new[[i]]@curve@hazardrates <- portfolio[[i]]@curve@hazardrates * (1 + epsilon)
+ }
+ return( portfolio.new )
+}
+
+SPmatrix <- function(portfolio, index){
+ cs <- couponSchedule(nextIMMDate(today()), index$maturity, "Q", "FIXED", index$coupon)
+ SP <- matrix(0, length(portfolio), length(cs$dates))
+ for(i in 1:length(portfolio)){
+ SP[i,] <- 1 - PD(portfolio[[i]]@curve)[1:length(cs$dates)]
+ }
+ return( SP )
+}
+
diff --git a/R/cds_functions_generic.R b/R/cds_functions_generic.R new file mode 100644 index 00000000..e6e5a654 --- /dev/null +++ b/R/cds_functions_generic.R @@ -0,0 +1,664 @@ +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))
+}
diff --git a/R/cds_utils.R b/R/cds_utils.R new file mode 100644 index 00000000..afb8f400 --- /dev/null +++ b/R/cds_utils.R @@ -0,0 +1,146 @@ +library("RQuantLib")
+
+today <- function() {
+ Sys.Date()
+}
+
+convertTenor <- function(tenor) {
+ ## convert tenors of the form '1y', '2y', etc...
+ ## and '1m', '2m'... into yearfrac
+ month <- regexpr("([0-9]+)m", tenor, perl = T)
+ year <- regexpr("([0-9]+)y", tenor, perl = T)
+ if ( month != -1 ) {
+ a <- attr(month, "capture.start")
+ b <- a + attr(month, "capture.length") - 1
+ l <- as.numeric(substr(tenor, a, b))
+ return ( 30 * l )
+ }else if ( year != -1) {
+ a <- as.numeric(attr(year, "capture.start"))
+ b <- a + attr(year, "capture.length") - 1
+ l <- as.numeric(substr(tenor, a,b))
+ return ( 365 * l )
+ }else{
+ stop("format not recognized")
+ }
+}
+
+addTenor <- function(date, tenor) {
+ month <- regexpr("([0-9]+)m", tenor, perl = T)
+ year <- regexpr("([0-9]+)y", tenor, perl = T)
+ if ( month != -1 ) {
+ a <- attr(month, "capture.start")
+ b <- a + attr(month, "capture.length") - 1
+ l <- as.numeric(substr(tenor, a, b))
+ return ( seq(date, length=2, by=paste(l,"month"))[2])
+ }else if ( year != -1) {
+ a <- as.numeric(attr(year, "capture.start"))
+ b <- a + attr(year, "capture.length") - 1
+ l <- as.numeric(substr(tenor, a,b))
+ return ( seq(date, length=2, by=paste(l,"year"))[2] )
+ }else{
+ stop("format not recognized")
+ }
+}
+
+creditSchedule <- function(startdate, enddate) {
+ ## generate a credit coupon payment schedule from startdate
+ ## to enddate based on the IMM dates
+ month <- c(3, 6, 9, 12)
+ day <- 20
+ startmonth <- as.numeric(format(startdate, "%m"))
+ startday <- as.numeric(format(startdate, "%d"))
+ startyear <- as.numeric(format(startdate, "%Y"))
+ if (startday > day){
+ followingmonth <- which(month>startmonth)
+ if (length(followingmonth)==0){
+ startyear <- startyear+1
+ startmonth <- 3
+ }else{
+ startmonth <- followingmonth[1]
+ }
+ }else{
+ startmonth <- which(month>=startmonth)[1]
+ }
+ newdate <- as.Date(paste(startyear, month[startmonth], day, sep="-"))
+ seq(newdate, enddate, by="3 months")
+}
+
+couponSchedule <- function(nextpaydate=NULL, maturity, frequency, coupontype, currentcoupon,
+ margin, startdate=today()){
+ ## computes the coupon schedule
+ ## inputs:
+ ## nextpaydate: first payment date of the coupon schedule
+ ## maturity: last payment date of the schedule
+ ## frequency: letter specifying the frequency amon "Q", "M", "B", or "S"
+ ## if startdate is provided, this generates the forward coupon schedule starting from that date.
+ bystring <- switch(frequency,
+ Q = "3 months",
+ M = "1 month",
+ B = "2 months",
+ S = "6 months",
+ A = "12 months")
+ if(is.null(bystring)){
+ stop("unknown frequency")
+ }
+ if(is.null(nextpaydate)){
+ dates <- rev(seq(maturity, today(), by =paste0("-", bystring)))
+ }else{
+ dates <- seq(nextpaydate, maturity, by = bystring)
+ }
+ if(dates[length(dates)]<maturity){
+ dates <- c(dates, maturity)
+ }
+ dates <- dates[ dates >= today()]
+ DC <- switch(frequency,
+ S = DiscountCurve(L6m$params, L6m$tsQuotes, yearFrac(L6m$params$tradeDate, dates)),
+ Q = DiscountCurve(L3m$params, L3m$tsQuotes, yearFrac(L3m$params$tradeDate, dates)),
+ M = DiscountCurve(L1m$params, L1m$tsQuotes, yearFrac(L1m$params$tradeDate, dates)),
+ B = DiscountCurve(L2m$params, L2m$tsQuotes, yearFrac(L2m$params$tradeDate, dates)),
+ A = DiscountCurve(L12m$params, L12m$tsQuotes, yearFrac(L12m$params$tradeDate, dates)))
+
+ if(coupontype=="FLOAT" && !is.na(margin)){ #if is.na(margin) probably letter of credit
+ coupons <- pmax(currentcoupon, DC$forwards + margin)
+ }else{
+ coupons <- rep(currentcoupon, length(dates))
+ }
+ coupons <- diff(c(0, yearFrac(startdate, dates, "act/360"))) * coupons
+ if(startdate!=DC$params$tradeDate){
+ df <- cumprod(exp(-DC$forwards * diff(c(0, yearFrac(startdate, dates)))))
+ }else{
+ df <- DC$discounts
+ }
+ return( data.frame(dates=dates, coupons=coupons, df = df) )
+}
+
+nextIMMDate <- function(date) {
+ startyear <- as.numeric(format(date, "%Y"))
+ nextimmdates <- seq(as.Date(paste(startyear, 3, 20, sep="-")), length=5, by="3 months")
+ return( nextimmdates[nextimmdates >= date][1] )
+}
+
+cdsMaturity <- function(tenor, date=today()){
+ enddate <- addTenor(date, tenor)
+ month <- c(3, 6, 9, 12)
+ day <- 20
+ startmonth <- as.numeric(format(enddate, "%m"))
+ startday <- as.numeric(format(enddate, "%d"))
+ startyear <- as.numeric(format(enddate, "%Y"))
+ if (startday > day) {
+ followingmonth <- which(month>startmonth)
+ if (length(followingmonth)==0){
+ startyear <- startyear+1
+ startmonth <- 3
+ }else{
+ startmonth <- followingmonth[1]
+ }
+ }else{
+ startmonth <- which(month >= startmonth)[1]
+ }
+ return( as.Date(paste(startyear, month[startmonth], day, sep="-")))
+}
+
+yearFrac <- function(date1, date2, daycount="act/365") {
+ switch(daycount,
+ "act/365"=as.numeric( (as.Date(date2) - as.Date(date1)) / 365),
+ "act/360"=as.numeric( (as.Date(date2) - as.Date(date1)) / 360) )
+}
diff --git a/R/clopricer.R b/R/clopricer.R new file mode 100644 index 00000000..d428c2e8 --- /dev/null +++ b/R/clopricer.R @@ -0,0 +1,87 @@ +rootdir <- "T:/Analytics/CLO/"
+setwd(rootdir)
+options(stringsAsFactors=F)
+
+stripext <- function(name,ext){
+ sub(paste(".",ext,sep=""),"",name)
+}
+
+today <- function(){
+ as.Date(Sys.time())
+}
+
+listcusips <- function(date, flag=""){
+ clodir <- paste(format(date, "%Y%m%d"), flag, sep="")
+ cusips <- dir(path=paste(".",clodir,sep="/"), pattern=".report$")
+ return(stripext(cusips,"report"))
+}
+
+readcusip <- function(cusip,date,flag=""){
+ clodir <- paste(format(date,"%Y%m%d"),flag,sep="")
+ temp <- read.table(paste(".",clodir,paste(cusip,"report",sep="."),sep="/"),sep="\t",header=T,fill=T)
+ temp <- temp[temp$Scenario!="Scenario",1:25]
+ a <-as.matrix(temp$Scenario)
+ a <- sub("SCENARIO","",as.character(a))
+ temp <- data.frame(a,temp[,2:25])
+ tempnames <- colnames(temp)
+ colnames(temp) <- c("Scenario",tempnames[2:25])
+ return(temp[order(as.numeric(temp$Scenario)),])
+}
+
+readcusip2 <- function(cusip,date){
+ clodir <- format(date,"%Y%m%d")
+ temp <- read.table(paste(".",clodir,paste(cusip,"report",sep="."),sep="/"),sep="\t",header=T,fill=T)
+ temp <- temp[temp$Scenario!="Scenario",1:25]
+ a <- t(matrix(unlist(strsplit(as.character(temp$Scenario),"_")),2,nrow(temp)))
+ a[,1] <- sub("SCENARIO","",as.character(a[,1]))
+ colnames(a) <- c("Scenario","ScenarioType")
+ temp <- data.frame(a,temp[,2:25])
+ return(temp[order(as.numeric(temp$Scenario)),])
+}
+getcusip <- function(cusip,date,colname,scenariotype){
+ cusipdata <- readcusip2(cusip,date)
+ as.numeric(cusipdata[cusipdata$ScenarioType==scenariotype,][[colname]])
+}
+
+getcusipinfo <- function(cusip,date,scenariotype,infotype,indextype,indextenor){
+ clodir <- paste(format(date,"%Y%m%d"),"Info",sep="/")
+ cusipscen <- paste(cusip,scenariotype,sep="_")
+ temp <- read.table(paste(".",clodir,paste(cusipscen,infotype,sep="."),sep="/"),sep="\t",header=T,fill=T)
+ return(temp[(temp$IndexType==indextype)&(temp$IndexTenor==indextenor),])
+}
+
+A <- c()
+for(cusip in cusipsBB){
+ A <- cbind(A,as.numeric(readcusip(cusip,today())$Price))
+}
+colnames(A) <- cusipsBB
+
+pricesAAA <- c()
+for(cusip in cusipsAAA){
+ pricesAAA <- rbind(pricesAAA,readcusip(cusip,today())[,c("Price","Duration","Yield","WAL")])
+}
+
+pricesAA <- c()
+for(cusip in cusipsAA){
+ pricesAA <- rbind(pricesAA,readcusip(cusip,today())[,c("Price","Duration","Yield","WAL")])
+}
+
+pricesA <- c()
+for(cusip in cusipsA){
+ pricesA <- rbind(pricesA,readcusip(cusip,today())[,c("Price","Duration","Yield","WAL")])
+}
+
+pricesBBB <- c()
+for(cusip in cusipsBBB){
+ pricesBBB <- rbind(pricesBBB,readcusip(cusip,today())[,c("Price","Duration","Yield","WAL")])
+}
+
+pricesBB <- c()
+for(cusip in cusipsBB){
+ pricesBB <- rbind(pricesBB,readcusip(cusip,today())[,c("Price","Duration","Yield","WAL")])
+}
+
+pricesNA <- c()
+for(cusip in cusipsNA){
+ pricesNA <- rbind(pricesNA,readcusip(cusip,today())[,c("Price","Duration","Yield","WAL")])
+}
diff --git a/R/deal_pricer.R b/R/deal_pricer.R new file mode 100644 index 00000000..0fb821de --- /dev/null +++ b/R/deal_pricer.R @@ -0,0 +1,22 @@ +library(RPostgreSQL) +drv <- dbDriver("PostgreSQL") +dbCon <- dbConnect(drv, dbname="ET", user="et_user", password="Serenitas;1") +dealnames <- dbGetQuery(dbCon, "select dealname from clo_universe") +pricingcoverage <- c() +for (dealname in dealnames$dealname){ + r <- dbGetQuery(dbCon, paste( + paste("select sum(currentbalance* coalesce(b.bid,c.price))/sum(currentbalance) as wap,", + "sum(currentbalance) as pricedbalance from et_collateral a", + "left join latest_markit_prices b on a.loanxid=b.loanxid", + "left join bloomberg_corp c on a.cusip=c.cusip", + "where a.dealname='"), + dealname, + "' and coalesce(b.bid,c.price) is not Null",sep="")) + s <- dbGetQuery(dbCon, + paste("select sum(currentbalance) as totalbalance from et_collateral where dealname='",dealname,"'",sep="")) + pricingcoverage <- rbind(pricingcoverage, c(dealname,r$wap,r$pricedbalance/s$totalbalance)) +} +pricingcoverage <- data.frame(dealname=dealnames$dealname,wap=as.numeric(pricingcoverage[,2]), coverage = as.numeric(pricingcoverage[,3])) + +dbDisconnect(dbCon) + diff --git a/R/dumpdata.R b/R/dumpdata.R new file mode 100644 index 00000000..d0f4bc73 --- /dev/null +++ b/R/dumpdata.R @@ -0,0 +1,4 @@ +y <-
+c(1.28682479818149, 0.99111388484709, 0.617435861498514, -0.471105689808945,
+-1.18148209395743, 0.81499608398862, -0.923401317102897, 0.481822756404342,
+-1.09823529971577, 1.07505266897986)
diff --git a/R/ema.R b/R/ema.R new file mode 100644 index 00000000..b4af957e --- /dev/null +++ b/R/ema.R @@ -0,0 +1,3 @@ + ema<-function(x, beta = 0.1, init=x[1]){
+filter(beta*x,filter = 1-beta,method = "recursive", init= init)
+}
diff --git a/R/etdb.R b/R/etdb.R new file mode 100644 index 00000000..93b72e27 --- /dev/null +++ b/R/etdb.R @@ -0,0 +1,3 @@ +library(RPostgreSQL)
+drv <- dbDriver("PostgreSQL")
+dbCon <- dbConnect(drv, dbname="ET", user="et_user", password="Serenitas1")
diff --git a/R/filterruns.R b/R/filterruns.R new file mode 100644 index 00000000..b6d93f73 --- /dev/null +++ b/R/filterruns.R @@ -0,0 +1,25 @@ +lcds <- c()
+blcds <- c()
+ #non upfront case
+upfront <- grep("[0-9]+ *- *[0-9]+\\+500",x)
+quotes <- x[upfront]
+for(i in 1:length(quotes)){
+ matches <- gregexpr("[0-9]+ *- *[0-9]+\\+500",quotes[i])
+ row <- substring(quotes[i],matches[[1]],matches[[1]]+attr(matches[[1]],"match.length")-1)
+ row <- gsub("\\+500","",row)
+ row <- gsub(" ","",row)
+ lcds <- rbind(lcds,c(as.numeric(strsplit(row,"-")[[1]]),500,500))
+ blcds <- rbind(blcds,c(as.numeric(strsplit(row,"-")[[2]]),500,500))
+}
+quotes <- x[-upfront]
+running <- grep("[0-9]+ *- *[0-9]+",quotes)
+quotes <- quotes[running]
+for(i in 1:length(quotes)){
+ matches <- gregexpr("[0-9]+ *- *[0-9]+",quotes[i])
+ row <- substring(quotes[i],matches[[1]],matches[[1]]+attr(matches[[1]],"match.length")-1)
+ row <- gsub("\\+500","",row)
+ row <- gsub(" ","",row)
+ lcds <- rbind(lcds,c(0,0,as.numeric(strsplit(row,"-")[[1]])))
+ blcds <- rbind(blcds,c(0,0,as.numeric(strsplit(row,"-")[[2]])))
+}
+
diff --git a/R/index_definitions.R b/R/index_definitions.R new file mode 100644 index 00000000..71c61688 --- /dev/null +++ b/R/index_definitions.R @@ -0,0 +1,24 @@ +hy19 <- list(maturity=as.Date("2017-12-20"),
+ coupon=0.05,
+ factor=1,
+ loss=0,
+ recovery=0.3)
+
+hy18 <- list(maturity=as.Date("2017-06-20"),
+ coupon=0.05,
+ factor=0.99,
+ loss=0.0082375,
+ recovery=0.3)
+
+hy17 <- list(maturity=as.Date("2016-12-20"),
+ coupon=0.05,
+ factor=0.96,
+ loss=0.027075,
+ recovery=0.3,
+ defaulted=c("DYNHLD", "EK", "GM-ResCLLC", "PMI"))
+
+hy15 <- list(maturity=as.Date("2015-12-20"),
+ coupon=0.05,
+ factor=0.96,
+ loss=0.026375,
+ recovery=0.3)
diff --git a/R/interpweights.R b/R/interpweights.R new file mode 100644 index 00000000..03de58d7 --- /dev/null +++ b/R/interpweights.R @@ -0,0 +1,94 @@ +interpweights <- function(w, v1, v2){
+ #Given L=(w,v1), compute neww such that newL=(new,v2)=L in distribution
+ cumw <- cumsum(w)
+ neww <- splinefun(v1,cumw,method= "monoH.FC")(v2,deriv=1)
+ #neww <- diff(newcumw)
+ interpweights <- neww/sum(neww)
+ return(interpweights)
+}
+
+adjust_scenario <- function(scenario, epsilon){
+ 1-(1-scenario)^(1/(1+epsilon))
+}
+
+adjust_weights <- function(weights, scenario, epsilon){
+ interpweights(weights,scenario,adjust_scenario(scenario,epsilon))
+}
+
+obj <- function(epsilon, vecpv, prob, support, cte){
+ newprob <- adjust_weights(prob, support, epsilon)
+ return( 1 - crossprod(newprob, vecpv) - cte)
+}
+
+optimize <- function(min, max, vecpv, prob, support, cte){
+ mid <- (min + max)/2
+ objective <- obj(mid, vecpv, prob, support, cte)
+ while( abs(objective)>1e-6){
+ if(objective>0){
+ min <- mid
+ }else{
+ max <- mid
+ }
+ mid <- (min+max)/2
+ objective <- obj(mid, vecpv, prob, support, cte)
+ }
+ return( mid )
+}
+
+interpweightsadjust <- function(w, v1, v2, vecpv){
+ interpweightsadjust <- interpweights(w, v1, v2)
+ epsilon <- optimize(-0.5, 0.5, vecpv, interpweightsadjust, v2, 1)
+ return( adjust_weights(interpweightsadjust, v2, epsilon) )
+}
+
+transformweightslike <- function(p1, v1, p2, v2, p, v){
+ cump2 <- cumsum(p2)
+ cump1 <- cumsum(p1)
+ P1 <- splinefun(v1,cump1,method= "monoH.FC")
+ dP1 <- function(x){P1(x,deriv=1)}
+ pomme <- interpweights(p2,v2,v)
+ pomme <- cumsum(pomme)
+ r <- rep(0,length(pomme))
+ for(i in 1:length(pomme)){
+ r[i] <- inverse(P1,dP1,pomme[i])
+ }
+ return(r)
+}
+
+clipw <- function(x){
+ write(x,file="clipboard",sep="\n")
+}
+
+clipr <- function(){
+ scan(file="clipboard")
+}
+
+sclipr <- function(){
+ scan(file="clipboard",what="character")
+}
+
+inverse <- function(f,Df,x, x0=x){
+ #inverse a function by the newton's method.
+ x1 <- x0-f(x0)/(Df(x0)-1)
+ counter <- 0
+ while(abs(x1-x0)>1e-6&&counter<500){
+ x0 <- x1
+ x1 <- x0-(f(x0)-x)/Df(x0)
+ counter <- counter+1
+ }
+ return(x1)
+}
+setwd("W:/corpCDOS/NextGen")
+lcdx10 <- read.table("lcdx10",sep=",",header=T)
+lcdx12 <- read.table("lcdx12",sep=",",header=T)
+hy10_5y <- read.table("hy10_5y",sep=",",header=T)
+hy10_7y <- read.table("hy10_7y",sep=",",header=T)
+
+lcdx10_pv <- read.table("lcdx10_pv",sep=",",header=T)
+lcdx12_pv <- read.table("lcdx12_pv",sep=",",header=T)
+hy10_5y_pv <- read.table("hy10_5y_pv",sep=",",header=T)
+hy10_7y_pv <- read.table("hy10_7y_pv",sep=",",header=T)
+clipw(interpweightsadjust(lcdx10$dist,lcdx10$support,lcdx10_pv$support,lcdx10_pv$pv))
+clipw(interpweightsadjust(hy10_5y$dist,hy10_5y$support,hy10_5y_pv$support,hy10_5y_pv$pv))
+clipw(interpweightsadjust(hy10_7y$dist,hy10_7y$support,hy10_7y_pv$support,hy10_7y_pv$pv))
+clipw(interpweightsadjust(lcdx12$dist,lcdx12$support,lcdx12_pv$support,lcdx12_pv$pv))
diff --git a/R/interpweights_2.R b/R/interpweights_2.R new file mode 100644 index 00000000..189915f5 --- /dev/null +++ b/R/interpweights_2.R @@ -0,0 +1,56 @@ +interpweights <- function(w, v1, v2){
+ #Given L=(w,v1), compute neww such that newL=(new,v2)=L in distribution
+ cumw <- cumsum(w)
+ neww <- splinefun(v1,cumw,method= "monoH.FC")(v2,deriv=1)
+ #neww <- diff(newcumw)
+ interpweights <- neww/sum(neww)
+ return(interpweights)
+}
+
+adjust_scenario <- function(scenario, epsilon){
+ 1-(1-scenario)^(1/(1+epsilon))
+}
+
+adjust_weights <- function(weights, scenario, epsilon){
+ interpweights(weights,scenario,adjust_scenario(scenario,epsilon))
+}
+
+obj <- function(epsilon, vecpv, prob,support, cte){
+ newprob <- adjust_weights(prob, support, epsilon)
+ return( 1 - crossprod(newprob, vecpv) - cte)
+}
+
+optimize <- function(min, max, vecpv, prob, support, cte){
+ mid <- (min + max)/2
+ objective <- obj(mid, vecpv, prob, support, cte)
+ while( abs(objective)>1e-6){
+ if(objective>0){
+ min <- mid
+ }else{
+ max <- mid
+ }
+ mid <- (min+max)/2
+ objective <- obj(mid, vecpv, prob, support, cte)
+ }
+ return( mid )
+}
+
+interpweightsadjust <- function(w, v1, v2, vecpv){
+ interpweightsadjust <- interpweights(w, v1, v2)
+ epsilon <- optimize(-0.5, 0.5, vecpv, interpweightsadjust, v2, 1)
+ return( adjust_weights(interpweightsadjust, v2, epsilon) )
+}
+
+transformweightslike <- function(p1, v1, p2, v2, p, v){
+ cump2 <- cumsum(p2)
+ cump1 <- cumsum(p1)
+ P1 <- splinefun(v1,cump1,method= "monoH.FC")
+ dP1 <- function(x){P1(x,deriv=1)}
+ pomme <- interpweights(p2,v2,v)
+ pomme <- cumsum(pomme)
+ r <- rep(0,length(pomme))
+ for(i in 1:length(pomme)){
+ r[i] <- inverse(P1,dP1,pomme[i])
+ }
+ return(r)
+}
diff --git a/R/intex_deals_functions.R b/R/intex_deals_functions.R new file mode 100644 index 00000000..a85b9f24 --- /dev/null +++ b/R/intex_deals_functions.R @@ -0,0 +1,253 @@ +library(RQuantLib)
+library(statmod)
+root = "//WDSENTINEL/share/CorpCDOs/"
+source(file.path(root, "R", "yieldCurve.R"))
+source(file.path(root, "R", "cds_functions_generic.R"))
+source(file.path(root, "R", "etdb.R"))
+source(file.path(root, "R", "tranche_functions.R"))
+load(file.path(root, "R", "bloomberg_data.RData"))
+
+cdorating <- function(cusip){
+ return( sub("[0-9]","", dataMtge[dataMtge$CUSIP %in% cusip,]$RTG_MDY_INITIAL ))
+}
+
+getcollateral <- function(dealname, date=Sys.Date()){
+ sqlstring <- sprintf("select * from et_aggdealinfo_historical('%s', '%s')", dealname, date)
+ collatdata <- dbGetQuery(dbCon, sqlstring)
+ return(collatdata)
+}
+
+getdealdata <- function(dealnames){
+ sqlstring <- sprintf("select * from latest_clo_universe where dealname in ('%s')",
+ paste(dealnames, collapse="','"))
+ return( dbGetQuery(dbCon, sqlstring) )
+}
+
+recovery <- function(collateral) {
+ ## return assumed recovery based on assumptions from recovery.assumptions
+ if(!is.na(collateral$secondlien) && collateral$secondlien){
+ collateral$assettype <- "SecondLien"
+ }
+ recovery <- with(global.params,
+ as.numeric(recovery.assumptions[collateral$assettype]))
+ if( !is.na(collateral$covlite) && collateral$covlite) {
+ recovery <- recovery - global.params$recovery.assumptions$Adj_Covlite
+ }
+ if( !is.na(collateral$iscdo) && collateral$iscdo ){
+ recovery <- 0
+ }
+ ## price is too low need to lower the assumed recovery
+ if(!is.na(collateral$price) && recovery > collateral$price/100 - 0.1){
+ recovery <- max(collateral$price/100-0.2, 0)
+ }
+ return(recovery)
+}
+
+dealnamefromcusip <- function(cusips){
+ ## wrapper around the sql procedure, not the fastest probably
+ r <- NULL
+ for(i in 1:length(cusips)){
+ sqlstr <- sprintf("select * from dealname_from_cusip('%s')", cusips[i])
+ r <- c(r, as.character(dbGetQuery(dbCon, sqlstr)))
+ }
+ return( r )
+}
+
+cusipsfromdealnames <- function(dealnames){
+ unlist(strsplit(getdealdata(dealnames)$"Deal Cusip List", ","))
+}
+
+fithazardrate.fast <- function(collateral, eps=1e-6){
+ lambda <- 0.05
+ cs <- couponSchedule(collateral$nextpaydate, collateral$maturity,
+ collateral$frequency, collateral$fixedorfloat,
+ collateral$grosscoupon * 0.01, collateral$spread*0.01)
+ R <- recovery(collateral)
+ while(abs(bondprice(lambda, cs, R) * 100 - collateral$price) > eps){
+ lambda <- lambda - (bondprice(lambda, cs, R) - 0.01*collateral$price)/dbondprice(lambda, cs, R)
+ }
+ return( lambda )
+}
+
+vanillabondprice <- function(h, collateral, prepay=TRUE) {
+ R <- recovery(collateral)
+ cs <- couponSchedule(collateral$nextpaydate, collateral$maturity,
+ collateral$frequency, collateral$fixedorfloat,
+ collateral$grosscoupon*0.01, collateral$spread*0.01)
+ if(prepay){
+ dpc <- new("defaultprepaycurve", dates=cs$dates, hazardrates=rep(h,length(cs$dates)),
+ prepayrates=rep(k(h), length(cs$dates)))
+ }else{
+ dpc <- new("defaultprepaycurve", dates=cs$dates, hazardrates=rep(h,length(cs$dates)),
+ prepayrates=numeric(0))
+ }
+ return( bondprice(cs, dpc, R) )
+}
+
+dvanillabondprice <- function(hazardrate, collateral) {
+ R <- recovery(collateral)
+ cs <- couponSchedule(collateral$nextpaydate, collateral$maturity,
+ collateral$frequency, collateral$fixedorfloat,
+ collateral$grosscoupon*0.01, collateral$spread*0.01)
+ return( bondprice(hazardrate, cs, R) )
+}
+
+fithazardrate <- function(collateral){
+ R <- recovery(collateral)
+ cs <- couponSchedule(collateral$nextpaydate, collateral$maturity,
+ collateral$frequency, collateral$fixedorfloat,
+ collateral$grosscoupon*0.01, collateral$spread*0.01)
+ f <- function(lambda){
+ u <- bondprice(lambda, cs, R )
+ return( (u * 100-collateral$price)^2 )
+ }
+ return( optimize(f, c(0,1), tol=1e-6)$minimum )
+}
+
+
+maturity <- function(creditcurve){
+ if(class(creditcurve)=="creditcurve"){
+ dates <- creditcurve@curve@dates
+ if(length(dates)){
+ return( dates[length(dates)] )
+ }else{
+ return( creditcurve@startdate )
+ }
+ }else{
+ stop("not of class creditcurve")
+ }
+}
+
+stackcurve <- function(SC, line.item, global.params, startdate){
+ newdates <- seq(startdate, line.item$maturity, by="3 months")
+ if(line.item$assettype=="Loan"){
+ hvec <- global.params$shape(yearFrac(today(), newdates[-1])) * global.params$defaultloanhazardrate
+ kvec <- global.params$alpha * exp(-global.params$beta * hvec)
+ }else if(line.item$assettype=="Bond" || line.item$iscdo){
+ hvec <- global.params$shape(yearFrac(today(), newdates[-1])) * global.params$defaultbondhazardrate
+ kvec <- rep(0, length(hvec))
+ }
+ SC@curve@hazardrates <- c(SC@curve@hazardrates, hvec)
+ SC@curve@prepayrates <- c(SC@curve@prepayrates, kvec)
+ SC@curve@dates <- c(SC@curve@dates, newdates[-1])
+ return(SC)
+}
+
+buildSC.matured <- function(SC, line.item, reinvdate, dealmaturity, global.params, startdate){
+ if(startdate<=reinvdate){ #reinvest
+ line.item$maturity <- min(dealmaturity, startdate + global.params$rollingmaturity)
+ SC <- stackcurve(SC, line.item, global.params, startdate)
+ }else{ #no reinvestment
+ SC@curve@dates <- startdate
+ SC@curve@hazardrates <- 0
+ SC@curveprepayrates <- 0
+ }
+ return( SC )
+}
+
+buildSC <- function(line.item, reinvdate, dealmaturity, global.params, startdate){
+ ## cat(line.item$issuername, "\n")
+ if(!is.na(line.item$iscdo) && line.item$iscdo && is.na(line.item$price)){
+ ##we have prices for some cdos e.g. 210795PS3
+ orig.moody <- cdorating(line.item$cusip)
+ if(length(orig.moody)==0){
+ orig.moody <- "NR"
+ }
+ line.item$price <- as.numeric(global.params$cdoprices[orig.moody])
+ }
+ ##build survival curve
+ SC <- new("creditcurve", recovery=recovery(line.item), startdate=startdate,
+ issuer=line.item$issuername)
+ SC@curve <- new("defaultprepaycurve", dates=as.Date(character(0)))
+ ## defaulted asset
+ if(!is.na(line.item$defaultedflag) && line.item$defaultedflag){
+ if(!is.na(line.item$price)){
+ line.item$currentbalance <- line.item$currentbalance * line.item$price/100
+ }else{
+ line.item$currentbalance <- line.item$currentbalance * recovery(line.item)
+ }
+ SC@startdate <- startdate + global.params$defaultedlag
+ line.item$maturity <- min(dealmaturity, SC@startdate + global.params$rollingmaturity)
+ ## automatic reinvest
+ SC<- stackcurve(SC, line.item, global.params, SC@startdate)
+ }else if(line.item$maturity<=startdate){#matured asset
+ SC <- buildSC.matured(SC, line.item, reinvdate, dealmaturity, global.params, startdate)
+ }else if(is.na(line.item$price)){ #missing price
+ SC <- stackcurve(SC, line.item, global.params, SC@startdate)
+ }else{ #normal case
+ if(line.item$assettype=="Bond"){ #no prepay rate
+ alpha <- 0
+ }else{
+ alpha <- global.params$alpha
+ }
+ try <- bondhazardrate.shaped(line.item, global.params$shape,
+ recovery(line.item), alpha)
+ if(is.null(try)){
+ SC <- stackcurve(SC, line.item, global.params, SC@startdate)
+ }else{
+ SC@curve <- try
+ }
+ }
+ if(maturity(SC) <= reinvdate){ #we reinvest
+ newstartdate <- line.item$maturity
+ line.item$maturity <- min(dealmaturity, newstartdate + global.params$rollingmaturity)
+ SC <- stackcurve(SC, line.item, global.params, newstartdate)
+ }
+ return( list(SC=SC, notional=line.item$currentbalance) )
+}
+
+buildSC.portfolio <- function(dealname, global.params, startdate=today()) {
+ dealdata <- getdealdata(dealname)
+ collatdata <- getcollateral(dealname)
+ notionalvec <- c()
+ SCvec <- c()
+ betavec <- c()
+ for(i in 1:nrow(collatdata)){
+ line.item <- collatdata[i,]
+ if( is.na(line.item$maturity) ){
+ stop("empty maturity")
+ }
+ ##most likely equity, doesn't impact the risk anyway
+ if(line.item$currentbalance < 1){
+ next
+ }
+ temp <- buildSC(line.item, dealdata$"Reinv End Date", dealdata$maturity, global.params, startdate)
+ notionalvec <- c(notionalvec, temp$notional)
+ SCvec <- c(SCvec, temp$SC)
+ betavec <- c(betavec, if(!is.na(line.item$iscdo) && line.item$iscdo) 1 else
+ global.params$defaultcorr)
+ }
+ return( list(notional=notionalvec, SC=SCvec, beta=betavec) )
+}
+
+cdrfromscenarios <- function(scenarios, dates){
+ ## compute the forward cdr rates to pass to intex
+ ## so that we match the default curves in scenarios
+ cdr <- matrix(0, nrow(scenarios), ncol(scenarios))
+ for(i in 1:nrow(scenarios)){
+ cdr[i,] <- 100*(1-exp(diff(c(0, log(1-scenarios[i,])))))/diff(c(0, yearFrac(today(), dates)))
+ }
+ return( cdr )
+}
+
+recoveryfromscenarios <- function(scenariosd, scenariosr){
+ ## compute the forward recovery rate based on the term
+ ## structure of recovery scenarios
+ ## we run into trouble for very stressed scenarios
+ ## this code should cap the scenarios at 0 it this happens
+ intexrecov <- matrix(0, n.scenarios, ncol(scenariosr))
+ for(i in 1:n.scenarios){
+ current <- 1
+ intexrecov[i,1] <- scenariosr[i,1]
+ for(t in 2:ncol(scenariosr)){
+ w <- scenariosd[i,current]/scenariosd[i,t]
+ if(scenariosr[i,t]-w*scenariosr[i,current]>=0){
+ intexrecov[i,t] <- (scenariosr[i,t]-w*scenariosr[i,current])/(1-w)
+ current <- current+1
+ }else{
+ intexrecov[i,t] <- 0
+ }
+ }
+ }
+ return(intexrecov)
+}
diff --git a/R/l1tf.R b/R/l1tf.R new file mode 100644 index 00000000..eab5208d --- /dev/null +++ b/R/l1tf.R @@ -0,0 +1,113 @@ +l1tf <- function(y,lambda){
+
+# PARAMETERS
+ALPHA = 0.01; # backtracking linesearch parameter (0,0.5]
+BETA = 0.5; # backtracking linesearch parameter (0,1)
+MU = 2; # IPM parameter: t update
+MAXITER = 40; # IPM parameter: max iteration of IPM
+MAXLSITER = 20; # IPM parameter: max iteration of line search
+TOL = 1e-4; # IPM parameter: tolerance
+
+# DIMENSIONS
+n = length(y); # length of signal x
+m = n-2; # length of Dx
+
+
+# OPERATOR MATRICES
+I2 = diag(n-2);
+O2 = matrix(0,n-2,1);
+D = cbind(I2,O2,O2) + cbind(O2-2*I2,O2) + cbind(O2,O2,I2);
+
+DDT = D%*%t(D);
+Dy = D%*%y;
+
+# VARIABLES
+z = rep(1,m); % dual variable
+mu1 = rep(1,m); % dual of dual variable
+mu2 = rep(1,m); % dual of dual variable
+
+t = 1e-10;
+pobj = Inf;
+dobj = 0;
+step = Inf;
+f1 = z-lambda;
+f2 = -z-lambda;
+
+for(iters in 0:MAXITER){
+
+ DTz = t(t(z)%*%D);
+ DDTz = D%*%DTz;
+ w = Dy-(mu1-mu2);
+
+ pobj1 = 0.5*t(w)%*%solve(DDT,w)+lambda*sum(mu1+mu2);
+ pobj2 = 0.5*t(DTz)%*%DTz+lambda*sum(abs(Dy-DDTz));
+ pobj = min(pobj1,pobj2);
+ dobj = -0.5*t(DTz)%*%DTz+t(Dy)%*%z;
+
+ gap = pobj - dobj;
+
+ # STOPPING CRITERION
+ if (gap <= TOL){
+ x = y-t(D)%*%z;
+ return;
+ }
+
+ if (step >= 0.2){
+ t =max(2*m*MU/gap, 1.2*t);
+ }
+
+ # CALCULATE NEWTON STEP
+
+ rz = DDTz - w;
+ S = DDT-sparse(1:m,1:m,mu1/f1+mu2/f2);
+ r = -DDTz + Dy + (1/t)/f1 - (1/t)/f2;
+ dz = solve(S,r);
+ dmu1 = -(mu1+((1/t)+dz*mu1)/f1);
+ dmu2 = -(mu2+((1/t)-dz*mu2)/f2);
+
+ resDual = rz;
+ resCent = rbind(-mu1*f1-1/t, -mu2*f2-1/t);
+ residual= rbind(resDual,resCent);
+
+ # BACKTRACKING LINESEARCH
+ negIdx1 = (dmu1 < 0);
+ negIdx2 = (dmu2 < 0);
+ step = 1;
+ if any(negIdx1){
+ step = min( step, 0.99*min(-mu1[negIdx1)/dmu1[negIdx1]) );
+ }
+ if (any(negIdx2){
+ step = min( step, 0.99*min(-mu2[negIdx2]/dmu2[negIdx2]) );
+ }
+
+ for( liter in 1:MAXLSITER){
+ newz = z + step*dz;
+ newmu1 = mu1 + step*dmu1;
+ newmu2 = mu2 + step*dmu2;
+ newf1 = newz - lambda;
+ newf2 = -newz - lambda;
+
+ # UPDATE RESIDUAL
+ newResDual = DDT*newz - Dy + newmu1 - newmu2;
+ newResCent = rbind(-newmu1.*newf1-1/t, -newmu2*newf2-1/t);
+ newResidual = rbind(newResDual, newResCent);
+
+ if ( max(max(newf1),max(newf2)) < 0 && ...
+ norm(newResidual) <= (1-ALPHA*step)*norm(residual) )
+ break;
+ }
+ step = BETA*step;
+ }
+ # UPDATE PRIMAL AND DUAL VARIABLES
+ z = newz; mu1 = newmu1; mu2 = newmu2; f1 = newf1; f2 = newf2;
+ }
+
+# The solution may be close at this point, but does not meet the stopping
+# criterion (in terms of duality gap).
+x = y-t(D)%*%z;
+if (iters >= MAXITER){
+ status = "maxiter exceeded";
+ cat(status,"\n");
+ return;
+ }
+
diff --git a/R/latestprices.R b/R/latestprices.R new file mode 100644 index 00000000..6394e695 --- /dev/null +++ b/R/latestprices.R @@ -0,0 +1,44 @@ +library(RODBC)
+conn <- odbcConnect("MLP-PROD")
+deltas.hist <- sqlQuery(conn,"select * from ET_CusipDeltasHist where IndexType='LCDX' and ReinvSpreadScenario='MID'")
+deltas.live <- sqlQuery(conn,"select * from ET_CusipDeltas where IndexType='LCDX' and ReinvSpreadScenario='MID'")
+deltas <- rbind(deltas.live,deltas.hist)
+deltas$PriceDate <- as.Date(deltas$PriceDate,format="%m/%d/%Y")
+deltas <- deltas[order(deltas$PriceDate,decreasing=T),]
+cusiplist <- unique(deltas$Cusip)
+deltas <- deltas[match(cusiplist,deltas$Cusip),]
+
+walduration <- sqlQuery(conn,"select * from ET_CusipSumProduct_WALDuration")
+prices_lcdx <- sqlQuery(conn,"select * from ET_CusipSumProduct_AWPrice where IndexType='LCDX'")
+prices_t1 <- sqlQuery(conn,"select * from ET_CusipSumProduct_AWPrice where IndexType='T1'")
+prices_t2 <- sqlQuery(conn,"select * from ET_CusipSumProduct_AWPrice where IndexType='T2'")
+prices_hy <- sqlQuery(conn,"select * from ET_CusipSumProduct_AWPrice where IndexType='HY' and IndexSeries='10' and IndexTenor='7y'")
+odbcClose(conn)
+selectlatest <- function(data){
+ data$PriceDate <- as.Date(data$PriceDate,format="%m/%d/%Y")
+ data <- data[order(data$PriceDate,data$UpdateDate,decreasing=T),]
+ data_mid <- data[data$ReinvSpreadScenario=="MID",]
+ data_mid <- data_mid[match(unique(data_mid$Cusip),data_mid$Cusip,),]
+ data_high <- data[data$ReinvSpreadScenario=="HIGH",]
+ data_high <- data_high[match(unique(data_high$Cusip),data_high$Cusip,),]
+ data_low <- data[data$ReinvSpreadScenario=="LOW",]
+ data_low <- data_low[match(unique(data_low$Cusip),data_low$Cusip,),]
+ data_noreinv <- data[data$ReinvSpreadScenario=="NOREINV",]
+ data_noreinv <- data_noreinv[match(unique(data_noreinv$Cusip),data_noreinv$Cusip,),]
+ data <- rbind(data_mid,data_noreinv,data_low,data_high)
+ data[order(data$Cusip),]
+}
+
+prices_t1 <- selectlatest(prices_t1)
+prices_t2 <- selectlatest(prices_t2)
+prices_hy <- selectlatest(prices_hy)
+prices_lcdx <- selectlatest(prices_lcdx)
+walduration <- selectlatest(walduration)
+prices_lcdx <- prices_lcdx[prices_lcdx$Cusip%in%prices_t1$Cusip,]
+prices_hy <- prices_hy[prices_hy$Cusip%in%prices_t1$Cusip,]
+deltas <- deltas[deltas$Cusip%in%prices_t1$Cusip,]
+walduration <- walduration[walduration$Cusip%in%prices_t1$Cusip,]
+deltas <- deltas[pmatch(prices_t1$Cusip,deltas$Cusip,dup=T),]
+data <- cbind(prices_lcdx[,c("Cusip","ReinvSpreadScenario","AWJuniorWgts_Price","AWSeniorWgts_Price")],prices_t1[,"AWJuniorWgts_Price"],prices_t2[,"AWJuniorWgts_Price"],prices_hy[,c("AWJuniorWgts_Price","AWSeniorWgts_Price")],walduration[,c("AW_WAL","AW_Duration","PriceDate")],deltas[,c("FullTranche","Tranche1","Tranche2","Tranche3","Tranche4")])
+colnames(data)<-c("Cusip","ReinvScen","LCDX12 5y Jr","LCDX12 5y Sr","T1","T2","HY10 7y Jr","HY10 7y Sr","Wal","Duration","Price Date","LCDX index","LCDX 0-8","LCDX 8-15","LCDX 15-30","LCDX 30-100")
+write.table(data,file="W:/CorpCDOs/latestprices_Prod.txt",row.names=F,col.names=T,sep=",")
diff --git a/R/load_bloomberg_data.R b/R/load_bloomberg_data.R new file mode 100644 index 00000000..b3dff18a --- /dev/null +++ b/R/load_bloomberg_data.R @@ -0,0 +1,49 @@ +library(RPostgreSQL)
+library(Rbbg)
+source('etdb.R')
+load("bloomberg_data.RData")
+cusips <- dbGetQuery(dbCon, "select distinct cusip from et_collateral")
+bbgCon <- blpConnect(throw.ticker.errors=FALSE)
+
+fields.corp <- c("PX_LAST","LAST_UPDATE_DT","ISSUER","MATURITY","CPN","CPN_TYP",
+ "CPN_FREQ","FLT_SPREAD","LIBOR_FLOOR","LN_CURRENT_MARGIN",
+ "LN_COVENANT_LITE","SECOND_LIEN_INDICATOR","DEFAULTED", "PRICING_SOURCE")
+fields.mtge <- c("LAST_UPDATE_DT", "ISSUER","MATURITY","CPN","CPN_TYP","CPN_FREQ","FLT_SPREAD","RTG_MOODY","RTG_MDY_INITIAL")
+
+
+
+secCorp <- paste(corpcusips, "Corp")
+dataCorp <- bdp(bbgCon, secCorp, fields.corp)
+corpcusips <- substr(rownames(dataCorp[which(!is.na(dataCorp$ISSUER)),]),1,9)
+
+dataCorp <- dataCorp[which(!is.na(dataCorp$ISSUER)),]
+rownames(dataCorp) <- substr(rownames(dataCorp), 1, 9)
+dataCorp <- data.frame(CUSIP=rownames(dataCorp), dataCorp )
+fh = file(paste0("//WDSENTINEL/share/CorpCDOs/data/bloomberg/bloomberg_datacorp_",
+Sys.Date(),".csv"),"wb")
+write.csv(dataCorp, file=fh, row.names=F)
+close(fh)
+
+
+chunkSize <- 5000
+dataMtge <- c()
+for(i in 1:(length(mtgecusips)/chunkSize)){
+ cusipSample <- mtgecusips[((i-1)*chunkSize+1):(i*chunkSize)]
+ secMtge <- paste(cusipSample, "Mtge")
+ dataMtge <- rbind(dataMtge, bdp(bbgCon, secMtge, fields.mtge))
+}
+i <- i+1
+secMtge <- paste(mtgecusips[((i-1)*chunkSize+1):length(mtgecusips)], "Mtge")
+dataMtge <- rbind(dataMtge, bdp(bbgCon, secMtge, fields.mtge))
+mtgecusips <- substr(rownames(dataMtge[which(!is.na(dataMtge$ISSUER)),]),1,9)
+dataMtge <- dataMtge[which(!is.na(dataMtge$ISSUER)),]
+rownames(dataMtge) <- substr(rownames(dataMtge), 1, 9)
+dataMtge <- data.frame(CUSIP=rownames(dataMtge), dataMtge )
+fh = file(paste0("//WDSENTINEL/share/CorpCDOs/data/bloomberg/bloomberg_datamtge_",
+Sys.Date(),".csv"),"wb")
+write.csv(dataMtge, file=fh, row.names=F)
+close(fh)
+
+
+## save(corpcusips, mtgecusips, dataCorp, dataMtge, file="bloomberg_data.RData")
+save(corpcusips, mtgecusips, dataCorp, dataMtge, file="bloomberg_data.RData")
diff --git a/R/loadcashflows.R b/R/loadcashflows.R new file mode 100644 index 00000000..26ea12f1 --- /dev/null +++ b/R/loadcashflows.R @@ -0,0 +1,38 @@ +cusips <- c("000743AA2","00083VAE1","009368AD3","00936BAA2","00936XAA4","03761KAG3","03761LAE6","039549AC4","30605LAA7","82626RAA0","92849RAA0")
+cfdir <- "T:/Analytics/Runs/20100608"
+test <- read.table(paste(cusip[1],"cashflow.txt",sep="_"),header=T)
+index_scen <- function(data,i){
+ return( which(data$"Scenario"==paste("SCENARIO",i,sep="")) )
+}
+
+as.Date(as.character(test[index_scen(test,1),]$Date),format="%Y%m%d")
+
+library(zoo)
+today <- function(){
+ as.Date(Sys.time())
+}
+df <- function(spreadcurve){
+ df <- exp(-cumsum(as.numeric(diff(c(spreadcurve$curvedate,time(spreadcurve$curve)))/365)*spreadcurve$curve))
+ list(curvedate=spreadcurve$curvedate,curve=zoo(df,time(spreadcurve$curve)))
+}
+test <- zoo(c(100,200,300),c(today+365*1:3))
+yc.usd <-list(curvedate=today(),curve=zoo(c(100,200,300),today()+365*1:3))
+
+addcurves <- function(curveA,curveB){
+ if(curveA$curvedate==curveB$curvedate){
+ curvemerged <- na.locf(na.locf(merge(curveA$curve,curveB$curve),na.rm=F,fromLast=T))
+ return(list(curvedate=curveA$curvedate,curve=curvemerged[,1]+curvemerged[,2]))
+ }else{
+ return(0)
+ }
+}
+multcurves <- function(curveA,curveB){
+ if(curveA$curvedate==curveB$curvedate){
+ curvemerged <- na.locf(na.locf(merge(curveA$curve,curveB$curve),na.rm=F,fromLast=T))
+ return(list(curvedate=curveA$curvedate,curve=curvemerged[,1]*curvemerged[,2]))
+ }else{
+ return(0)
+ }
+}
+spreadfrombondprice <- function(price,spread,maturity,yieldcurve){
+}
diff --git a/R/loan_universe.R b/R/loan_universe.R new file mode 100644 index 00000000..821022a3 --- /dev/null +++ b/R/loan_universe.R @@ -0,0 +1,70 @@ +source("cds_utils.R")
+source("cds_functions_generic.R")
+source("index_definitions.R")
+source("tranche_functions.R")
+source("yieldcurve.R")
+source("etdb.R")
+library(ggplot2)
+
+##source("optimization.R")
+MarkitData <- getMarkitIRData()
+L1m <- buildMarkitYC(MarkitData, dt = 1/12)
+L2m <- buildMarkitYC(MarkitData, dt = 1/6)
+L3m <- buildMarkitYC(MarkitData)
+L6m <- buildMarkitYC(MarkitData, dt = 1/2)
+setEvaluationDate(as.Date(MarkitData$effectiveasof))
+
+library(RQuantLib)
+loan.data <- dbGetQuery(dbCon, "select * from latest_markit_prices where lower(facility) in (select distinct lower(facility) as lf from latest_markit_prices where lower(facility) like 'tl%') order by maturity")
+flatshape <- splinefun(c(0,20),rep(1,2))
+loan.data.clean <- loan.data[loan.data$depth>=10,]
+loan.data.clean <- loan.data.clean[loan.data.clean$latestdate==as.Date("2012-09-11"),]
+loan.data.clean <- loan.data.clean[-(1:4),]
+SC.universe <- c()
+Rvec <- c()
+maturityvec <- as.Date(character(0))
+indexvec <- c()
+w1 <- c()
+w2 <- c()
+for(i in 1:nrow(loan.data.clean)){
+ if(i%%10==0){
+ cat(i, sep="\n")
+ }
+ collateral <- loan.data.clean[i,]
+ collateral$frequency <- "Q"
+ collateral$fixedorfloat <- "FLOAT"
+ collateral$grosscoupon <- 0
+ collateral$spread <- collateral$spread/100
+ collateral$price <- (collateral$bid+collateral$offer)/2
+ if(is.na(collateral$spread)){
+ next
+ }
+ ## negative hazard rate
+ if(collateral$price-100-
+ (yearFrac(today(), collateral$maturity) * collateral$spread)>=0){
+ next
+ }
+ if(collateral$price<=50){
+ next
+ }
+ recov <- min(collateral$price/100-0.2, 0.6)
+ sc <- bondhazardrate.shaped(collateral, flatshape, recov)
+ if(sc@h<=0){
+ next
+ }
+ SC.universe <- c(SC.universe, sc)
+ Rvec <- c(Rvec, recov)
+ maturityvec <- c(maturityvec, collateral$maturity)
+ indexvec <- c(indexvec, i)
+ w1 <- c(w1, collateral$amount)
+ w2 <- c(w2, collateral$amount * collateral$price/100)
+}
+
+hvec <- sapply(SC.universe, attr, "h")
+spreads <- hvec*(1-Rvec)
+T <- yearFrac(today(), maturityvec)
+test <- loess(spreads~T)
+
+loan.data2 <- loan.data.clean[!is.na(loan.data.clean$stm),]
+maturities <- yearFrac(today(),loan.data2$maturity)
+test2 <- loess(loan.data2$stm~maturities)
diff --git a/R/lossdistrib.c b/R/lossdistrib.c new file mode 100644 index 00000000..c873df9e --- /dev/null +++ b/R/lossdistrib.c @@ -0,0 +1,374 @@ +#include <R.h>
+#include <Rmath.h>
+#include <string.h>
+
+#define MIN(x, y) (((x) < (y)) ? (x) : (y))
+
+void lossdistrib(double *p, int *np, double *w, double *S, int *N, int* defaultflag, double *q) {
+ /* recursive algorithm with first order correction for computing
+ the loss distribution.
+ p vector of default probabilities
+ np length of p
+ S vector of severities (should be same length as p)
+ N number of ticks in the grid
+ defaultflat if true compute the default distribution
+ q the loss distribution */
+
+ int i, j, d1, d2, M;
+ double lu, d, p1, p2, sum;
+ double *qtemp;
+
+ lu = 1./(*N-1);
+ qtemp = Calloc(*N, double);
+ q[0] = 1;
+ M = 1;
+ for(i=0; i<(*np); i++){
+ d = (*defaultflag)? w[i]/lu : S[i] * w[i]/ lu;
+ d1 = floor(d);
+ d2 = ceil(d);
+ p1 = p[i] * (d2-d);
+ p2 = p[i] - p1;
+ memcpy(qtemp, q, MIN(M, *N) * sizeof(double));
+ for(j=0; j < MIN(M, *N); j++){
+ q[j] = (1-p[i]) * q[j];
+ }
+ for(j=0; j < MIN(M, *N-d2); j++){
+ q[d1+j] += p1 * qtemp[j];
+ q[d2+j] += p2 * qtemp[j];
+ };
+ M+=d2;
+ }
+
+ /* correction for weight loss */
+ if(M>*N){
+ sum = 0;
+ for(j=0; j<MIN(M, *N); j++){
+ sum += q[j];
+ }
+ q[*N-1] += 1-sum;
+ }
+ Free(qtemp);
+}
+
+void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N,
+ int *T, int *defaultflag, double *q) {
+ /* recursive algorithm with first order correction for computing
+ the loss distribution.
+ p vector of default probabilities
+ np length of p
+ S vector of severities (should be same length as p)
+ N number of ticks in the grid
+ T where to truncate the distribution
+ defaultflat if true computes the default distribution
+ q the loss distribution */
+
+ int i, j, d1, d2, M;
+ double lu, d, p1, p2, sum;
+ double *q1, *q2;
+
+ lu = 1./(*N-1);
+ q1 = Calloc( *T, double);
+ q2 = Calloc( *T, double);
+ q[0] = 1;
+ M = 1;
+ for(i=0; i<(*np); i++){
+ d = (*defaultflag)? w[i] / lu : S[i] * w[i] / lu;
+ d1 = floor(d);
+ d2 = ceil(d);
+ p1 = p[i] * (d2-d);
+ p2 = p[i] - p1;
+ for(j=0; j < MIN(M, *T); j++){
+ q1[j] = p1 * q[j];
+ q2[j] = p2 * q[j];
+ q[j] = (1-p[i]) * q[j];
+ }
+ for(j=0; j < MIN(M, *T-d1); j++){
+ q[d1+j] += q1[j];
+ };
+ for(j=0; j < MIN(M, *T-d2); j++){
+ q[d2+j] += q2[j];
+ };
+ M += d2;
+ }
+ Free(q1);
+ Free(q2);
+}
+
+void lossdistrib_joint( double *p, int *np, double *w, double *S, int *N, int *defaultflag, double *q) {
+ /* recursive algorithm with first order correction
+ computes jointly the loss and recovery distribution
+ p vector of default probabilities
+ np length of p
+ w vector of issuer weights (length np)
+ S vector of severities (should be same length as p)
+ N number of ticks in the grid
+ defaultflag if true computes the default distribution
+ q the joint probability distribution */
+
+ int i, j, k, m, n;
+ int Mx, My;
+ double lu, x, y, sum;
+ double alpha1, alpha2, beta1, beta2;
+ double w1, w2, w3, w4;
+ double *qtemp;
+
+ lu = 1./(*N-1);
+ qtemp = Calloc( (*N) * (*N), double);
+ q[0] = 1;
+ Mx=1;
+ My=1;
+ for(k=0; k<(*np); k++){
+ x = (*defaultflag)? w[k] /lu : S[k] * w[k] / lu;
+ y = (1-S[k]) * w[k] / lu;
+ i = floor(x);
+ j = floor(y);
+ alpha1 = i + 1 - x;
+ alpha2 = 1 - alpha1;
+ beta1 = j + 1 - y;
+ beta2 = 1 - beta1;
+ w1 = alpha1 * beta1;
+ w2 = alpha1 * beta2;
+ w3 = alpha2 * beta2;
+ w4 = alpha2 * beta1;
+
+ for(n=0; n<MIN(My, *N); n++){
+ memcpy(qtemp+n*(*N), q+n*(*N), MIN(Mx, *N) * sizeof(double));
+ }
+ for(n=0; n<MIN(My, *N); n++){
+ for(m=0; m<MIN(Mx, *N); m++){
+ q[m+(*N)*n] = (1-p[k])* q[m+(*N)*n];
+ }
+ }
+ for(n=0; n < MIN(My, *N-j-1); n++){
+ for(m=0; m < MIN(Mx, *N-i-1); m++){
+ q[i+m+(*N)*(j+n)] += w1 * p[k] * qtemp[m+(*N)*n];
+ q[i+m+(*N)*(j+1+n)] += w2 * p[k] * qtemp[m+(*N)*n];
+ q[i+1+m+(*N)*(j+1+n)] += w3 * p[k] * qtemp[m+(*N)*n];
+ q[i+1+m+(*N)*(j+n)] += w4 * p[k] *qtemp[m+(*N)*n];
+ }
+ }
+ Mx += i+1;
+ My += j+1;
+ }
+ /* correction for weight loss */
+ if(Mx>*N || My>*N){
+ sum = 0;
+ for(m=0; m < MIN(Mx, *N); m++){
+ for(n=0; n < MIN(My, *N); n++){
+ sum += q[m+n*(*N)];
+ }
+ }
+ q[MIN(*N, Mx)*MIN(My,*N)-1] += 1 - sum;
+ }
+ Free(qtemp);
+}
+
+void recovdist(double *dp, double *pp, int *n, double *w, double *S, int *N, double *q) {
+ /* recursive algorithm with first order correction for computing
+ the recovery distribution in case of prepayment.
+ dp vector of default probabilities
+ pp vector of prepay probabilities
+ n length of p
+ S vector of severities (should be same length as p)
+ w vector of weights
+ N number of ticks in the grid
+ q the loss distribution */
+
+ int i, j, d1l, d1u, d2l, d2u;
+ int M;
+ double lu, d1, d2, dp1, dp2, pp1, pp2, sum;
+ double *qtemp;
+
+ lu = 1./(*N - 1);
+ qtemp = Calloc( (*N), double);
+ q[0] = 1;
+ M=1;
+ for(i=0; i<(*n); i++){
+ d1 = w[i] * (1-S[i]) /lu;
+ d2 = w[i]/lu;
+ d1l = floor(d1);
+ d1u = d1l + 1;
+ d2l = floor(d2);
+ d2u = d2l + 1;
+ dp1 = dp[i] * (d1u - d1);
+ dp2 = dp[i] - dp1;
+ pp1 = pp[i] * (d2u - d2);
+ pp2 = pp[i] - pp1;
+ memcpy(qtemp, q, MIN(M, *N) * sizeof(double));
+ for(j = 0; j< MIN(M, *N); j++){
+ q[j] = (1-dp[i]-pp[i]) * q[j];
+ }
+ for(j=0; j < MIN(M, *N-d2u); j++){
+ q[d1l+j] += dp1 * qtemp[j];
+ q[d1u+j] += dp2 * qtemp[j];
+ q[d2l+j] += pp1 * qtemp[j];
+ q[d2u+j] += pp2 * qtemp[j];
+ };
+ M += d2u;
+ }
+ /* correction for weight loss */
+ if(M>*N){
+ sum = 0;
+ for(j=0; j<MIN(M, *N); j++){
+ sum += q[j];
+ }
+ q[*N-1] += 1-sum;
+ }
+ Free(qtemp);
+}
+
+void lossdistrib_prepay_joint(double *dp, double *pp, int *ndp, double *w,
+ double *S, int *N, int *defaultflag, double *q) {
+ int i, j1, j2, k, m, n;
+ double lu, x, y1, y2, sum;
+ double alpha1, alpha2, beta1, beta2;
+ double dpw1, dpw2, dpw3, dpw4;
+ double ppw1, ppw2, ppw3;
+ double temp;
+ double *qtemp;
+ int Mx, My;
+
+ lu = 1./(*N-1);
+ qtemp = Calloc((*N) * (*N), double);
+ q[0] = 1;
+ Mx=1;
+ My=1;
+ for(k=0; k<(*ndp); k++){
+ y1 = (1-S[k]) * w[k]/lu;
+ y2 = w[k]/lu;
+ x = (*defaultflag)? y2: y2-y1;
+ i = floor(x);
+ j1 = floor(y1);
+ j2 = floor(y2);
+ alpha1 = i + 1 - x;
+ alpha2 = 1 - alpha1;
+ beta1 = j1 + 1 - y1;
+ beta2 = 1 - beta1;
+ dpw1 = alpha1 * beta1 * dp[k];
+ dpw2 = alpha1 * beta2 * dp[k];
+ dpw3 = alpha2 * beta2 * dp[k];
+ dpw4 = alpha2 * beta1 * dp[k];
+
+ /* by default distribution, we mean names fractions of names that disappeared
+ either because of default or prepayment */
+ for(n=0; n<MIN(My, *N); n++){
+ memcpy(qtemp+n*(*N), q+n*(*N), MIN(Mx, *N) * sizeof(double));
+ }
+ for(n=0; n<MIN(My, *N); n++){
+ for(m=0; m<MIN(Mx, *N); m++){
+ q[m+(*N)*n] = (1-dp[k]-pp[k]) * q[m+(*N)*n];
+ }
+ }
+ if(*defaultflag){
+ ppw1 = alpha1 * alpha1 * pp[k];
+ ppw2 = alpha1 * alpha2 * pp[k];
+ ppw3 = alpha2 * alpha2 * pp[k];
+ for(n=0; n < MIN(My, *N-j2-1); n++){
+ for(m=0; m < MIN(Mx, *N-i-1); m++){
+ q[i+m+(*N)*(j1+n)] += dpw1 * qtemp[m+(*N)*n];
+ q[i+m+(*N)*(j1+1+n)] += dpw2 * qtemp[m+(*N)*n];
+ q[i+1+m+(*N)*(j1+1+n)] += dpw3 * qtemp[m+(*N)*n];
+ q[i+1+m+(*N)*(j1+n)] += dpw4 * qtemp[m+(*N)*n];
+
+ q[i+m+(*N)*(j2+n)] += ppw1 * qtemp[m+(*N)*n];
+ temp = ppw2 * qtemp[m+(*N)*n];
+ q[i+m+(*N)*(j2+1+n)] += temp;
+ q[i+1+m+(*N)*(j2+1+n)] += ppw3 * qtemp[m+(*N)*n];
+ q[i+1+m+(*N)*(j2+n)] += temp;
+ }
+ }
+ }else{
+ for(n=0; n < MIN(My, *N-j2-1); n++){
+ for(m=0; m < MIN(Mx, *N-i-1); m++){
+ q[i+m+(*N)*(j1+n)] += dpw1 * qtemp[m+(*N)*n];
+ q[i+m+(*N)*(j1+1+n)] += dpw2 * qtemp[m+(*N)*n];
+ q[i+1+m+(*N)*(j1+1+n)] += dpw3 * qtemp[m+(*N)*n];
+ q[i+1+m+(*N)*(j1+n)] += dpw4 * qtemp[m+(*N)*n];
+ q[m+(*N)*(j2+n)] += pp[k]*(j2+1-y2) * qtemp[m+(*N)*n];
+ q[m+(*N)*(j2+1+n)] += pp[k]*(y2-j2) * qtemp[m+(*N)*n];
+ }
+ }
+ }
+ Mx += i + 1;
+ My += j2 + 1;
+ }
+ /* correction for weight loss */
+ if(Mx>*N || My>*N){
+ sum = 0;
+ for(m=0; m < MIN(Mx, *N); m++){
+ for(n=0; n < MIN(My, *N); n++){
+ sum += q[m+n*(*N)];
+ }
+ }
+ q[MIN(*N, Mx)*MIN(My,*N)-1] += 1 - sum;
+ }
+ Free(qtemp);
+}
+
+double shockprob(double p, double rho, double Z, int give_log){
+ return( pnorm( (qnorm(p, 0, 1, 1, 0) -sqrt(rho) * Z)/sqrt(1 - rho), 0, 1, 1, give_log));
+}
+
+double shockseverity(double S, double Z, double rho, double p){
+ return( exp(shockprob(S * p, rho, Z, 1) - shockprob(p, rho, Z, 1)) );
+
+}
+
+void addandmultiply(double *X, double alpha, double *Y, int n) {
+ int i;
+ for(i = 0; i<n; i++){
+ Y[i] += alpha*X[i];
+ }
+}
+
+void BClossdist(double *SurvProb, int *dim1, int *dim2, double *issuerweights,
+ double *recov, double *Z, double *w, int *n, double *rho, int *N,
+ int *defaultflag, double *L, double *R) {
+ /*
+ computes the loss and recovery distribution over time with a flat gaussiancorrelation
+ inputs:
+ Survprob: matrix of size dim1 x dim2. dim1 is the number of issuers and dim2 number of time steps
+ recov: vector of recoveries (length dim1)
+ issuerweights: vector of issuer weights (length dim2)
+ Z: vector of factor values (length n)
+ w: vector of factor weights (length n)
+ rho: correlation beta
+ N: number of ticks in the grid
+ defaultflag: if true, computes the default distribution
+ outputs:
+ L: matrix of size (N, dim2)
+ R: matrix of size (N, dim2)
+ */
+ int t, i, j;
+ double g;
+ double *gshocked, *Rshocked, *Sshocked, *Lw, *Rw;
+
+ gshocked = Calloc((*dim1), double);
+ Rshocked = Calloc((*dim1), double);
+ Sshocked = Calloc((*dim1), double);
+ Lw = malloc((*N) * sizeof(double));
+ Rw = malloc((*N) * sizeof(double));
+
+ for(t=0; t < (*dim2); t++) {
+ for(i=0; i < *n; i++){
+ for(j=0; j < (*dim1); j++){
+ g = 1 - SurvProb[j + (*dim1) * t];
+ gshocked[j] = shockprob(g, *rho, Z[i], 0);
+ Sshocked[j] = shockseverity(1-recov[j], Z[i], *rho, g);
+ Rshocked[j] = 1 - Sshocked[j];
+ }
+ /* reset Lw and Rw to 0 */
+ memset(Lw, 0, *N * sizeof(double));
+ memset(Rw, 0, *N * sizeof(double));
+ lossdistrib(gshocked, dim1, issuerweights, Sshocked, N, defaultflag, Lw);
+ lossdistrib(gshocked, dim1, issuerweights, Rshocked, N, defaultflag, Rw);
+ addandmultiply(Lw, w[i], L + t * (*N), *N);
+ addandmultiply(Rw, w[i], R + t * (*N), *N);
+ }
+ }
+ Free(gshocked);
+ Free(Rshocked);
+ Free(Sshocked);
+ free(Lw);
+ free(Rw);
+}
diff --git a/R/mapping.R b/R/mapping.R new file mode 100644 index 00000000..2dfafe64 --- /dev/null +++ b/R/mapping.R @@ -0,0 +1,140 @@ +library(RPostgreSQL)
+options(stringsAsFactors=FALSE)
+drv <- dbDriver("PostgreSQL")
+dbCon <- dbConnect(drv, dbname="ET", user="et_user", password="Serenitas;1", host="192.168.1.108")
+prices <- dbGetQuery(dbCon, "select * from latest_markit_prices")
+collat <- dbGetQuery(dbCon, "select * from et_collateral")
+userLoanxid <- dbGetQuery(dbCon, "select * from loanx_user_mapping")
+userCUSIP <- dbGetQuery(dbCon, "select * from cusip_user_mapping")
+bbgmap <- read.csv(file="Z:/CorpCDOs/data/lx-bbgmap.csv", header = TRUE,colClasses = "character")
+RMBSmap <- read.csv(file="Z:/CorpCDOs/data/MTGECusips.csv", header = FALSE,colClasses = "character")
+
+prices$mid <-(prices$bid + prices$offer)/2
+prices$spread <- prices$spread/100
+
+from <- c(' LLC.?$', ', INC.?$', ' Inc.?$',' GMBH.?$', ' SA$', ' LTD$', ' PLC$',
+ ' CORP$', ' S A.?$', ' S.A.?$', ' BV$', ' AB$', ' LP$', ' N.V.?$',
+ ' AG$', ' CO.?$', ' HLDG.?$', ' HLDGS?$', ' S P A$', ' S.P.A.?$', ' A/S$', ' Holdings$'
+ ' B.V.$', '(', ')')
+prices$scrubedissuer <- prices$issuer
+collat$scrubedissuer <- collat$issuername
+for (i in 1:length(from)){
+ prices$scrubedissuer <- gsub(from[i], '', prices$scrubedissuer, ignore.case = TRUE)
+ collat$scrubedissuer <- gsub(from[i], '', collat$scrubedissuer, ignore.case = TRUE)
+}
+
+#collat1 <- collat
+collat1 <- collat[collat$dealname == "abrlf",]
+#collat1$src[!is.na(collat1$price)] <- "IDMatch"
+#collat1$src[collat1$cusip %in% RMBSmap$V1] <- "RMBS"
+#collat1$iscdo[collat1$cusip %in% CLOmap$cusip] <- TRUE
+
+for (i in 1:length(collat1$issuername)){
+ #Match CUSIP vs. LoanxID
+ matched0 <- grep(collat1$cusip[i], bbgmap$cusip, ignore.case=TRUE)
+ if (is.na(collat1$loanxid[i])&!is.na(collat1$cusip[i])&length(matched0)==1) {
+ collat1$loanxid[i] <- bbgmap$loanxid[matched0]
+ matched0a <- grep(collat1$loanxid[i], prices$loanxid, ignore.case=TRUE)
+ if(length(matched0a)==1) collat1$price[i] <- prices$mid[matched0a]
+ collat1$src[i]<- "CUSIPMapped"
+ }
+ if (is.na(collat1$loanxid[i])){
+ collat1$et_loanxid[i] <- userLoanxid$loanxid[collat1$scrubedissuer[i] == userLoanxid$issuername & collat1$maturity[i] == userLoanxid$maturity & collat1$spread[i] == userLoanxid$spread]
+ collat1$src[i]<- "UserMapped"
+ }
+ if (is.na(collat1$loanxid[i])&is.na(collat1$et_loanxid[i])) {
+ #matchN1 <- grep(collat1$scrubedissuer[i], prices$scrubedissuer, ignore.case=TRUE)
+ #matchN1 <- prices$scrubedissuer %in% matchN1
+ maxchange <- as.integer(nchar(collat1$scrubedissuer[i])*.25)#25% of string length
+ z <- data.frame(x=as.numeric(adist(collat1$scrubedissuer[i], prices$scrubedissuer, ignore.case = TRUE)),y=prices$scrubedissuer)
+ z <- z[order(z$x),]
+ z <- z[z$x<=maxchange,]$y
+ matchN1 <- prices$scrubedissuer %in% z
+
+ matchS1 <- prices$spread == collat1$spread[i] #Match Spread
+ matchS2 <- abs(prices$spread - collat1$spread[i]) <=.5 #Match Spread+/-50bps
+ matchS3 <- abs(prices$spread - collat1$spread[i]) <=1 #Match Spread+/-100bps
+
+ matchD1 <- abs(prices$maturity - collat1$maturity[i]) <=2 #Match Maturity 2 day diff
+ matchD2 <- abs(prices$maturity - collat1$maturity[i]) <=30 #Match Maturity <=30 days diff
+ matchD3 <- abs(prices$maturity - collat1$maturity[i]) <=200 #Match Maturity <=200 days diff
+
+ CompMatchV1NA <- matchN1 & matchS1 & matchD1
+ CompMatchV2NA <- matchN1 & matchS2 & matchD1
+ CompMatchV3NA <- matchN1 & matchS2 & matchD2
+ CompMatchV4NA <- matchN1 & matchS3 & matchD1
+ CompMatchV5NA <- matchN1 & matchS1 & matchD3
+ CompMatchV6NA <- matchN1 & matchS1 & matchD3
+
+ CompMatchV1NA[is.na(CompMatchV1NA)] <- FALSE
+ CompMatchV2NA[is.na(CompMatchV2NA)] <- FALSE
+ CompMatchV3NA[is.na(CompMatchV3NA)] <- FALSE
+ CompMatchV4NA[is.na(CompMatchV4NA)] <- FALSE
+ CompMatchV5NA[is.na(CompMatchV5NA)] <- FALSE
+
+ CompMatchV1 <- CompMatchV1NA
+ CompMatchV2 <- CompMatchV2NA
+ CompMatchV3 <- CompMatchV3NA
+ CompMatchV4 <- CompMatchV4NA
+ CompMatchV5 <- CompMatchV5NA
+
+ updatetable <- FALSE
+
+ #Matches Found, Randomly picks the first one
+ if (any(CompMatchV1)){
+ collat1$et_loanxid[i] <- prices$loanxid[which(CompMatchV1)][1]
+ updatetable <- TRUE
+ collat1$price[i] <- prices$mid[which(CompMatchV1)][1]
+ collat1$src[i]<- "MatchedNameSpreadDate"
+ } else if (any(CompMatchV2)) {
+ collat1$et_loanxid[i] <- prices$loanxid[which(CompMatchV2)][1]
+ updatetable <- TRUE
+ collat1$price[i] <- prices$mid[which(CompMatchV2)][1]
+ collat1$src[i]<- "MatchedNameSpread+-50bpsDate"
+ } else if (any(CompMatchV3)) {
+ collat1$et_loanxid[i] <- prices$loanxid[which(CompMatchV3)][1]
+ updatetable <- TRUE
+ collat1$price[i] <- prices$mid[which(CompMatchV3)][1]
+ collat1$src[i]<- "MatchedNameSpread+-50bpsDate+-30Days"
+ } else if (any(CompMatchV4)) {
+ collat1$et_loanxid[i] <- prices$loanxid[which(CompMatchV4)]
+ updatetable <- TRUE
+ collat1$price[i] <- prices$mid[which(CompMatchV4)][1]
+ collat1$src[i]<- "MatchedNameSpread+-100bpsDate"
+ } else if (any(CompMatchV5)) {
+ collat1$et_loanxid[i] <- prices$loanxid[which(CompMatchV5)]
+ updatetable <- TRUE
+ collat1$price[i] <- prices$mid[which(CompMatchV5)][1]
+ collat1$src[i]<- "MatchedNameSpreadDate+-200Days"
+ }
+ if (updatetable == true) {
+ str <- paste0("a.dealname = '", collat1$dealname[i], "' and a.name = '", collat1$name[i], "' and a.maturity = '", collat1$maturity[i], "'")
+ dbGetQuery(dbCon, paste0("UPDATE et_collateral a SET et_loanxid = '", collat1$et_loanxid[i],"' WHERE ",str))
+ }
+ }
+}
+
+#collat1$loanxid[collat1$cusip[i] %in% bbgmap$cusip]
+#rm(collat)
+#gc()
+
+
+
+#write.csv(collat, file="Z:/edwin/collat.csv")
+write.csv(collat1, file="Z:/edwin/collat1.csv")
+write.csv(prices, file="Z:/edwin/prices.csv")
+
+#scrubed <- prices[,c("scrubedissuer", "issuer")]
+#scrubed1 <- collat[,c("scrubedissuer", "issuername")]
+#write.csv(scrubed, file="Z:/edwin/scrubedissusername.csv")
+#write.table(scrubed,"clipboard", sep = "\t")
+
+#pD1 <- merge(collat, prices, by.x=c("loanxid"), by.y=c("loanxid"))
+#pD1$src <- "LXID"
+pD2 <- merge(pD1, prices, by.x=c("scrubedissuer.x", "spread.x","maturity.x"),
+ by.y=c("scrubedissuer", "spread","maturity"))
+pD2$src <- "ExactMatch"
+
+pD3[1:50,c("issuername", "loanxid.y", "loanxid.x")]
+## pD2 <- pD1[c("issuername","loanxid","mid")]
+
diff --git a/R/mapping_fast.R b/R/mapping_fast.R new file mode 100644 index 00000000..1692f3fd --- /dev/null +++ b/R/mapping_fast.R @@ -0,0 +1,62 @@ +library(RPostgreSQL)
+library(data.table)
+library(hash)
+options(stringsAsFactors=FALSE)
+drv <- dbDriver("PostgreSQL")
+dbCon <- dbConnect(drv, dbname="ET", user="et_user", password="Serenitas;1", host="192.168.1.108")
+markitdata <- dbGetQuery(dbCon, "select * from latest_markit_prices")
+CLOmap <- dbGetQuery(dbCon, "select * from dealcusipmapping")
+intexdata <- dbGetQuery(dbCon, sprintf("select * from et_historical_collateral('%s')",Sys.Date()))
+
+bbgmap <- read.csv(file="//WDSENTINEL/share/CorpCDOs/data/lx-bbgmap.csv", header = TRUE, colClasses = "character")
+bbgmap <- data.table(bbgmap, key="cusip")
+RMBSmap <- read.csv(file="//WDSENTINEL/share/CorpCDOs/data/MTGECusips.csv", header = FALSE, colClasses = "character")
+markitdata$mid <-(markitdata$bid + markitdata$offer)/2
+markitdata$spread <- markitdata$spread/100
+
+toscrub <- c(' LLC.?$', ', INC.?$', ' Inc.?$',' GMBH.?$', ' SA$', ' LTD$', ' PLC$',
+ ' CORP$', ' S A.?$', ' S.A.?$', ' BV$', ' AB$', ' LP$', ' N.V.?$',
+ ' AG$', ' CO.?$', ' HLDG.?$', ' HLDGS?$', ' S P A$', ' S.P.A.?$',
+ ' A/S$', ' Holdings$', ' B.V.$',"\\(", "\\)")
+intexdata$scrubedissuer <- tolower(intexdata$issuername)
+markitdata$scrubedissuer <- tolower(markitdata$issuer)
+for (i in 1:length(toscrub)){
+ markitdata$scrubedissuer <- gsub(toscrub[i], '', markitdata$scrubedissuer, ignore.case = TRUE)
+ intexdata$scrubedissuer <- gsub(toscrub[i], '', intexdata$scrubedissuer, ignore.case = TRUE)
+}
+markitdata <- data.table(markitdata, key="scrubedissuer")
+
+intexdatatomap <- intexdata[is.na(intexdata$loanxid),]
+intexdatatomap <- data.frame(intexdatatomap,src=rep(NA,nrow(intexdatatomap)))
+
+a <- unique(intexdatatomap$scrubedissuer)
+b <- unique(markitdata$scrubedissuer)
+levdist <- adist(a, b, costs=c("ins"=1, "del"=1, "sub"=3))
+ratio <- (outer(nchar(a), nchar(b), "+")-levdist)/outer(nchar(a), nchar(b), "+")
+index <- hash(a,1:length(a))
+rownames(ratio) <- a
+
+for (i in 1:nrow(intexdatatomap)){
+ #Match CUSIP vs. LoanxID
+ matched0 <- bbgmap[intexdatatomap$cusip[i]]$loanxid
+ if(!is.na(matched0[1])){
+ intexdatatomap$loanxid[i] <- matched0[1]
+ intexdatatomap$src[i] <- "CUSIPMapped"
+ next
+ }
+
+ issuermatch <- b[ratio[index[[intexdatatomap$scrubedissuer[i]]],]>=0.9]
+ if( length(issuermatch) == 0){
+ next
+ }else{
+ issuermatch <- markitdata[issuermatch]
+ issuermatch <- issuermatch[!is.na(issuermatch$spread) & !is.na(issuermatch$maturity),]
+ spreadmatch <- abs(issuermatch$spread-intexdatatomap$spread[i])<1
+ maturitymatch <- abs(issuermatch$maturity-intexdatatomap$maturity[i])<90
+ issuermatch <- issuermatch[spreadmatch & maturitymatch,]
+ if(nrow(issuermatch)>0){
+ intexdatatomap$loanxid[i] <- issuermatch$loanxid[1]
+ intexdatatomap$src[i] <- "loanxidmatched"
+ }
+ }
+}
diff --git a/R/optimization.R b/R/optimization.R new file mode 100644 index 00000000..8f662b4e --- /dev/null +++ b/R/optimization.R @@ -0,0 +1,440 @@ +#various optimization functions
+library(Rglpk)
+
+"blkdiag" <-
+function (m1, m2, p.tr = 0, p.ll = 0)
+{
+ ## p.tr and p.ll are padding values
+ topleft <- m1
+ topright <- matrix(p.tr, nrow(m1), ncol(m2))
+ colnames(topright) <- colnames(m2)
+ lowleft <- matrix(p.ll, nrow(m2), ncol(m1))
+ lowright <- m2
+ rbind(cbind(topleft, topright), cbind(lowleft, lowright))
+}
+
+linfit<-function(A, q, gamma){
+ #minimize ||Ax-q||_1 + gamma||DAx||_1
+ #1^Tx=1
+ #x>=0
+
+ n<-nrow(P)
+ p<-ncol(P)
+ b<-rep(0,n)
+ r<-length(q)
+ obj<-c(rep(0,p), rep(1,r), rep(gamma,r-1))
+ D<-diff(diag(r))
+ ineqconst<-rbind(cbind(A, -diag(r), matrix(0,r,r-1)),
+ cbind(A, diag(r), matrix(0,r,r-1)),
+ cbind(D%*%A, matrix(0,r-1,r), -diag(r-1)),
+ cbind(D%*%A, matrix(0,r-1,r), diag(r-1)))
+ const.mat <- ineqconst
+ const.dir <- c(rep("<=",r), rep(">=",r), rep("<=",r), rep(">=",r))
+ const.rhs <- c(1,q,q,rep(0,r),rep(0,r))
+ model <- Rglpk_solve_LP(obj,const.mat,const.dir,const.rhs)
+ list(weight=model$solution[1:p], model=model)
+}
+
+linfit1<-function(P, b, A, q, gamma){
+ #minimize ||Ax-q||_1 + gamma||DAx||_1
+ # s.t. Px = b
+ #1^Tx=1
+ #x>=0
+
+ n<-nrow(P)
+ p<-ncol(P)
+ r<-length(q)
+ obj<-c(rep(0,p),rep(1,r), rep(gamma,r-1))
+ D<-diff(diag(r))
+ linconst<-rbind(cbind(P, matrix(0,n,2*r-1)),
+ c(rep(1,p), rep(0,2*r-1)))
+ ineqconst<-rbind(cbind(A,-diag(r), matrix(0,r,r-1)),
+ cbind(A,diag(r), matrix(0,r,r-1)),
+ cbind(D%*%A, matrix(0,r-1,r), -diag(r-1)),
+ cbind(D%*%A, matrix(0,r-1,r), diag(r-1)))
+ const.mat<-rbind(linconst, ineqconst)
+ const.dir<-c(rep("==",n+1), rep("<=",r), rep(">=",r), rep("<=",r), rep(">=",r))
+ const.rhs<-c(b, 1, q, q, rep(0,r), rep(0,r))
+ model <- Rglpk_solve_LP(obj, const.mat, const.dir, const.rhs)
+ list(weight=model$solution[1:p], model=model)
+}
+
+linfit2<-function(P, A, q, gamma){
+ #minimize ||Ax-p||_1 + gamma||DAx||_1
+ # s.t. Px=0
+ #1^Tx=1
+ #x>=0
+
+ n <- nrow(P)
+ p <- ncol(P)
+ b <- rep(0,n)
+ r <- length(q)
+ obj <- c(rep(0,p),rep(1,r), rep(gamma,r-2))
+ D <- diff(diff(diag(r)))
+ linconst <- rbind(cbind(P,matrix(0,n,2*r-2)),
+ c(rep(1,p),rep(0,2*r-2)))
+ ineqconst <- rbind(cbind(A,-diag(r),matrix(0,r,r-2)),
+ cbind(A,diag(r),matrix(0,r,r-2)),
+ cbind(D%*%A,matrix(0,r-2,r),-diag(r-2)),
+ cbind(D%*%A,matrix(0,r-2,r),diag(r-2)))
+ const.mat <- rbind(linconst,ineqconst)
+ const.dir <- c(rep("==",n+1), rep("<=",r), rep(">=",r), rep("<=",r), rep(">=",r))#"
+ const.rhs <- c(b,1,q,q,rep(0,r),rep(0,r))
+ model <- Rglpk_solve_LP(obj, const.mat, const.dir, const.rhs)
+ list(weight=model$solution[1:p], model=model)
+}
+
+
+minfit<-function(P){
+ #minimize ||Px||_1
+ #subject to x>=0
+ # 1^Tx=1
+
+ require(Rglpk)
+ n<-nrow(P)
+ p<-ncol(P)
+ b<-rep(0,n)
+ obj<-c(rep(0,p),1)
+ linconst<-t(c(rep(1,p),0))
+ ineqconst<-rbind(cbind(P,-rep(1,n)),
+ cbind(P,rep(1,n)))
+ const.mat<-rbind(linconst,ineqconst)
+ const.dir<-c("==",rep("<=",n),rep(">=",n))#"
+ const.rhs<-c(1,rep(0,2*n))
+ model <- Rglpk_solve_LP(obj,const.mat,const.dir,const.rhs)
+ list(weight=model$solution[1:p],model=model)
+}
+
+quadfit<-function(A, P, q, b=rep(0,nrow(P))){
+ #solves for x, min||Ax-q||^2
+ # s.t. Px=b
+ # x>=0
+ # sum(x)=1
+ alpha<-0.4
+ beta<-0.8
+ t<-1/alpha
+ b<-c(b,1)
+ P<-rbind(P,rep(1,ncol(P)))
+ attach(svd(A))
+ #init x and nu
+ x<-rep(1,ncol(P))/ncol(P)
+ nu<-rep(1,nrow(P))
+ for( i in 1:10){
+ t <- alpha*t
+ decr <- Inf
+ while(decr>1e-6){
+ #solve the KKT equation
+ rdual <- t(A)%*%A%*%x-2*t(t(q)%*%A)-t/x+t(P)%*%nu
+ rprimal <- P%*%x-b
+ #inverse of H=A^TA-t*diag(1/w^2) via woodbury
+ vainv <- sweep(t(v),2,x^2,"*")#"
+ Hinv <- 1/t*diag(as.numeric(x)^2)-1/t^2*sweep(v,1,x^2,"*")%*%solve(diag(1/d)^2+1/t*vainv%*%v,vainv)#"
+ S <- -P%*%Hinv%*%t(P)
+ Dnu_pd <- solve(S,P%*%Hinv%*%rdual-rprimal)
+ Dx_pd <- -Hinv%*%(t(P)%*%Dnu_pd+rdual)
+ #backtracking search
+ phi <- 1
+ newx <- x+phi*Dx_pd
+ newnu <- nu+phi*Dnu_pd
+ newrdual <- t(A)%*%A%*%newx-2*t(t(q)%*%A)-t/newx+t(P)%*%newnu
+ newrprimal <- P%*%newx-b
+
+ while((sqrt(sum(newrdual^2)+sum(newrprimal^2))>(1-alpha*phi)*sqrt(sum(rdual^2)+sum(rprimal^2))) || (sum(newx<0)>0)){
+ phi <- phi*beta
+ newx <- x+phi*Dx_pd
+ newnu <- nu+phi*Dnu_pd
+ newrdual <- t(A)%*%A%*%newx-2*t(t(q)%*%A)-t/newx+t(P)%*%newnu
+ newrprimal <- P%*%newx-b
+ }
+ x <- newx
+ nu <- newnu
+ decr <- sqrt(sum(newrdual^2)+sum(newrprimal^2))
+ cat(decr,"\n")
+ }
+ }
+ return(x)
+}
+
+
+KLfit<-function(P, q, b=rep(0,nrow(P))){
+ #solves for x min KL(x,q)
+ # s.t. Px=b
+ # x>=0
+ # sum(x)=1
+ alpha <- 0.4
+ beta <- 0.8
+ b <- c(b,1)
+ P <- rbind(P,rep(1,ncol(P)))
+ #init x and nu
+ x <- rep(1,ncol(P))/ncol(P)
+ nu <- rep(1,nrow(P))
+ decr <- Inf
+ eps <- 1e-12
+ niter <- 1
+ q <- (q+eps)/sum(q+eps)
+ while(decr>eps & niter<500){
+ #solve the KKT equation
+ rdual <- 1+log(x)-log(q)+t(P)%*%nu
+ rprimal <- P%*%x-b
+ S <- -P%*%diag(as.numeric(x))%*%t(P)
+ Dnu_pd <- solve(S,P%*%diag(as.numeric(x))%*%rdual-rprimal)
+ Dx_pd <- -diag(as.numeric(x))%*%(t(P)%*%Dnu_pd+rdual)
+ #backtracking search
+ phi <- 1
+ newx <- x+phi*Dx_pd
+ newnu <- nu+phi*Dnu_pd
+
+ while(sum(newx<0)>0){
+ phi <- phi*beta
+ newx <- x+phi*Dx_pd
+ newnu <- nu+phi*Dnu_pd
+ }
+ newrdual <- 1+log(newx)-log(q)+t(P)%*%newnu
+ newrprimal <- P%*%newx-b
+ while((sqrt(sum(newrdual^2)+sum(newrprimal^2))>(1-alpha*phi)*sqrt(sum(rdual^2)+sum(rprimal^2)))){
+ phi <- phi*beta
+ newx <- x+phi*Dx_pd
+ newnu <- nu+phi*Dnu_pd
+ newrdual <- 1+log(newx)-log(q)+t(P)%*%newnu
+ newrprimal <- P%*%newx-b
+ }
+ x <- newx
+ nu <- newnu
+ decr <- sqrt(sum(newrdual^2)+sum(newrprimal^2))
+ niter <- niter+1
+ #cat(decr,"\n")
+ }
+ return(list(obj=sum(x*log(x/q)), weight=x, status=decr))
+}
+
+KLfitmod<-function(P,q,b=rep(0,nrow(P)),lambda,H,x=rep(1,ncol(P))/ncol(P)){
+ #solves for x min KL(x,q)+lambda*<x,H^2>
+ # s.t. Px=b
+ # x>=0
+ # sum(x)=1
+ alpha<-0.4
+ beta<-0.8
+ b <- c(b,1)
+ P <- rbind(P,rep(1,ncol(P)))
+ #init x and nu
+ nu <- rep(1,nrow(P))
+ decr <- Inf
+ t <- 1
+ while(decr>1e-6){
+ #solve the KKT equation
+ rdual <- 1+log(x)-log(q)+lambda*H^2+t(P)%*%nu
+ rprimal<-P%*%x-b
+ S <- -P%*%diag(as.numeric(x))%*%t(P)
+ Dnu_pd <- solve(S,P%*%diag(as.numeric(x))%*%rdual-rprimal)
+ Dx_pd <- -diag(as.numeric(x))%*%(t(P)%*%Dnu_pd+rdual)
+ #backtracking search
+ phi<-1
+ newx<-x+phi*Dx_pd
+ newnu<-nu+phi*Dnu_pd
+
+ while(sum(newx<0)>0){
+ phi<-phi*beta
+ newx<-x+phi*Dx_pd
+ newnu<-nu+phi*Dnu_pd
+ }
+ newrdual<-1+log(newx)-log(q)+lambda*H^2+t(P)%*%newnu
+ newrprimal<-P%*%newx-b
+ while((sqrt(sum(newrdual^2)+sum(newrprimal^2))>(1-alpha*phi)*sqrt(sum(rdual^2)+sum(rprimal^2)))){
+ phi<-phi*beta
+ newx<-x+phi*Dx_pd
+ newnu<-nu+phi*Dnu_pd
+ newrdual<-1+log(newx)-log(q)+lambda*H^2+t(P)%*%newnu
+ newrprimal<-P%*%newx-b
+ }
+ x<-newx
+ nu<-newnu
+ decr<-sqrt(sum(newrdual^2)+sum(newrprimal^2))
+ decr<-sum(newrdual^2)
+ #cat(decr,"\n")
+ }
+ return(x)
+}
+
+KLfit2<-function(P,q,b=rep(0,nrow(P))){
+ #solves for x min KL(q,x)
+ # s.t. Px=b
+ # x>=0
+ # sum(x)=1
+ alpha<-0.4
+ beta<-0.8
+ b<-c(b,1)
+ P<-rbind(P,rep(1,ncol(P)))
+ #init x and nu
+ x<-rep(1,ncol(P))/ncol(P)
+ nu<-rep(1,nrow(P))
+ decr<-Inf
+ t<-1
+ while(decr>1e-6){
+ #solve the KKT equation
+ rdual<--q/x+t(P)%*%nu
+ rprimal<-P%*%x-b
+ S <- -P%*%diag(as.numeric(x^2/q))%*%t(P)
+ Dnu_pd <- solve(S,P%*%diag(as.numeric(x^2/q))%*%rdual-rprimal)
+ Dx_pd<--diag(as.numeric(x^2/q))%*%(t(P)%*%Dnu_pd+rdual)
+ #backtracking search
+ phi<-1
+ newx<-x+phi*Dx_pd
+ newnu<-nu+phi*Dnu_pd
+
+ while(sum(newx<0)>0){
+ phi<-phi*beta
+ newx<-x+phi*Dx_pd
+ newnu<-nu+phi*Dnu_pd
+ }
+ newrdual<--q/newx+t(P)%*%newnu
+ newrprimal<-P%*%newx-b
+ while((sqrt(sum(newrdual^2)+sum(newrprimal^2))>(1-alpha*phi)*sqrt(sum(rdual^2)+sum(rprimal^2)))){
+ phi<-phi*beta
+ newx<-x+phi*Dx_pd
+ newnu<-nu+phi*Dnu_pd
+ newrdual<--q/newx+t(P)%*%newnu
+ newrprimal<-P%*%newx-b
+ }
+ x<-newx
+ nu<-newnu
+ decr<-sqrt(sum(newrdual^2)+sum(newrprimal^2))
+ decr<-sum(newrdual^2)
+ cat(decr,"\n")
+ }
+ return(x)
+}
+
+solveKKT<-function(H,A,g,h,inv){
+ #solves the KKT system
+ #inv(H,v) should solve H^-1g
+ S <- -t(A)%*%inv(H,A)
+ if(ncol(A)==1){
+ Dy <- 1/S*(h-t(A)%*%inv(H,g))
+ }else{
+ Dy <- solve(S,h-t(A)%*%inv(H,g))
+ }
+ Dx <- inv(H,g-A%*%Dy)
+ return(list(Dx=as.vector(Dx),Dy=as.vector(Dy)))
+}
+
+lmconst<-function(Y,X,A,b,pos){
+ #solves for \beta min ||Y-X\beta||
+ # s.t. A\beta = b
+ # \beta>=0 (if pos = TRUE)
+ alpha<-0.4
+ beta<-0.8
+ flag<-1
+ #init x and nu
+ x<-rep(1,ncol(A))
+ nu<-rep(1,nrow(A))
+
+ t<-pos
+ XTX<-crossprod(X,X)
+ while(t >1e-10 || flag){
+ t <- alpha*t
+ flag <- 0
+ decr<-Inf
+ niter<-1
+ while(decr>1e-6 & niter<100){
+ #solve the KKT equation
+ rdual<--crossprod(X,Y-X%*%x)+crossprod(A,nu)-t/x
+ rprimal<-A%*%x-b
+ H<-XTX+t*diag(1/x^2)
+ KKT<-solveKKT(H,t(A),-rdual,-rprimal,invchol)
+ #KKT<-rbind(cbind(H,t(A)),cbind(A,matrix(0,nrow(A),nrow(A))))
+ #temp<--solve(KKT,rbind(rdual,rprimal))
+ #Dx_pd<-temp[1:nrow(H)]
+ #Dnu_pd<-temp[-(1:nrow(H))]
+ Dnu_pd<-KKT$Dy
+ Dx_pd<-KKT$Dx
+ #backtracking search
+ phi<-1
+ newx<-x+phi*Dx_pd
+ newnu<-nu+phi*Dnu_pd
+ if( pos==1 ){
+ while(min(newx)<0){
+ phi<-phi*beta
+ newx<-x+phi*Dx_pd
+ newnu<-nu+phi*Dnu_pd
+ }
+ }
+ newrdual<--crossprod(X,Y-X%*%newx)+crossprod(A,newnu)-t/newx
+ newrprimal<-A%*%newx-b
+ while((sqrt(sum(newrdual^2)+sum(newrprimal^2))>(1-alpha*phi)*sqrt(sum(rdual^2)+sum(rprimal^2)))){
+ phi<-phi*beta
+ newx<-x+phi*Dx_pd
+ newnu<-nu+phi*Dnu_pd
+ newrdual<--crossprod(X,Y-X%*%newx)+crossprod(A,newnu)-t/newx
+ newrprimal<-A%*%newx-b
+ }
+ x<-newx
+ nu<-newnu
+ decr<-sqrt(sum(newrdual^2)+sum(newrprimal^2))
+ #cat(decr,"\n")
+ niter<-niter+1
+ }
+ }
+ return(list(obj=norm(Y-X%*%x),beta=x))
+}
+
+#cleanup cdf
+cleancdf<-function(x,lambda){
+ #minimize for a ||a-x||_1+lambda*||diff(x)||_1
+ #subject to diff(x)>=0
+ r<-length(x)
+ obj<-c(rep(0,r),rep(1,r), rep(lambda,r-1))
+ D<-diff(diag(r))
+ const.mat<-rbind(cbind(diag(r),diag(r),matrix(0,r,r-1)),
+ cbind(D,matrix(0,r-1,r),matrix(0,r-1,r-1)),
+ cbind(D,matrix(0,r-1,r),-diag(r-1)),
+ cbind(diag(r),-diag(r),matrix(0,r,r-1) ))
+
+ const.dir<-c(rep(">=",2*r-1),rep("<=",2*r-1))#"
+ const.rhs<-c(x,rep(0,r-1),rep(0,r-1),x)
+ model <- Rglpk_solve_LP(obj,const.mat,const.dir,const.rhs)
+ model$solution[1:r]
+}
+
+posls<-function(Y,X,w=rep(1,nrow(X))){
+ #minimize ||Y-X\beta||_w
+ #subject to beta>=0
+ beta <- rep(1,ncol(X))
+ t <- 1
+ phi <- 0.6
+ for(i in 1:30){
+ dbeta<-rep(Inf,ncol(X))
+ while(norm(dbeta)>=1e-10){
+ g<--crossprod(X,w*(Y-X%*%beta))-t/beta
+ H<-crossprod(X,sweep(X,1,w,"*"))+diag(as.vector(t/beta^2))
+ dbeta<-solve(H,g)
+ while(min(beta-dbeta)<0){
+ dbeta<-phi*dbeta
+ }
+ beta<-beta-dbeta
+ cat(norm(dbeta),"\n")
+ }
+ t<-phi*t
+ }
+ return( beta )
+}
+
+invchol <- function(H,b){
+ R <- chol(H)
+ z <- forwardsolve(t(R),b)
+ x <- backsolve(R,z)
+ return(x)
+}
+
+#solve the KKT equation by block elemination
+solveswm <- function(A,B,C,b){
+ #solve Dx=b for D=A+BC with A diagonal and B and C have low ranks
+ z <- 1/A*b
+ E <- diag(nrow(C))+C%*%sweep(B,1,A,"/")#"
+ w <- solve(E,C%*%z)
+ x <- z-sweep(B%*%w,1,A,"/")
+ return( x )
+}
+
+norm <- function(x){
+ sqrt(sum(x^2))
+}
diff --git a/R/parse_intex.R b/R/parse_intex.R new file mode 100644 index 00000000..b9c52e07 --- /dev/null +++ b/R/parse_intex.R @@ -0,0 +1,18 @@ +root = "//WDSENTINEL/share/CorpCDOs"
+source(file.path(root, "R", "intex_deals_functions.R"))
+source(file.path(root, "R", "etdb.R"))
+dealnames <- c("limes", "stonln1")
+cusips <- cusipsfromdealnames(dealnames)
+
+deals.universe <- dbGetQuery(dbCon, "select distinct dealname from clo_universe order by dealname asc")$dealname
+cusips.universe <- cusipsfromdealnames(deals.universe)
+n.scenarios <- 100
+offset <- 2
+r <- data.frame()
+for(cusip in cusips){
+ data <- read.table(paste(cusip, "-PY.txt", sep=""), sep="\t", header=T, nrow=3)
+ price <- sum(as.numeric(sub("\\((.*)\\)", "-\\1", data[1,1:n.scenarios+offset])), na.rm=T)/n.scenarios
+ wal <- sum(as.numeric(sub("\\((.*)\\)", "-\\1", data[2,1:n.scenarios+offset])), na.rm=T)/n.scenarios
+ duration <- sum(as.numeric(sub("\\((.*)\\)", "-\\1", data[3,1:n.scenarios+offset])), na.rm=T)/n.scenarios
+ r <- rbind(r, data.frame(cusip, price, wal, duration))
+}
diff --git a/R/plot_distributions.R b/R/plot_distributions.R new file mode 100644 index 00000000..0b369065 --- /dev/null +++ b/R/plot_distributions.R @@ -0,0 +1,35 @@ +## some plots
+matplot(Z, cbind(w, w.mod), type="l", ylab="probability density",
+ main="Factor distribution (Gaussian, and Market implied)")
+lossdist <- array(0, dim=c(Ngrid, n.int))
+for(i in 1:n.int){
+ pshocked <- sapply(p[,ncol(p)], shockprob, rho=rho, Z=Z[i])
+ S <- 1 - Rstoch[i,,ncol(p)]
+ lossdist[,i] <- lossrecovdist(pshocked, 0, issuerweights, S, Ngrid)$L
+}
+
+lossdist.orig <- BClossdistC(SurvProb, issuerweights, recov, rho, Ngrid)
+matplot(seq(0,1,0.01), cbind(lossdist.orig$L[,ncol(p)], lossdist%*%w.mod), type="l", xlab="loss percentage",
+ ylab="probability density", main="market implied loss distribution")
+#3d surface of the loss distribution
+Rstoch <- array(0, dim=c(ncol(SurvProb), n.int, n.credit))
+for(t in 1:ncol(SurvProb)){
+ for(i in 1:n.credit){
+ Rstoch[t,,i] <- stochasticrecov(recov[i], 0, Z, w.mod, rho, defaultprob[i,t], p[i,t])
+ }
+}
+
+clusterExport(cl, list("p", "shockprob", "rho", "Z", "lossdistribC.joint", "Rstoch", "lu"))
+
+parf <- function(i){
+ pshocked <- apply(p, 2, shockprob, rho=rho, Z=Z[i])
+ return( lossdistribC.joint(pshocked[,ncol(p)], issuerweights, 1-Rstoch[ncol(p),i,], lu) )
+}
+
+dist <- parSapply(cl, 1:n.int, parf)
+dist <- array(dist, dim=c(1/lu+1, 1/lu+1, 100))
+distw <- array(0, dim=c(1/lu+1, 1/lu+1))
+for(i in 1:n.int){
+ distw <- distw+dist[,,i] * w.mod[i]
+}
+persp(distw, theta=-20, phi=13, col="lightgreen",ticktype="detailed", shade=0.6, expand=0.8, border=NA, ltheta=10)
diff --git a/R/reinvestingdistribution.R b/R/reinvestingdistribution.R new file mode 100644 index 00000000..e2c42e3e --- /dev/null +++ b/R/reinvestingdistribution.R @@ -0,0 +1,56 @@ +source("tranche_functions.R")
+
+p <- runif(100)
+S <- rep(0.6, 100)
+lu <- 1
+
+Q <- lossdistribC.joint(p, S, lu)
+#col sums is the recovery distribution
+#row sums is the loss distribution
+
+## some checks
+crossprod(seq(0,1,lu), colSums(Q))-crossprod(p,S)
+crossprod(seq(0,1,lu), rowSums(Q))-crossprod(p,1-S)
+
+## compute the distribution of a reinvested portfolio
+## beta is 1 over the average price at which the loans are bought
+## so higher beta, means more par building
+betavec <- c(0.9,1,1.1,1.5)
+resultx <- c()
+resulty <- c()
+for(beta in betavec){
+ support2 <- outer(1-seq(0,1,0.01), (1-seq(0,1,0.01))^beta,"/")
+ resultx <- cbind(resultx, density(x=support2[-length(support2)],weights=Q[-length(Q)],from=0,to=5)$x)
+ resulty <- cbind(resulty, density(x=support2[-length(support2)],weights=Q[-length(Q)],from=0,to=5)$y)
+}
+
+
+n.int <- 100
+quadrature <- gauss.quad.prob(n.int, "normal")
+w <- quadrature$weights
+Z <- quadrature$nodes
+rho <- 0.45
+Qvec <- array(0, c(100,101,101))
+for(i in 1:n.int){
+ pshocked <- shockprob(p, rho, Z[i])
+ Sshocked <- shockseverity(S, 1, Z[i], rho, p)
+ Qvec[i,,] <- lossdistribC.joint(pshocked, Sshocked, lu)
+}
+
+## plot the density
+for(i in 40:65){
+ persp(z=Qvec[i,,], xlab="recovery", ylab="loss", zlab="density")
+ Sys.sleep(0.3)
+}
+
+## gaussian copula joint distribution of the loss and recovery
+Qresult <- matrix(0, 101, 101)
+for(i in 1:100){
+ Qresult <- Qresult+w[i] * Qvec[i,,]
+}
+
+## 3d tour of the distribution
+for(theta in seq(0, 360, 10)){
+ persp(z=Qresult, xlab="Recovery", ylab="Loss", zlab="density", theta=theta)
+ Sys.sleep(0.5)
+}
diff --git a/R/test.cds_functions.R b/R/test.cds_functions.R new file mode 100644 index 00000000..c30651ba --- /dev/null +++ b/R/test.cds_functions.R @@ -0,0 +1,46 @@ +library("RUnit")
+root = "//WDSENTINEL/share/CorpCDOs/R"
+source(file.path(root, "yieldCurve.R"))
+source(file.path(root, "cds_functions.R"))
+
+#unit tests
+cs <- couponSchedule(as.Date("2012-09-28"), as.Date("2016-04-22"),"Q", "FLOAT", 0.075,0.0703)
+#test that the curve version gives the same result as the flat version for flat curves
+h <- rep(0.05, nrow(cs))
+dpc <- new("defaultprepaycurve", dates=cs$dates, hazardrates=h)
+eps <- 1e-6
+test.flat <- function() {
+ checkEquals(contingentleg.flat(cs, 0.05, 0.7), contingentleg(cs, dpc, 0.7))
+}
+#test derivative of the flat version
+test.derivativeflat <- function() {
+ checkEquals((contingentleg.flat(cs, 0.05+eps, 0.7)-contingentleg.flat(cs, 0.05-eps,0.7))/(2*eps),
+ dcontingentleg.flat(cs, 0.05, 0.7))
+}
+#test derivative of the curved version
+hvec <- c(rep(0.05,5), rep(0.07, 10))
+hvecplus <- c(rep(0.05,5), rep(0.07 + eps,10))
+hvecminus <- c(rep(0.05, 5), rep(0.07 - eps,10))
+index <- c(rep(0,5), rep(1,10))
+dpc <- new("defaultprepaycurve", dates=cs$dates, hazardrates=hvec)
+dpcplus <- new("defaultprepaycurve", dates=cs$dates, hazardrates=hvecplus)
+dpcminus <- new("defaultprepaycurve", dates=cs$dates, hazardrates=hvecminus)
+
+test.derivativecurve <- function() {
+ checkEquals((contingentleg(cs, dpcplus, 0.7) - contingentleg(cs, dpcminus, 0.7))/(2*eps),
+ dcontingentleg(cs, dpc, 0.7, index))
+}
+## test that the prepay version with 0 prepay is the same as the curved version
+dpc <- new("defaultprepaycurve", dates=cs$dates, hazardrates=hvec, prepayrates=rep(0,nrow(cs)))
+test.prepay <- function() {
+ checkEquals(contingentleg(cs, dpc, 0.7),
+ contingentleg.prepay(cs, dpc, 0.7))
+}
+## test derivative of prepay version
+dpcplus <- new("defaultprepaycurve", dates=cs$dates, hazardrates=hvecplus, prepayrates=k(hvecplus))
+dpcminus <- new("defaultprepaycurve", dates=cs$dates, hazardrates=hvecminus, prepayrates=k(hvecminus))
+dpc <- new("defaultprepaycurve", dates=cs$dates, hazardrates=hvec, prepayrates=k(hvec))
+(contingentleg.prepay(cs, dpcplus, 0.7)-contingentleg.prepay(cs, dpcminus, 0.7))/(2*bps)
+dcontingentleg.prepay(cs, dpc, 0.7, index, 15)
+(couponleg.prepay(cs, dpcplus)-couponleg.prepay(cs, dpcminus, 0.7))/(2*bps)
+dcouponleg.prepay(cs, dpc, index, 15)
diff --git a/R/time_of_default.R b/R/time_of_default.R new file mode 100644 index 00000000..369621ab --- /dev/null +++ b/R/time_of_default.R @@ -0,0 +1,99 @@ +spreads<-scan(file="clipboard")
+z1<-rnorm(100)
+
+rho<-0.4
+l<-c()
+for(i in 1:10000){
+ z<-rho*z1+sqrt(1-rho^2)*rnorm(1)
+ u<-pnorm(z)
+ l<-c(l,sum(u))
+}
+
+l<-c()
+for(i in 1:10000){
+ u<-runif(100)
+ l<-c(l,sum(u)/100)
+}
+t<-seq(0,5,1/4)
+f<-function(lambda,t){
+ 1-exp(-lambda*t)
+}
+A<-outer(spreads/(1-0.35),t,f)
+
+
+rho<-0.4
+L<-c()
+z1<-rnorm(100)
+for(i in 1:10000){
+ z1<-rnorm(100)
+ z<-sqrt(1-rho^2)*z1+rho*rnorm(1)
+ u<-pnorm(z)
+ L<-rbind(L,colSums(A>u)/100*(1-0.35))
+}
+
+cdf<-c()
+for(l in 0:100){
+ cdf<-c(cdf,sum(L40[,21]<=l/100*(1-0.35))/100000)
+}
+
+
+d<-seq(-5,5,1)
+p<-abs(1/10*1/d)
+p[6]<-0.1
+pnormvec<-function(z,d){
+ f<-function(d,z){
+ pnorm(z-d)
+ }
+ outer(z,d,"f")
+}
+
+L<-c()
+Z<-matrix(rnorm(1000000),100,10000)
+pv1<-c()
+pv2<-c()
+pv3<-c()
+pv4<-c()
+for(i in 1:10000){
+ z <- Z[,i]
+ D <- sample(d,1,prob=p)
+ ztilde <- z+D
+ u <- pnormvec(z,d)%*%p
+ L <- rbind(L,colSums(A>as.vector(u))/100*(1-0.35))
+ pv1<-c(pv1,pl(t,L[i,],r,0,0.08)-0.05*cl(t,L[i,],r,0,0.08))
+ pv2<-c(pv2,pl(t,L[i,],r,0.08,0.15)-0.05*cl(t,L[i,],r,0.08,0.15))
+ pv3<-c(pv3,pl(t,L[i,],r,0.15,0.3)-0.05*cl(t,L[i,],r,0.15,0.3))
+ pv4<-c(pv4,pl(t,L[i,],r,0.3,1)-0.05*cl(t,L[i,],r,0.3,1))
+}
+pv<-100-100*c(mean(pv1/0.08),mean(pv2/0.07),mean(pv3/0.15),mean(pv4/0.7))
+
+u<-c()
+vecD<-c()
+for(i in 1:1000000){
+ D <- sample(d,1,prob=p)
+ #z<-rnorm(1)
+ #u <- c(u,pnormvec(z+D,d)%*%p)
+ vecD<-c(vecD,D)
+}
+
+tranchesize<-function(L, l, u){
+ if( u==1){
+ u - l - pmin(L,u-l) -pmin( pmax(L-l,0), u-l)
+ }else{
+ u - l - pmin( 0, u - l) - pmin( pmax(L-l,0), u-l)
+ }
+}
+
+protectionsize<-function(L, l, u){
+ pmin( pmax(L-l,0), u-l)
+}
+
+cl<-function(t, L, r, l, u){
+ dT <- diff(t)
+ e <- tranchesize(L,l,u)
+ e <- (e[-length(e)]+e[-1])/2
+ return( sum(dT*exp(-r[-1]*t[-1])*e) )
+}
+pl<-function(t, L, r, l, u){
+ result<-sum(exp(-r[-1]*t[-1])*diff(protectionsize(L,l,u)))
+ return(result)
+}
diff --git a/R/tranche_functions.R b/R/tranche_functions.R new file mode 100644 index 00000000..e1f0ab91 --- /dev/null +++ b/R/tranche_functions.R @@ -0,0 +1,798 @@ +library(statmod)
+
+## todo:
+## -investigate other ways to interpolate the random severities on the grid
+## I'm thinking that at eah severity that we add to the distribution, round it down
+## and keep track of the missing mass: namely if X_i=S_i w.p p_i, then add
+## X_i=lu*floor(S_i/lu) with probability p_i and propagate
+## X_{i+1}=S_{i+1}+(S_i-lu*floor(S_i/lu)) with the right probability
+## - investigate truncated distributions more (need to compute loss and recov distribution
+## separately, for the 0-10 equity tranche, we need the loss on the 0-0.1 support and
+## recovery with 0.1-1 support, so it's not clear that there is a big gain.
+## - do the correlation adjustments when computing the deltas since it seems to be
+## the market standard
+
+lossdistrib <- function(p){
+ ## basic recursive algorithm of Andersen, Sidenius and Basu
+ n <- length(p)
+ q <- rep(0, (n+1))
+ q[1] <- 1
+ for(i in 1:n){
+ q[-1] <- p[i]*q[-(n+1)]+(1-p[i])*q[-1]
+ q[1] <- (1-p[i])*q[1]
+ }
+ return(q)
+}
+
+lossdistrib.fft <- function(p){
+ ## computes loss distribution using the fft
+ ## complexity is of order O(n*m)+O(m*log(m))
+ ## where m is the size of the grid and n the number of probabilities.
+ ## this is slower than the recursive algorithm
+ theta <- 2*pi*1i*(0:n)/(n+1)
+ Phi <- 1 - p + p%o%exp(theta)
+ v <- apply(Phi, 2, prod)
+ return(1/(n+1)*Re(fft(v)))
+}
+
+lossdistrib2 <- function(p, w, S, N, defaultflag=FALSE){
+ ## recursive algorithm with first order correction
+ ## p vector of default probabilities
+ ## w vector of weigths
+ ## S vector of severities
+ ## N number of ticks in the grid
+ ## defaultflag if true computes the default distribution
+ n <- length(p)
+ lu <- 1/(N-1)
+ q <- rep(0, N)
+ q[1] <- 1
+ for(i in 1:n){
+ if(defaultflag){
+ d <- w[i] /lu
+ }else{
+ d <- S[i] * w[i] / lu
+ }
+ d1 <- floor(d)
+ d2 <- ceiling(d)
+ p1 <- p[i]*(d2-d)
+ p2 <- p[i] - p1
+ q1 <- c(rep(0,d1), p1*q[1:(N-d1)])
+ q2 <- c(rep(0,d2), p2*q[1:(N-d2)])
+ q <- q1 + q2 + (1-p[i])*q
+ }
+ q[length(q)] <- q[length(q)]+1-sum(q)
+ return(q)
+}
+
+lossdistrib2.truncated <- function(p, w, S, N, cutoff=N){
+ ## recursive algorithm with first order correction
+ ## p vector of default probabilities
+ ## w vector of weigths
+ ## S vector of severities
+ ## N number of ticks in the grid (for best accuracy should
+ ## be a multiple of the number of issuers)
+ ## cutoff where to stop computing the exact probabilities
+ ## (useful for tranche computations)
+
+ ## this is actually slower than lossdistrib2. But in C this is
+ ## twice as fast.
+ ## for high severities, M can become bigger than N, and there is
+ ## some probability mass escaping.
+ n <- length(p)
+ lu <- 1/(N-1)
+ q <- rep(0, truncated)
+ q[1] <- 1
+ M <- 1
+ for(i in 1:n){
+ d <- S[i] * w[i] / lu
+ d1 <- floor(d)
+ d2 <- ceiling(d)
+ p1 <- p[i]*(d2-d)
+ p2 <- p[i] - p1
+ q1 <- p1*q[1:min(M, cutoff-d1)]
+ q2 <- p2*q[1:min(M, cutoff-d2)]
+ q[1:min(M, cutoff)] <- (1-p[i])*q[1:min(M, cutoff)]
+ q[(d1+1):min(M+d1, cutoff)] <- q[(d1+1):min(M+d1, cutoff)]+q1
+ q[(d2+1):min(M+d2, cutoff)] <- q[(d2+1):min(M+d2, cutoff)]+q2
+ M <- M+d2
+ }
+ return(q)
+}
+
+recovdist <- function(dp, pp, w, S, N){
+ ## computes the recovery distribution for a sum of independent variables
+ ## R=\sum_{i=1}^n X_i
+ ## where X_i = 0 w.p 1-dp[i]-pp[i]
+ ## = w[i]*(1-S[i]) w.p dp[i]
+ ## = w[i] w.p pp[i]
+ ## each non zero value v is interpolated on the grid as
+ ## the pair of values floor(v/lu) and ceiling(v/lu) so that
+ ## X_i has four non zero values
+ n <- length(dp)
+ q <- rep(0, N)
+ q[1] <- 1
+ lu <- 1/(N-1)
+ for(i in 1:n){
+ d1 <- w[i]*(1-S[i])/lu
+ d1l <- floor(d1)
+ d1u <- ceiling(d1)
+ d2 <- w[i] / lu
+ d2l <- floor(d2)
+ d2u <- ceiling(d2)
+ dp1 <- dp[i] * (d1u-d1)
+ dp2 <- dp[i] - dp1
+ pp1 <- pp[i] * (d2u - d2)
+ pp2 <- pp[i] - pp1
+ q1 <- c(rep(0, d1l), dp1 * q[1:(N-d1l)])
+ q2 <- c(rep(0, d1u), dp2 * q[1:(N-d1u)])
+ q3 <- c(rep(0, d2l), pp1 * q[1:(N-d2l)])
+ q4 <- c(rep(0, d2u), pp2 *q[1:(N-d2u)])
+ q <- q1+q2+q3+q4+(1-dp[i]-pp[i])*q
+ }
+ return(q)
+}
+
+lossdist.joint <- function(p, w, S, N, defaultflag=FALSE){
+ ## recursive algorithm with first order correction
+ ## to compute the joint probability distribution of the loss and recovery
+ ## inputs:
+ ## p: vector of default probabilities
+ ## w: vector of issuer weights
+ ## S: vector of severities
+ ## N: number of tick sizes on the grid
+ ## defaultflag: if true computes the default distribution
+ ## output:
+ ## q: matrix of joint loss, recovery probability
+ ## colSums(q) is the recovery distribution marginal
+ ## rowSums(q) is the loss distribution marginal
+ n <- length(p)
+ lu <- 1/(N-1)
+ q <- matrix(0, N, N)
+ q[1,1] <- 1
+ for(k in 1:n){
+ if(defaultflag){
+ x <- w[k] / lu
+ }else{
+ x <- S[k] * w[k]/lu
+ }
+ y <- (1-S[k]) * w[k]/lu
+ i <- floor(x)
+ j <- floor(y)
+ weights <- c((i+1-x)*(j+1-y), (i+1-x)*(y-j), (x-i)*(y-j), (j+1-y)*(x-i))
+ psplit <- p[k] * weights
+ qtemp <- matrix(0, N, N)
+ qtemp[(i+1):N,(j+1):N] <- qtemp[(i+1):N,(j+1):N] + psplit[1] * q[1:(N-i),1:(N-j)]
+ qtemp[(i+1):N,(j+2):N] <- qtemp[(i+1):N,(j+2):N] + psplit[2] * q[1:(N-i), 1:(N-j-1)]
+ qtemp[(i+2):N,(j+2):N] <- qtemp[(i+2):N,(j+2):N] + psplit[3] * q[1:(N-i-1), 1:(N-j-1)]
+ qtemp[(i+2):N, (j+1):N] <- qtemp[(i+2):N, (j+1):N] + psplit[4] * q[1:(N-i-1), 1:(N-j)]
+ q <- qtemp + (1-p[k])*q
+ }
+ q[length(q)] <- q[length(q)]+1-sum(q)
+ return(q)
+}
+
+lossdist.prepay.joint <- function(dp, pp, w, S, N, defaultflag=FALSE){
+ ## recursive algorithm with first order correction
+ ## to compute the joint probability distribition of the loss and recovery
+ ## inputs:
+ ## dp: vector of default probabilities
+ ## pp: vector of prepay probabilities
+ ## w: vector of issuer weights
+ ## S: vector of severities
+ ## N: number of tick sizes on the grid
+ ## defaultflag: if true computes the default
+ ## outputs
+ ## q: matrix of joint loss and recovery probability
+ ## colSums(q) is the recovery distribution marginal
+ ## rowSums(q) is the loss distribution marginal
+ n <- length(dp)
+ lu <- 1/(N-1)
+ q <- matrix(0, N, N)
+ q[1,1] <- 1
+ for(k in 1:n){
+ y1 <- (1-S[k]) * w[k]/lu
+ y2 <- w[k]/lu
+ j1 <- floor(y1)
+ j2 <- floor(y2)
+ if(defaultflag){
+ x <- y2
+ i <- j2
+ }else{
+ x <- y2-y1
+ i <- floor(x)
+ }
+
+ ## weights <- c((i+1-x)*(j1+1-y1), (i+1-x)*(y1-j1), (x-i)*(y1-j1), (j1+1-y1)*(x-i))
+ weights1 <- c((i+1-x)*(j1+1-y1), (i+1-x)*(y1-j1), (x-i)*(y1-j1), (j1+1-y1)*(x-i))
+ dpsplit <- dp[k] * weights1
+
+ if(defaultflag){
+ weights2 <- c((i+1-x)*(j2+1-y2), (i+1-x)*(y2-j2), (x-i)*(y2-j2), (j2+1-y2)*(x-i))
+ ppsplit <- pp[k] * weights2
+ }else{
+ ppsplit <- pp[k] * c(j2+1-y2, y2-j2)
+ }
+ qtemp <- matrix(0, N, N)
+ qtemp[(i+1):N,(j1+1):N] <- qtemp[(i+1):N,(j1+1):N] + dpsplit[1] * q[1:(N-i),1:(N-j1)]
+ qtemp[(i+1):N,(j1+2):N] <- qtemp[(i+1):N,(j1+2):N] + dpsplit[2] * q[1:(N-i), 1:(N-j1-1)]
+ qtemp[(i+2):N,(j1+2):N] <- qtemp[(i+2):N,(j1+2):N] + dpsplit[3] * q[1:(N-i-1), 1:(N-j1-1)]
+ qtemp[(i+2):N,(j1+1):N] <- qtemp[(i+2):N, (j1+1):N] + dpsplit[4] * q[1:(N-i-1), 1:(N-j1)]
+ if(defaultflag){
+ qtemp[(i+1):N,(j2+1):N] <- qtemp[(i+1):N,(j2+1):N] + ppsplit[1] * q[1:(N-i),1:(N-j2)]
+ qtemp[(i+1):N,(j2+2):N] <- qtemp[(i+1):N,(j2+2):N] + ppsplit[2] * q[1:(N-i), 1:(N-j2-1)]
+ qtemp[(i+2):N,(j2+2):N] <- qtemp[(i+2):N,(j2+2):N] + ppsplit[3] * q[1:(N-i-1), 1:(N-j2-1)]
+ qtemp[(i+2):N,(j2+1):N] <- qtemp[(i+2):N, (j2+1):N] + ppsplit[4] * q[1:(N-i-1), 1:(N-j2)]
+ }else{
+ qtemp[, (j2+1):N] <- qtemp[,(j2+1):N]+ppsplit[1]*q[,1:(N-j2)]
+ qtemp[, (j2+2):N] <- qtemp[,(j2+2):N]+ppsplit[2]*q[,1:(N-j2-1)]
+ }
+ q <- qtemp + (1-pp[k]-dp[k]) * q
+ }
+ q[length(q)] <- q[length(q)] + 1 - sum(q)
+ return(q)
+}
+
+lossdistC <- function(p, w, S, N, defaultflag=FALSE){
+ ## C version of lossdistrib2, roughly 50 times faster
+ if(!is.loaded("lossdistrib")){
+ dyn.load(paste0("lossdistrib", .Platform$dynlib.ext))
+ }
+ .C("lossdistrib", as.double(p), as.integer(length(p)),
+ as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), q = double(N))$q
+}
+
+lossdistC.truncated <- function(p, w, S, N, T=N){
+ ## C version of lossdistrib2, roughly 50 times faster
+ if(!is.loaded("lossdistrib_truncated")){
+ dyn.load(paste0("lossdistrib", .Platform$dynlib.ext))
+ }
+ dyn.load("lossdistrib.dll")
+ .C("lossdistrib_truncated", as.double(p), as.integer(length(p)),
+ as.double(w), as.double(S), as.integer(N), as.integer(T), q = double(T))$q
+}
+
+recovdistC <- function(dp, pp, w, S, N){
+ ## C version of recovdist
+ if(!is.loaded("recovdist")){
+ dyn.load(paste0("lossdistrib", .Platform$dynlib.ext))
+ }
+ .C("recovdist", as.double(dp), as.double(pp), as.integer(length(dp)),
+ as.double(w), as.double(S), as.integer(N), q = double(N))$q
+}
+
+lossdistC.joint <- function(p, w, S, N, defaultflag=FALSE){
+ ## C version of lossdistrib.joint, roughly 20 times faster
+ if(!is.loaded("lossdistrib_joint")){
+ dyn.load(paste0("lossdistrib", .Platform$dynlib.ext))
+ }
+ .C("lossdistrib_joint", as.double(p), as.integer(length(p)), as.double(w),
+ as.double(S), as.integer(N), as.logical(defaultflag), q = matrix(0, N, N))$q
+}
+
+lossdistC.prepay.joint <- function(dp, pp, w, S, N, defaultflag=FALSE){
+ ## C version of lossdist.prepay.joint
+ if(!is.loaded("lossdistrib_prepay_joint")){
+ dyn.load(paste0("lossdistrib", .Platform$dynlib.ext))
+ }
+ r <- .C("lossdistrib_prepay_joint", as.double(dp), as.double(pp), as.integer(length(dp)),
+ as.double(w), as.double(S), as.integer(N), as.logical(defaultflag), q=matrix(0, N, N))$q
+ gc()
+ return(r)
+}
+
+lossrecovdist <- function(defaultprob, prepayprob, w, S, N, defaultflag=FALSE, useC=TRUE){
+ if(missing(prepayprob)){
+ if(useC){
+ L <- lossdistC(defaultprob, w, S, N, defaultflag)
+ R <- lossdistC(defaultprob, w, 1-S, N)
+ }else{
+ L <- lossdistrib2(defaultprob, w, S, N, defaultflag)
+ R <- lossdistrib2(defaultprob, w, 1-S, N)
+ }
+ }else{
+ if(useC){
+ L <- lossdistC(defaultprob+defaultflag*prepayprob, w, S, N, defaultflag)
+ R <- recovdistC(defaultprob, prepayprob, w, S, N)
+ }else{
+ L <- lossdistrib2(defaultprob+defaultflag*prepayprob, w, S, N, defaultflag)
+ R <- recovdist(defaultprob, prepayprob, w, S, N)
+ }
+ }
+ return(list(L=L, R=R))
+}
+
+lossrecovdist.term <- function(defaultprob, prepayprob, w, S, N, defaultflag=FALSE, useC=TRUE){
+ ## computes the loss and recovery distribution over time
+ L <- array(0, dim=c(N, ncol(defaultprob)))
+ R <- array(0, dim=c(N, ncol(defaultprob)))
+ if(missing(prepayprob)){
+ for(t in 1:ncol(defaultprob)){
+ temp <- lossrecovdist(defaultprob[,t], , w, S[,t], N, defaultflag, useC)
+ L[,t] <- temp$L
+ R[,t] <- temp$R
+ }
+ }else{
+ for(t in 1:ncol(defaultprob)){
+ temp <- lossrecovdist(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag, useC)
+ L[,t] <- temp$L
+ R[,t] <- temp$R
+ }
+ }
+ return(list(L=L, R=R))
+}
+
+lossrecovdist.joint.term <- function(defaultprob, prepayprob, w, S, N, defaultflag=FALSE, useC=TRUE){
+ ## computes the joint loss and recovery distribution over time
+ Q <- array(0, dim=c(ncol(defaultprob), N, N))
+ if(useC){
+ if(missing(prepayprob)){
+ for(t in 1:ncol(defaultprob)){
+ Q[t,,] <- lossdistC.joint(defaultprob[,t], w, S[,t], N, defaultflag)
+ }
+ }else{
+ for(t in 1:ncol(defaultprob)){
+ Q[t,,] <- lossdistC.prepay.joint(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag)
+ }
+ }
+ }else{
+ if(missing(prepayprob)){
+ for(t in 1:ncol(defaultprob)){
+ Q[t,,] <- lossdist.joint(defaultprob[,t], w, S[,t], N, defaultflag)
+ }
+ }else{
+ for(t in 1:ncol(defaultprob)){
+ Q[t,,] <- lossdist.prepay.joint(defaultprob[,t], prepayprob[,t], w, S[,t], N, defaultflag)
+ }
+ }
+ }
+ gc()
+ return(Q)
+}
+
+dist.transform <- function(dist.joint){
+ ## compute the joint (D, R) distribution
+ ## from the (L, R) distribution using D = L+R
+ distDR <- array(0, dim=dim(dist.joint))
+ Ngrid <- dim(dist.joint)[2]
+ u <- seq(0, 1, length=Ngrid)
+ v <- seq(0, 1, length=Ngrid)
+ for(t in 1:ncol(dp)){
+ for(i in 1:Ngrid){
+ for(j in 1:Ngrid){
+ index <- i+j
+ if(index <= Ngrid){
+ distDR[t,index,j] <- distDR[t,index,j] + dist.joint[t,i,j]
+ }
+ }
+ }
+ }
+ return( distDR )
+}
+
+shockprob <- function(p, rho, Z, log.p=F){
+ ## computes the shocked default probability as a function of the copula factor
+ pnorm((qnorm(p)-sqrt(rho)*Z)/sqrt(1-rho), log.p=log.p)
+}
+
+shockseverity <- function(S, Stilde=1, Z, rho, p){
+ #computes the severity as a function of the copula factor Z
+ Stilde * exp( shockprob(S/Stilde*p, rho, Z, TRUE) - shockprob(p, rho, Z, TRUE))
+}
+
+dshockprob <- function(p,rho,Z){
+ dnorm((qnorm(p)-sqrt(rho)*Z)/sqrt(1-rho))*dqnorm(p)/sqrt(1-rho)
+}
+
+dqnorm <- function(x){
+ 1/dnorm(qnorm(x))
+}
+
+fit.prob <- function(Z, w, rho, p0){
+ ## if the weights are not perfectly gaussian, find the probability p such
+ ## E_w(shockprob(p, rho, Z)) = p0
+ if(p0==0){
+ return(0)
+ }else{
+ eps <- 1e-12
+ dp <- (crossprod(shockprob(p0,rho,Z),w)-p0)/crossprod(dshockprob(p0,rho,Z),w)
+ p <- p0
+ while(abs(dp) > eps){
+ dp <- (crossprod(shockprob(p,rho,Z),w)-p0)/crossprod(dshockprob(p,rho,Z),w)
+ phi <- 1
+ while ((p-phi*dp)<0 || (p-phi*dp)>1){
+ phi <- 0.8*phi
+ }
+ p <- p - phi*dp
+ }
+ return(p)
+ }
+}
+
+stochasticrecov <- function(R, Rtilde, Z, w, rho, porig, pmod){
+ ptilde <- fit.prob(Z, w, rho, (1-R)/(1-Rtilde) * porig)
+ return(abs(1-(1-Rtilde) * exp(shockprob(ptilde, rho, Z, TRUE) - shockprob(pmod, rho, Z, TRUE))))
+}
+
+pos <- function(x){
+ pmax(x, 0)
+}
+
+trancheloss <- function(L, K1, K2){
+ pos(L - K1) - pos(L - K2)
+}
+
+trancherecov <- function(R, K1, K2){
+ pos(R - 1 + K2) - pos(R - 1 +K1)
+}
+
+tranche.cl <- function(L, R, cs, K1, K2, Ngrid=nrow(L), scaled=FALSE){
+ ## computes the couponleg of a tranche
+ ## if scaled is TRUE, scale it by the size of the tranche (K2-K1)
+ ## can make use of the fact that the loss and recov distribution are
+ ## truncated (in that case nrow(L) != Ngrid
+ if(K1==K2){
+ return( 0 )
+ }else{
+ support <- seq(0, 1, length=Ngrid)[1:nrow(L)]
+ size <- K2 - K1 - crossprod(trancheloss(support, K1, K2), L) -
+ crossprod(trancherecov(support, K1, K2), R)
+ sizeadj <- as.numeric(0.5 * (size + c(K2-K1, size[-length(size)])))
+ if(scaled){
+ return( 1/(K2-K1) * crossprod(sizeadj * cs$coupons, cs$df) )
+ }else{
+ return( crossprod(sizeadj * cs$coupons, cs$df) )
+ }
+ }
+}
+
+tranche.cl.scenarios <- function(l, r, cs, K1, K2, scaled=FALSE){
+ ## computes the couponleg of a tranche for one scenario
+ ## if scaled is TRUE, scale it by the size of the tranche (K2-K1)
+ ## can make use of the fact that the loss and recov distribution are
+ ## truncated (in that case nrow(L) != Ngrid
+ if(K1==K2){
+ return( 0 )
+ }else{
+ size <- K2 - K1 - trancheloss(l, K1, K2) - trancherecov(r, K1, K2)
+ sizeadj <- as.numeric(0.5 * (size + c(K2-K1, size[-length(size)])))
+ if(scaled){
+ return( 1/(K2-K1) * crossprod(sizeadj * cs$coupons, cs$df) )
+ }else{
+ return( crossprod(sizeadj * cs$coupons, cs$df) )
+ }
+ }
+}
+
+tranche.pl <- function(L, cs, K1, K2, Ngrid=nrow(L), scaled=FALSE){
+ ## computes the protection leg of a tranche
+ ## if scaled
+ if(K1==K2){
+ return(0)
+ }else{
+ support <- seq(0, 1, length=Ngrid)[1:nrow(L)]
+ cf <- K2 - K1 - crossprod(trancheloss(support, K1, K2), L)
+ cf <- c(K2 - K1, cf)
+ if(scaled){
+ return( 1/(K2-K1) * crossprod(diff(cf), cs$df))
+ }else{
+ return( crossprod(diff(cf), cs$df))
+ }
+ }
+}
+
+tranche.pl.scenarios <- function(l, cs, K1, K2, scaled=FALSE){
+ ## computes the protection leg of a tranche
+ ## if scaled
+ if(K1==K2){
+ return(0)
+ }else{
+ cf <- K2 - K1 - trancheloss(l, K1, K2)
+ cf <- c(K2 - K1, cf)
+ if(scaled){
+ return( 1/(K2-K1) * as.numeric(crossprod(diff(cf), cs$df)))
+ }else{
+ return( as.numeric(crossprod(diff(cf), cs$df)))
+ }
+ }
+}
+
+tranche.pv <- function(L, R, cs, K1, K2, Ngrid=nrow(L)){
+ return( tranche.pl(L, cs, K1, K2, Ngrid) + tranche.cl(L, R, cs, K1, K2, Ngrid))
+}
+
+tranche.pv.scenarios <- function(l, r, cs, K1, K2){
+ return( tranche.pl.scenarios(l, cs, K1, K2, TRUE) +
+ tranche.cl.scenarios(l, r, cs, K1, K2, TRUE))
+}
+
+
+adjust.attachments <- function(K, losstodate, factor){
+ ## computes the attachments adjusted for losses
+ ## on current notional
+ return( pmin(pmax((K-losstodate)/factor, 0),1) )
+}
+
+tranche.pvvec <- function(K, L, R, cs){
+ r <- rep(0, length(K)-1)
+ for(i in 1:(length(K)-1)){
+ r[i] <- 1/(K[i+1]-K[i]) * tranche.pv(L, R, cs, K[i], K[i+1])
+ }
+ return( r )
+}
+
+BClossdist <- function(SurvProb, issuerweights, recov, rho, N=length(recov)+1,
+ n.int=100){
+ quadrature <- gauss.quad.prob(n.int, "normal")
+ Z <- quadrature$nodes
+ w <- quadrature$weights
+ LZ <- matrix(0, N, n.int)
+ RZ <- matrix(0, N, n.int)
+ L <- matrix(0, N, ncol(SurvProb))
+ R <- matrix(0, N, ncol(SurvProb))
+ for(t in 1:ncol(SurvProb)){
+ g <- 1 - SurvProb[, t]
+ for(i in 1:length(Z)){
+ g.shocked <- shockprob(g, rho, Z[i])
+ S.shocked <- shockseverity(1-recov, 1, Z[i], rho, g)
+ temp <- lossrecovdist(g.shocked, 0, issuerweights, S.shocked, N)
+ LZ[,i] <- temp$L
+ RZ[,i] <- temp$R
+ }
+ L[,t] <- LZ%*%w
+ R[,t] <- RZ%*%w
+ }
+ list(L=L, R=R)
+}
+
+BClossdistC <- function(SurvProb, issuerweights, recov, rho,
+ N=length(issuerweights)+1, n.int=100, defaultflag){
+ if(!is.loaded("BClossdist")){
+ dyn.load(paste0("lossdistrib", .Platform$dynlib.ext))
+ }
+ quadrature <- gauss.quad.prob(n.int, "normal")
+ Z <- quadrature$nodes
+ w <- quadrature$weights
+ L <- matrix(0, N, dim(SurvProb)[2])
+ R <- matrix(0, N, dim(SurvProb)[2])
+ r <- .C("BClossdist", SurvProb, as.integer(dim(SurvProb)[1]), as.integer(dim(SurvProb)[2]),
+ as.double(issuerweights), as.double(recov), as.double(Z), as.double(w),
+ as.integer(n.int), as.double(rho), as.integer(N), as.logical(defaultflag), L=L, R=R)
+ return(list(L=r$L,R=r$R))
+}
+
+BCtranche.pv <- function(portfolio, index, coupon, K1, K2, rho1, rho2, N=length(portfolio)+1){
+ ## computes the protection leg, couponleg, and bond price of a tranche
+ ## in the base correlation setting
+ if(K1==0){
+ if(rho1!=0){
+ stop("equity tranche must have 0 lower correlation")
+ }
+ }
+ SurvProb <- SPmatrix(portfolio, index)
+ cs <- couponSchedule(nextIMMDate(Sys.Date()), index$maturity, "Q", "FIXED", coupon, 0)
+ recov <- sapply(portfolio, attr, "recovery")
+ issuerweights <- rep(1/length(portfolio), length(portfolio))
+ K <- adjust.attachments(c(K1,K2), index$loss, index$factor)
+ dK <- K[2] - K[1]
+ dist2 <- BClossdistC(SurvProb, issuerweights, recov, rho2, N)
+ if(rho1!=0){
+ dist1 <- BClossdistC(SurvProb, issuerweights, recov, rho1, N)
+ }
+ cl2 <- tranche.cl(dist2$L, dist2$R, cs, 0, K[2])
+ cl1 <- tranche.cl(dist1$L, dist1$R, cs, 0, K[1])
+ pl2 <- tranche.pl(dist2$L, cs, 0, K[2])
+ pl1 <- tranche.pl(dist1$L, cs, 0, K[1])
+ return(list(pl=(pl2-pl1)/dK, cl=(cl2-cl1)/dK,
+ bp=100*(1+(pl2-pl1+cl2-cl1)/dK)))
+}
+
+BCtranche.delta <- function(portfolio, index, coupon, K1, K2, rho1, rho2, N=length(portolio)+1){
+ ## computes the tranche delta (on current notional) by doing a proportional
+ ## blip of all the curves
+ ## if K2==1, then computes the delta using the lower attachment only
+ ## this makes sense for bottom-up skews
+ eps <- 1e-4
+ portfolioplus <- portfolio
+ portfoliominus <- portfolio
+ for(i in 1:length(portfolio)){
+ portfolioplus[[i]]@curve@hazardrates <- portfolioplus[[i]]@curve@hazardrates * (1 + eps)
+ portfoliominus[[i]]@curve@hazardrates <- portfoliominus[[i]]@curve@hazardrates * (1- eps)
+ }
+ dPVindex <- indexpv(portfolioplus, index)$bp - indexpv(portfoliominus, index)$bp
+ if(K2==1){
+ dPVtranche <- BCtranche.pv(portfolioplus, index, coupon, 0, K1, 0, rho1, N)$bp -
+ BCtranche.pv(portfoliominus, index, coupon, 0, K1, 0, rho1, N)$bp
+ K1adj <- adjust.attachments(K1, index$loss, index$factor)
+ delta <- (1 - dPVtranche/(100*dPVindex) * K1adj)/(1-K1adj)
+ }else{
+ dPVtranche <- BCtranche.pv(portfolioplus, index, coupon, K1, K2, rho1, rho2, N)$bp -
+ BCtranche.pv(portfoliominus, index, coupon, K1, K2, rho1, rho2, N)$bp
+ delta <- dPVtranche/(100*dPVindex)
+ }
+ ## dPVindex <- BCtranche.pv(portfolioplus, index, coupon, 0, 1, 0, 0.5, lu)$bp-
+ ## BCtranche.pv(portfoliominus, index, coupon, 0, 1, 0, 0.5, lu)$bp
+ return( delta )
+}
+
+BCstrikes <- function(portfolio, index, coupon, K, rho, N=101) {
+ ## computes the strikes as a percentage of expected loss
+ EL <- c()
+ for(i in 2:length(K)){
+ EL <- c(EL, -BCtranche.pv(portfolio, index, coupon, K[i-1], K[i], rho[i-1], rho[i], N)$pl)
+ }
+ Kmodified <- adjust.attachments(K, index$loss, index$factor)
+ return(cumsum(EL*diff(Kmodified))/sum(EL*diff(Kmodified)))
+}
+
+delta.factor <- function(K1, K2, index){
+ ## compute the factor to convert from delta on current notional to delta on original notional
+ ## K1 and K2 original strikes
+ factor <- (adjust.attachments(K2, index$loss, index$factor)
+ -adjust.attachments(K1, index$loss, index$factor))/(K2-K1)
+ return( factor )
+}
+
+MFupdate.prob <- function(Z, w, rho, defaultprob){
+ ## update the probabilites based on a non gaussian factor
+ ## distribution so that the pv of the cds stays the same.
+ p <- matrix(0, nrow(defaultprob), ncol(defaultprob))
+ for(i in 1:nrow(defaultprob)){
+ for(j in 1:ncol(defaultprob)){
+ p[i,j] <- fit.prob(Z, w, rho, defaultprob[i,j])
+ }
+ }
+ return( p )
+}
+
+MFlossrecovdist <- function(w, Z, rho, defaultprob, defaultprobmod, issuerweights, recov,
+ Ngrid=2*length(issuerweights)+1, defaultflag=FALSE){
+ ## computes the loss and recovery distribution using the modified factor distribution
+ n.credit <- length(issuerweights)
+ n.int <- length(w)
+ Rstoch <- array(0, dim=c(n.int, n.credit, ncol(defaultprob)))
+ for(t in 1:ncol(defaultprob)){
+ for(i in 1:n.credit){
+ Rstoch[,i,t] <- stochasticrecov(recov[i], 0, Z, w, rho, defaultprob[i,t], defaultprobmod[i,t])
+ }
+ }
+ parf <- function(i){
+ pshocked <- apply(defaultprobmod, 2, shockprob, rho=rho, Z=Z[i])
+ S <- 1 - Rstoch[i,,]
+ dist <- lossrecovdist.term(pshocked, 0, issuerweights, S, Ngrid, defaultflag)
+ }
+ L <- matrix(0, Ngrid, ncol(defaultprob))
+ R <- matrix(0, Ngrid, ncol(defaultprob))
+ for(i in 1:length(w)){
+ dist <- parf(i)
+ L <- L + dist$L * w[i]
+ R <- R + dist$R * w[i]
+ }
+ return( list(L=L, R=R) )
+}
+
+MFlossrecovdist.prepay <- function(w, Z, rho, defaultprob, defaultprobmod, prepayprob, prepayprobmod,
+ issuerweights, recov, Ngrid=2*length(issuerweights)+1, defaultflag=FALSE){
+ ## computes the loss and recovery distribution using the modified factor distribution
+ n.credit <- length(issuerweights)
+ n.int <- length(w)
+ Rstoch <- array(0, dim=c(n.int, n.credit, ncol(defaultprob)))
+ for(t in 1:ncol(defaultprob)){
+ for(i in 1:n.credit){
+ Rstoch[,i,t] <- stochasticrecov(recov[i], 0, Z, w, rho, defaultprob[i,t], defaultprobmod[i,t])
+ }
+ }
+ parf <- function(i){
+ dpshocked <- apply(defaultprobmod, 2, shockprob, rho=rho, Z=Z[i])
+ ppshocked <- apply(prepayprobmod, 2, shockprob, rho=rho, Z=-Z[i])
+ S <- 1 - Rstoch[i,,]
+ dist <- lossrecovdist.term(dpshocked, ppshocked, issuerweights, S, Ngrid, defaultflag)
+ }
+ L <- matrix(0, Ngrid, ncol(defaultprob))
+ R <- matrix(0, Ngrid, ncol(defaultprob))
+ for(i in 1:length(w)){
+ dist <- parf(i)
+ L <- L + dist$L * w[i]
+ R <- R + dist$R * w[i]
+ }
+ return( list(L=L, R=R) )
+}
+
+MFlossdist.joint <- function(cl, w, Z, rho, defaultprob, defaultprobmod, issuerweights, recov,
+ Ngrid=2*length(issuerweights)+1, defaultflag=FALSE){
+ ## rowSums(Q) is the loss/default distribution
+ ## colSums(Q) is the recovery distribution
+ ## so that recovery is the y axis and L/D is the x axis
+ ## if we use the persp function, losses is the axes facing us,
+ ## and R is the axis going away from us.
+ n.credit <- length(issuerweights)
+ n.int <- lenth(w)
+ Rstoch <- array(0, dim=c(n.int, n.credit, ncol(defaultprob)))
+ for(t in 1:ncol(defaultprob)){
+ for(i in 1:n.credit){
+ Rstoch[,i,t] <- stochasticrecov(recov[i], 0, Z, w, rho, defaultprob[i,t], defaultprobmod[i,t])
+ }
+ }
+ parf <- function(i){
+ pshocked <- apply(defaultprobmod, 2, shockprob, rho=rho, Z=Z[i])
+ S <- 1 - Rstoch[i,,]
+ dist <- lossrecovdist.joint.term(pshocked, 0, issuerweights, S, Ngrid, defaultflag)
+ gc()
+ return(dist)
+ }
+ temp <- parSapply(cl, 1:length(w), parf)
+ clusterCall(cl, gc)
+ Q <- array(0, dim=c(ncol(defaultprob), Ngrid, Ngrid))
+ for(i in 1:length(w)){
+ Q <- Q + w[i]*array(temp[,i], dim=c(ncol(defaultprob), Ngrid, Ngrid))
+ }
+ return( Q )
+}
+
+MFlossdist.prepay.joint <- function(cl, w, Z, rho, defaultprob, defaultprobmod,
+ prepayprob, prepayprobmod, issuerweights, recov,
+ Ngrid=2*length(issuerweights)+1, defaultflag=FALSE, n.chunks=4){
+ ## rowSums is the loss distribution
+ ## colSums is the recovery distribution
+ ## so that recovery is the y axis and L is the x axis
+ ## if we use the persp function, losses is the axes facing us,
+ ## and R is the axis going away from us.
+ n.credit <- length(issuerweights)
+ n.int <- length(w)
+ Rstoch <- array(0, dim=c(n.int, n.credit, ncol(defaultprob)))
+ for(t in 1:ncol(defaultprob)){
+ for(i in 1:n.credit){
+ Rstoch[,i,t] <- stochasticrecov(recov[i], 0, Z, w, rho, defaultprob[i,t], defaultprobmod[i,t])
+ }
+ }
+ parf <- function(i){
+ dpshocked <- apply(defaultprobmod, 2, shockprob, rho=rho, Z=Z[i])
+ ppshocked <- apply(prepayprobmod, 2, shockprob, rho=rho, Z=-Z[i])
+ S <- 1 - Rstoch[i,,]
+ dist <- lossrecovdist.joint.term(dpshocked, ppshocked, issuerweights, S, Ngrid, defaultflag)
+ return(dist)
+ }
+ Q <- array(0, dim=c(ncol(defaultprob), Ngrid, Ngrid))
+ ## when using the cluster this takes a crazy amount of memory
+ ## we cut it in chunks to keep it manageable
+ chunksize <- length(w)/n.chunks
+ for(chunk in 1:n.chunks){
+ if(length(cl)>0){
+ clusterExport(cl, list("Rstoch", "defaultprobmod", "prepayprobmod"), envir=environment())
+ temp <- parSapply(cl, (1+chunksize*(chunk-1)):(chunksize*chunk), parf)
+ ## need to call gc twice for some reason, otherwise doesn't release everything
+ clusterCall(cl, gc)
+ clusterCall(cl, gc)
+ }else{
+ temp <- sapply((1+chunksize*(chunk-1)):chunksize*chunk, parf)
+ }
+ for(i in 1:chunksize){
+ Q <- Q + w[i+chunksize*(chunk-1)]*array(temp[,i], dim=c(ncol(defaultprob), Ngrid, Ngrid))
+ }
+ ## same here
+ gc()
+ gc()
+ }
+ return( Q )
+}
+
+MFtranche.pv <- function(cl, cs, w, rho, defaultprob, defaultprobmod, issuerweights, recov,
+ Kmodified, Ngrid=length(issuerweights)+1){
+ ## computes the tranches pv using the modified factor distribution
+ ## p is the modified probability so that
+ n.credit <- length(issuerweights)
+ Rstoch <- array(0, dim=c(n.int, n.credit, ncol(defaultprob)))
+ for(t in 1:ncol(defaultprob)){
+ for(i in 1:n.credit){
+ Rstoch[,i,t] <- stochasticrecov(recov[i], 0, Z, w, rho, defaultprob[i,t], defaultprobmod[i,t])
+ }
+ }
+ parf <- function(i){
+ pshocked <- apply(defaultprobmod, 2, shockprob, rho=rho, Z=Z[i])
+ S <- 1 - Rstoch[i,,]
+ dist <- lossrecovdist.term(pshocked, 0, issuerweights, S, Ngrid)
+ return( tranche.pvvec(Kmodified, dist$L, dist$R, cs))
+ }
+ clusterExport(cl, list("Rstoch", "p"), envir=environment())
+ result <- parSapply(cl, 1:length(w), parf)
+ return( list(pv=100*(1+result%*%w), pv.w=result))
+}
diff --git a/R/transactions.R b/R/transactions.R new file mode 100644 index 00000000..e36f9c31 --- /dev/null +++ b/R/transactions.R @@ -0,0 +1,26 @@ +transdir <- "U:/cmo_cdu/trans"
+transsave <- "W:/CorpCDOs/Transaction Data/"
+dealnames<-dir()
+stripext <- function(name,extlen=3){
+ substr(name,1,nchar(name)-extlen-1)
+}
+for(i in 1:length(dealnames)){
+ dealnames[i] <- stripext(dealnames[i])
+}
+
+for(i in 1:length(dealnames)){
+ temp <- read.table(paste(dealnames[i],"cdt",sep="."),skip=23,header=T,fill=T,sep="\t",quote="")
+ temp$DATE <- as.Date(strptime(temp$DATE,"%Y%m%d"))
+ assign(dealnames[i],temp)
+}
+
+purchased <- function(dealname){
+ deal <- get(dealname)
+ subset <- which(deal$ACTIVITY=="Purchase")
+ zoo(cbind(deal$AMOUNT[subset],deal$PRICE[subset]),deal$DATE[subset])
+}
+sold <- function(dealname){
+ deal <- get(dealname)
+ subset <- which(deal$ACTIVITY=="Sale")
+ zoo(cbind(deal$AMOUNT[subset],deal$PRICE[subset]),deal$DATE[subset])
+}
diff --git a/R/yieldcurve.R b/R/yieldcurve.R new file mode 100644 index 00000000..a21937b9 --- /dev/null +++ b/R/yieldcurve.R @@ -0,0 +1,63 @@ +library(XML)
+
+getMarkitIRData <- function(date=Sys.Date()) {
+ ## downloads the latest available interest rates data from Markit
+ ## before date and returns the parsed file into a list
+ require(XML)
+ i <- 0
+ temp <- tempfile()
+ while( TRUE ) {
+ lastdate <- format(date-i,"%Y%m%d")
+ filename <- paste("InterestRates_USD_", lastdate, sep="")
+ download.file(paste("http://www.markit.com/news/", filename, ".zip", sep=""), temp, quiet=T)
+ con <- file(temp, "r")
+ firstline <- readLines(con, 1)
+ # Markit returns a plain text file if there is no data.
+ if(firstline == "Interest Rates not available, please check date entered") {
+ i <- i + 1
+ close(con)
+ } else {
+ cat("downloaded data for:", lastdate,"\n")
+ close(con)
+ # we unzip it
+ unzip(temp)
+ unlink(temp)
+ break
+ }
+ }
+ return( xmlToList(paste(filename,".xml", sep="")) )
+}
+
+buildMarkitYC <- function(MarkitData, dt=0.25){
+ settledate <- as.Date(MarkitData$deposits$spotdate)
+
+ params <- list(tradeDate=as.Date(MarkitData$effectiveasof,format="%Y-%m-%d"),
+ settleDate=settledate,
+ dt=dt,
+ interpWhat="discount",
+ interpHow="loglinear")
+
+ tsQuotes <- list(d1m=as.numeric(MarkitData$deposits[5]$curvepoint$parrate),
+ #d2m=as.numeric(MarkitData$deposits[6]$curvepoint$parrate),
+ d3m=as.numeric(MarkitData$deposits[7]$curvepoint$parrate),
+ ## d6m=as.numeric(MarkitData$deposits[8]$curvepoint$parrate),
+ ## d9m=as.numeric(MarkitData$deposits[9]$curvepoint$parrate),
+ ## d1y=as.numeric(MarkitData$deposits[10]$curvepoint$parrate),
+ s2y=as.numeric(MarkitData$swaps[8]$curvepoint$parrate),
+ s3y=as.numeric(MarkitData$swaps[9]$curvepoint$parrate),
+ #s4y=as.numeric(MarkitData$swaps[10]$curvepoint$parrate),
+ s5y=as.numeric(MarkitData$swaps[11]$curvepoint$parrate),
+ #s6y=as.numeric(MarkitData$swaps[12]$curvepoint$parrate),
+ #s7y=as.numeric(MarkitData$swaps[13]$curvepoint$parrate),
+ #s8y=as.numeric(MarkitData$swaps[14]$curvepoint$parrate),
+ #s9y=as.numeric(MarkitData$swaps[15]$curvepoint$parrate),
+ s10y=as.numeric(MarkitData$swaps[16]$curvepoint$parrate),
+ #s12y=as.numeric(MarkitData$swaps[17]$curvepoint$parrate),
+ s15y=as.numeric(MarkitData$swaps[18]$curvepoint$parrate),
+ s20y=as.numeric(MarkitData$swaps[19]$curvepoint$parrate),
+ #s25y=as.numeric(MarkitData$swaps[20]$curvepoint$parrate),
+ s30y=as.numeric(MarkitData$swaps[21]$curvepoint$parrate))
+ YC.USD <- list(params=params, tsQuotes=tsQuotes)
+ return( YC.USD )
+}
+
|