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 HY17
## 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

hy17portfolio <- c()
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))
## }

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,5:9])*0.01,
                         running=rep(nondefaulted$running[i]*bps, 5))
    SC@curve <- cdshazardrate(quotes, nondefaulted$recovery[i]/100)
    hy17portfolio <- c(hy17portfolio, SC)
}

hy17$indexref <- 1.035
hy17portfolio.tweaked <- tweakcurves(hy17portfolio, hy17)

SurvProb <- SPmatrix(hy17portfolio.tweaked, hy17)

## calibrate the tranches using base correlation
K <- c(0, 0.15, 0.25, 0.35, 1)
Kmodified <- adjust.attachments(K, hy17$loss, hy17$factor)
tranche.upf <- c(49.625, 94.75, 107.125, 114.875)
tranche.running <- c(0.05, 0.05, 0.05, 0.05)


lu <- 0.01
recov <- sapply(hy17portfolio.tweaked, attr, "recovery")

rhovec <- c()
cs <- couponSchedule(nextIMMDate(today()), hy17$maturity,"Q", "FIXED", 0.05, 0)

for(i in 2:length(Kmodified)){
    f <- function(rho, ...){
        temp <- BClossdistC(SurvProb, recov, rho, lu, 100)
        return( abs(tranche.upf[i-1]-1/(Kmodified[i]-Kmodified[i-1])*
                    (tranche.bp(temp$L, temp$R, cs, 0, Kmodified[i])*Kmodified[i]-
                     tranche.bp(oldtemp$L, oldtemp$R, cs, 0, Kmodified[i-1])*Kmodified[i-1])) )
    }
    rho <- optimize(f, SurvProb, recov, lu, tranche.upf, Kmodified, cs, oldtemp)$minimum
    oldtemp <- BClossdistC(SurvProb, rho, recov, lu)
    rhovec <- c(rhovec, rho)
}

#calibrate by modifying the factor distribution
bottomup <- 1:3
topdown <- 2:4
n.int <- 100
n.credit <- 96
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", "rho", "Z", "lossrecovdist.term",
                       "lossrecovdist", "lossdistribC", "lu",
                       "tranche.bpvec", "tranche.bp", "tranche.pl", "tranche.cl",
                       "trancheloss", "trancherecov", "pos", "Kmodified", "cs"))

parf <- function(i){
    pshocked <- apply(p, 2, shockprob, rho=rho, Z=Z[i])
    S <- 1 - Rstoch[i,,]
    dist <- lossrecovdist.term(pshocked, 0, S, lu)
    return( tranche.bpvec(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(result[topdown,], w, tranche.upf[topdown])

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

    ## update the new probabilities
    for(i in 1:n.credit){
        for(j in 1:ncol(p)){
            p[i,j] <- fit.prob(Z, program$weight, rho, defaultprob[i,j])
        }
    }

    errvec <- c(errvec, err)
    w.mod <- program$weight
    cat(err,"\n")
}

## 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(1/lu+1, 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, S, lu)$L
}

lossdist.orig <- BClossdistC(SurvProb, rho, recov, lu)
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])
    }
}
lu <- 0.01
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)], 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]
}

## for(i in seq(0,360, by=10)){
##     persp(distw, theta=i)
##     Sys.sleep(0.5)
## }

persp(distw, theta=-20, phi=13, col="lightgreen",ticktype="detailed", shade=0.3)
x <- seq(0,1,lu)
y <- seq(0,1,lu)

lossdist.bu <- lossdist
recovdist.bu <- recovdist

prepayprob <- matrix(0,100,20)
defaultprob <- matrix(0,100,20)
for( t in 1:20){
  prepayprob[,t] <- h/(h+lambda)*(1-exp(-(lambda+h)*yearfrac[t]))
  defaultprob[,t] <- lambda/(h+lambda)*(1-exp(-(lambda+h)*yearfrac[t]))
}

