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