path: root/R/build_scenarios.R

library(doParallel)
library(yaml)
hostname <- system("hostname", intern=TRUE)
if(hostname=="debian"){
    registerDoParallel(8)
}else{
    registerDoParallel(4)
}

args <- commandArgs(trailingOnly=TRUE)

if(.Platform$OS.type == "unix"){
  root.dir <- "/home/share/CorpCDOs"
}else{
  root.dir <- "//WDSENTINEL/share/CorpCDOs"
}

if(length(args) >= 1){
  workdate <- as.Date(args[1])
}else{
  workdate <- Sys.Date()
}

if(length(args) >=2){
  dealnames <- args[-1]
  reinvflags <- rep_len(TRUE, length(dealnames))
}else{
  data <- read.table(file.path(root.dir, "scripts", "scenarios.txt"),
                     colClasses=c("character", "logical"))
  dealnames <- data$V1
  reinvflags <- data$V2
  unlink(file.path(root.dir, "scripts", "scenarios.txt"))
}

source(file.path(root.dir, "code", "R", "intex_deal_functions.R"))
source(file.path(root.dir, "code", "R", "index_definitions.R"))
source(file.path(root.dir, "code", "R", "etdb.R"))
source(file.path(root.dir, "code", "R", "cds_functions_generic.R"))
source(file.path(root.dir, "code", "R", "cds_utils.R"))

get.reinvassets <- function(dealname, workdate){
  r <- list()
  sqlstr <- sprintf("select * from et_historicaldealinfo('%s', '%s') where ReinvFlag Is true",
                    dealname, workdate)
  data <- dbGetQuery(dbCon, sqlstr)
  if(nrow(data)>0){
      for(i in 1:nrow(data)){
          r[[data$issuername[i]]] <- list(coupontype=data$fixedorfloat[i], liborfloor=data$liborfloor[i])
      }
  }
  return( r )
}

convert.reinvtoperct <- function(d){
  newd <- list(REINV_TBA1=numeric(0), REINV_TBA2=numeric(0))
  if(is.null(d["REINV_TBA2"]) || d["REINV_TBA1"] == d["REINV_TBA2"]){
    if(d["REINV_TBA1"]=="Float"){
        return(c(1,0))
    }else{
        return(c(0,1))
    }
  }
  return(c(0.85, 0.15))
}

compute.reinvprices <- function(df, forwards, cdrmonthly, recoverymonthly,
                                spread, fixedrate, rollingmat, reinvlag){
  ## reinvlag is in months
  ## rollingmat is in months
  floatbp <- c()
  fixedbp <- c()
  for(i in 0:(ncol(cdrmonthly)-1)){
    index <- (i+1):min(i+rollingmat, ncol(cdrmonthly))
    indexlagged <- index + reinvlag
    if(i==0){
      currdf <- df[index]
      laggeddf <- df[indexlagged]
    }else{
      currdf <- df[index]/df[i]
      laggeddf <- df[indexlagged]/df[i]
    }
    floatcoupon <- (forwards[index]+ spread)/12
    fixedcoupon <- fixedrate/12
    currbalance <- 1 - cdrmonthly[,index]/100/12

    if(i < ncol(cdrmonthly)-1){
      currbalance <- t(apply(currbalance, 1, cumprod))
      recov <- -t(apply(cbind(1, currbalance), 1, diff)) *
        recoverymonthly[,index]
      fixedcl <- cbind(1,currbalance[,-length(index)])%*%(fixedcoupon*currdf)
      floatcl <- cbind(1,currbalance[,-length(index)])%*%(floatcoupon*currdf)
      pl <- currbalance[,dim(currbalance)[2]]*currdf[length(currdf)] +
        recov %*% laggeddf
    }else{
      recov <- -t(apply(cbind(1, currbalance), 1, diff)) *
        recoverymonthly[,index]
      fixedcl <- as.numeric(currbalance*fixedcoupon*currdf)
      floatcl <- as.numeric(currbalance*floatcoupon*currdf)
      pl <- as.numeric(currbalance * currdf + recov * laggeddf)
    }
    floatbp <- cbind(floatbp, pl+floatcl)
    fixedbp <- cbind(fixedbp, pl+fixedcl)
  }
  return( list(loan=floatbp, bond=fixedbp) )
}

calibration.date <- prevBusDay(workdate)
calibration <- read.table(file.path(root.dir, "Scenarios", "Calibration",
                                    paste0("calibration-", calibration.date,".csv")), sep=",", header=T)
Z <- calibration$Z
w <- calibration$w

Ngrid <- 201
MarkitData <- getMarkitIRData(calibration.date)
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))

support <- seq(0, 1, length = Ngrid)
recoverylag <- 90 ##days
## reinvestment parameters
## need to match parameters in build_portfolios.R
reinvspread <- 0.04
reinvfixed <- 0.07
basecase.rollingmaturity <- 84 ##months
reinvlag <- 3 ##months
n.scenarios <- 100
recov.adj <- 1


