1 files changed, 56 insertions, 0 deletions
diff --git a/reinvestingdistribution.R b/reinvestingdistribution.R
new file mode 100644
index 00000000..e2c42e3e
--- /dev/null
+++ b/reinvestingdistribution.R
@@ -0,0 +1,56 @@
+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)
+}