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[2:length(args)]
}else{
    data <- read.table(file.path(root.dir, "scripts", "scenarios.txt"))
    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"))

prevBusDay <- function(workdate = Sys.Date()){
    i <- 1
    while(!isBusinessDay(calendar = "UnitedStates/GovernmentBond",
                         workdate - i)){
        i <- i+1
    }
    return( workdate - i )
}

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
j <- 1
## reinvestment parameters
## need to match parameters in build_portfolios.R
reinvspread <- 0.025
reinvfixed <- 0.07
rollingmaturity <- 84 ##months
reinvlag <- 3 ##months
n.scenarios <- 100
recov.adj <- 0.925

for(deal.name in dealnames){
    load(file.path(root.dir, "Scenarios", paste("Portfolios", workdate, sep="_"),
                   paste(deal.name, "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(deal.name)
    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))
    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
    }
    recoverymonthly <- pmin(recoverymonthly,1)
    ## compute reinvestment price
    if(reinvflags[j]){
        DC <- DiscountCurve(L3m$params, L3m$tsQuotes, yearFrac(L3m$params$tradeDate, deal.datesmonthlylagged))
        df <- DC$discounts
        forwards <- DC$forwards
        reinvprices <- compute.reinvprices(df, forwards, cdrmonthly, recoverymonthly, reinvspread,
                                           reinvfixed, rollingmaturity, reinvlag)
    }
    ## 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] + 0.07577 #this should be principal account balance
    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]/(0.85*loanprices[1]+0.15*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]/(0.85*loanprices[t]+0.15*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(deal.name,"-cdr.csv")),
                row.names=F, col.names=F, sep=",")
    write.table(100 * recoverymonthly,
                file=file.path(save.dir, paste0(deal.name,"-recovery.csv")),
                row.names=F, col.names=F, sep=",")
    if(reinvflags[j]){
        write.table(100 * reinvprices$loan,
                    file = file.path(save.dir, paste0(deal.name, "-floatreinvprices.csv")),
                    row.names=F, col.names=F, sep=",")
        write.table(100 * reinvprices$bond,
                    file = file.path(save.dir, paste0(deal.name, "-fixedreinvprices.csv")),
                    row.names=F, col.names=F, sep=",")
    }
    save(scenariosd, scenariosr, dist.joint, file=file.path(save.dir, paste0(deal.name, ".RData")),
         compress="xz")

    cat("generated scenarios for:", deal.name, "\n")
    j <- j+1
}