for(j in seq_along(dealnames)){
  load(file.path(root.dir, "Scenarios", paste("Portfolios", workdate, sep="_"),
                 paste(dealnames[j], "RData", sep=".")))
  if(is.na(deal.data$"Reinv End Date")){
    reinvflags[j] <- FALSE
  }
  dp <- A$DP
  pp <- A$PP
  dpmod <- MFupdate.probC(Z, w, deal.portfolio$beta, dp)
  ppmod <- MFupdate.probC(-Z, w, deal.portfolio$beta, pp)
  dist.joint <- MFlossdist.prepay.joint(w, Z, deal.portfolio$beta, dp, dpmod, pp, ppmod,
                                        deal.weights, 1-S, Ngrid)

  distDR <- dist.transform(dist.joint)
  ## compute E(R|D)
  R <- matrix(0, Ngrid, dim(distDR)[1])
  for(t in 1:dim(distDR)[1]){
    R[,t] <- sweep(distDR[t,,], 1, rowSums(distDR[t,,]), "/") %*% support
    R[1,t] <- 0
    if(t >= 2){
      R[,t] <- pmax(R[,t], R[,t-1])
    }
  }

  ## compute scenariosd
  scenariosd <- matrix(0, n.scenarios, dim(distDR)[1])
  scenariosr <- matrix(0, n.scenarios, dim(distDR)[1])
  percentiles <- seq(0, 1, 1/n.scenarios)
  for(t in 1:dim(distDR)[1]){
    D <- rowSums(distDR[t,,])
    Dfun <- splinefun(c(0, cumsum(D)), c(0, support), "monoH.FC")
    Rfun <- approxfun(support, R[,t], rule=2)
    for(i in 1:n.scenarios){
      ## this is roughtly E(D|D is in ith percentile)
      ## using trapezoidal approximation
      scenariosd[i,t] <- 0.5 * (Dfun((i-1)*0.01)+Dfun(i*0.01))
      if(t>=2 && scenariosd[i,t] < scenariosd[i,t-1]){
        scenariosd[i,t] <- scenariosd[i,t-1]
      }
      scenariosr[i,t] <- Rfun(scenariosd[i,t])
      if(t>=2 && scenariosr[i,t] < scenariosr[i,t-1]){
        scenariosr[i,t] <- scenariosr[i,t-1]
      }
    }
  }

  ## we need to adjust the recovery because intex has some embedded amortization assumptions
  ## that we can't turn off (multiply by recov.adj)
  intexrecov <- matrix(0, n.scenarios, dim(distDR)[1])
  for(i in 1:dim(distDR)[1]){
    if(i==1){
      intexrecov[,i] <- recov.adj * (scenariosr[,i]/scenariosd[,1])
    }else{
      intexrecov[,i] <- recov.adj * (scenariosr[,i]-scenariosr[,i-1])/(scenariosd[,i]-scenariosd[,i-1])
    }
  }

  cdr <- cdrfromscenarios(scenariosd, deal.dates)
  ## linear approximation for monthly scenarios
  deal.data <- getdealdata(dealnames[j])
  deal.datesmonthly <- getdealschedule(deal.data, "1 month")
  deal.datesmonthlylagged <- getdealschedule(deal.data, "1 month", workdate, 92)
  cdrmonthly <- matrix(0, n.scenarios, length(deal.datesmonthly))
  recoverymonthly <- matrix(0, n.scenarios, length(deal.datesmonthly))
  scenariosrmonthly <- matrix(0, n.scenarios, length(deal.datesmonthly))
  scenariosdmonthly <- 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
    scenariosrmonthly[i,] <- approx(deal.dates, scenariosr[i,], deal.datesmonthly, rule=2)$y
    scenariosdmonthly[i,] <- approx(deal.dates, scenariosd[i,], deal.datesmonthly, rule=2)$y
  }
  recoverymonthly <- pmin(recoverymonthly,1)
  recoverymonthly[!is.finite(recoverymonthly)] <- 100
  ## compute reinvestment price
  if(!is.na(deal.data$"Reinv End Date") && deal.data$"Reinv End Date" <= workdate){
      ## we cap rolling maturity at the current weighted average maturity of the portfolio
      rollingmaturity <- (crossprod(sapply(deal.portfolio$SC, creditcurve.maturity),
                                    deal.portfolio$notional)/sum(deal.portfolio$notional)
                          - as.numeric(workdate))/365*12
  }else{
      rollingmaturity <- basecase.rollingmaturity
  }
  DC <- DiscountCurve(L3m$params, L3m$tsQuotes, yearFrac(L3m$params$tradeDate, deal.datesmonthlylagged))
  df <- DC$discounts
  forwards <- DC$forwards

  reinvassets <- get.reinvassets(dealnames[j], workdate)
  reinvprices <- list()
  if(reinvflags[j] && length(reinvassets)>0){
      for(assetname in names(reinvassets)){
          asset <- reinvassets[[assetname]]
          if(asset$coupontype=="FLOAT") {
              coupon <- reinvspread
          }else{
              coupon <- reinvfixed
          }
          ##reinvest tweak
          coupon <- coupon-0.0075

          reinvprices[[assetname]] <- foreach(date = iter(deal.datesmonthly), .combine=c) %dopar% {
              100 * forwardportfolioprice(deal.portfolio, date,
                                          min(date+rollingmaturity*30, deal.data$maturity),
                                          asset$coupontype, coupon, asset$liborfloor/100)
          }
      }
  }
    ## reinvprices <- compute.reinvprices(df, forwards, cdrmonthly, recoverymonthly, reinvspread,
    ##                                    reinvfixed, rollingmaturity, reinvlag)

    ## loanprices <- apply(reinvprices.tweak$loan, 2, mean)
    ## bondprices <- apply(reinvprices.tweak$bond, 2, mean)

    ## reinvassets <- convert.reinvtoperct(get.reinvassets(dealnames[j]))
    ## ## all amounts are in units of current collateral balance
    ## reinvdollar <- matrix(0, n.scenarios, length(deal.datesmonthly))
    ## maturingbalance <- matrix(0, n.scenarios, length(deal.datesmonthly))
    ## reinvdollar[,1] <- scenariosrmonthly[,1] + deal.data$"Principal Bal"/deal.data$"Curr Collat Bal"
    ## loanprices <- crossprod(reinvdollar[,1], reinvprices$loan[,1])/sum(reinvdollar[,1])
    ## bondprices <- c(crossprod(reinvdollar[,1], reinvprices$bond[,1])/sum(reinvdollar[,1]))
    ## maturingbalance[,min(1+rollingmaturity, length(deal.datesmonthly))] <-
    ##   reinvdollar[,1]/crossprod(reinvassets, c(loanprices[1], bondprices[1]))
    ## reinvbalance <- matrix(0, n.scenarios, length(deal.datesmonthly))
    ## for(t in 2:dim(cdrmonthly)[2]){
    ##   reinvdollar <- scenariosrmonthly[,t]-scenariosrmonthly[,t-1] +
    ##     reinvbalance * cdrmonthly[,t]/100 * recoverymonthly[,t] *
    ##       yearFrac(deal.datesmonthly[t-1], deal.datesmonthly[t])
    ##   if(t==(dim(cdrmonthly)[2]-1)){
    ##     reinvbalance[,t] <- maturingbalance[,length(deal.datesmonthly)]
    ##   }else if(t==dim(cdrmonthly)[2]){
    ##     reinvbalance[,t] <- 0
    ##   }else{
    ##     reinvbalance[,t] <- rowSums(maturingbalance[,(t+1):length(deal.datesmonthly)])
    ##   }
    ##   reinvdollar[,t] <- scenariosrmonthly[,t] - scenariosrmonthly[,t-1] +
    ##     reinvbalance[,t] * cdrmonthly[,t]/100 * recoverymonthly[,t] *
    ##       yearFrac(deal.datesmonthly[t-1], deal.datesmonthly[t])
    ##   loanprices <- c(loanprices, crossprod(reinvdollar[,t], reinvprices$loan[,t])/sum(reinvdollar[,t]))
    ##   bondprices <- c(bondprices, crossprod(reinvdollar[,t], reinvprices$bond[,t])/sum(reinvdollar[,t]))
    ##   maturingbalance[,min(t+rollingmaturity, length(deal.datesmonthly))] <-
    ##     reinvdollar[,t]/crossprod(reinvassets, c(loanprices[t], bondprices[t]))
    ##   maturingbalance[,min(t+rollingmaturity, length(deal.datesmonthly))] <-
    ##     maturingbalance[,min(t+rollingmaturity, length(deal.datesmonthly))] +
    ##       reinvdollar[,t]/loanprices[t]
    ## }

  save.dir <- file.path(root.dir, "Scenarios", paste("Intex curves", workdate, sep="_"), "csv")
  if(!file.exists(save.dir)){
    dir.create(save.dir, recursive = T)
  }

  write.table(cdrmonthly,
              file= file.path(save.dir, paste0(dealnames[j],"-cdr.csv")),
              row.names=F, col.names=F, sep=",")
  write.table(100 * recoverymonthly,
              file=file.path(save.dir, paste0(dealnames[j],"-recovery.csv")),
              row.names=F, col.names=F, sep=",", na="NaN")
  write.table(reinvprices, file = file.path(save.dir, paste0(dealnames[j], "-reinvprices.csv")),
              row.names=F, col.names=T, sep=",")
  
  configfile <- file.path(save.dir, paste0(dealnames[j], ".config"))
  config <- list(rollingmat = as.integer(rollingmaturity),
                 reinvflag = reinvflags[j])
  cat(as.yaml(config), file=configfile)
  save(scenariosd, scenariosr, dist.joint, file=file.path(save.dir, paste0(dealnames[j], ".RData")),
       compress="xz")

  cat("generated scenarios for:", dealnames[j], "\n")
}