#LCDX15 calibration
library(statmod)
n.int <- 100
n.credit <- 100
Z <- gauss.quad.prob(n.int, "normal")$nodes
w <- gauss.quad.prob(n.int, "normal")$weights
rho <- 0.35
T <- length(yearfrac)
Recov <- rep(0.35, n.credit)e
dt <- diff(c(0, yearfrac))
K <- c(0, 0.08, 0.15, 0.3, 1)
lu <- 0.01
errvec <- c()
for(l in 1:10){
  result <- c()
  test <- c()
  Rstoch <- array(0,dim=c(T,n.int,n.credit))
  for(t in 1:T){
    for(i in 1:n.credit){
      Rstoch[t,,i] <- stochasticrecov(Recov[i],0,Z,w.bak,rho,defaultprob[i,t], defaultprob[i,t])
    }
  }
  lossdist <- c()
  prepaydist <- c()
  result <- c()
  for(i in 1:n.int){
    dpshocked <- apply(defaultprob,2,shockprob,rho=rho,Z=Z[i])
    ppshocked <- apply(prepayprob,2,shockprob,rho=rho,Z=-Z[i])
    L <- c()
    R <- c()
    for(t in 1:length(yearfrac)){
      dist <- lossrecovdist(dpshocked[,t],ppsocked[,t],Rstoch[t,i,],lu)
      #dist <- lossrecovdist(dpshocked[,t],ppshocked[,t],0.35,lu)
      L <- cbind(L,dist$L)
      R <- cbind(R,dist$R)
    }
    lossdist <- cbind(lossdist,L[,T])
    prepaydist <- cbind(prepaydist,R[,T])
    result <- rbind(result,bpvec(K,L,R,df,dt))
  }

  pomme <- KLfit(t(result), w, tranchemarks)
  #pomme <- lmconst(tranchemarks,t(result),rbind(rep(1,100),Z),c(1,0),T)

  err <-0
  for(i in 1:dim(p.bak)[1]){
    for(j in 1:dim(p.bak)[2]){
      err <- err+abs(crossprod(shockprob(p[i,j],rho,Z),pomme$weight)-p.bak[i,j])
    }
  }

  ptilde <- defaultprob
  for(i in 1:dim(defaultprob)[1]){
    for(j in 1:dim(defaultprob)[2]){
      ptilde[i,j] <- fit.prob(Z,pomme$weight,rho,defaultprob[i,j])
    }
  }
  p <- ptilde
  errvec <- c(errvec,err)
  w.bak <- pomme$weight
  cat(err,"\n")
}

pl<-c()
for(i in 1:100){
  pd<-c(defaultprob[i,1],diff(defaultprob[i,]))
  pl<-c(pl,sum(pd*df*(1-0.35)))
}
cl<-c()
for(i in 1:100){
  cl<-c(cl,sum(df*dt*(1-defaultprob[i,]-prepayprob[i,])))
}

loss <- function(t,P,R,lambda,mu){
  lambda/(lambda+mu)*(1-exp(-(lambda+mu)*t))*(1-R/P)-mu/((lambda+mu)*P)*(1-exp(-(lambda+mu)*t))
}
losstest <- function(t){
  loss(t,P,R,lambda,mu)
}


##tranche pricing
library(statmod)
n.int <- 500
Z <- gauss.quad.prob(n.int,"normal")$nodes
w <- gauss.quad.prob(n.int,"normal")$weights
rho <- 0.35

Lt <- c()
Rt <- c()
for(t in 1:21){
  L <- c()
  R <- c()
  for(i in 1:n.int){
    pZ <- shockprob(1-p[,t],rho,Z[i])
    SZ <- shockseverity(1-p[,t],0.45,1,Z[i])
    temp <- lossrecovdist(pZ,0,1-SZ,0.0045)
    L <- cbind(L,temp$L)
    R <- cbind(R,temp$R)
  }
  Lt <- cbind(Lt, as.matrix(L)%*%w)
  Rt <- cbind(Rt, as.matrix(R)%*%w)
}

hazard.rates <- runif(100)*1000
yearfrac <- seq(0,7,0.25)
bps <- 1e-4
P <- exp(-outer(hazard.rates * bps, yearfrac))
R <- rep(0.6,100)
lu <- 0.01
test <- lossrecovdist(1-P[,29], rep(0,100), R, 0.01)


library(statmod)
n.int <- 500
Z <- gauss.quad.prob(n.int,"normal")$nodes
w <- gauss.quad.prob(n.int,"normal")$weights
rho <- 0.35


cl <- makeCluster(6)
x <- rnorm(1000000)
clusterCall(cl, function(x) for (i in 1:100) sum(x), x)