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