From 48eeb4e9c3a47c8d161e717f7f5d3db997967fac Mon Sep 17 00:00:00 2001
From: Guillaume Horel <guillaume.horel@serenitascapital.com>
Date: Wed, 4 Jun 2014 15:57:06 -0400
Subject: compute gammas as well, hopefully this is right

---
 R/tranche_functions.R | 36 ++++++++++++++++++++++++++++++------
 1 file changed, 30 insertions(+), 6 deletions(-)

(limited to 'R/tranche_functions.R')

diff --git a/R/tranche_functions.R b/R/tranche_functions.R
index 8b8b263..b21094f 100644
--- a/R/tranche_functions.R
+++ b/R/tranche_functions.R
@@ -658,7 +658,7 @@ BClossdistC <- function(defaultprob, issuerweights, recov, rho, Z, w,
 }
 
 BCtranche.pv <- function(index, K1, K2, rho1, rho2,
-                         Z, w, N=length(issuerweights)+1){
+                         Z, w, N=length(index$issuerweights)+1){
     ## computes the protection leg, couponleg, and bond price of a tranche
     ## in the base correlation setting
     if(K1==0){
@@ -679,7 +679,7 @@ BCtranche.pv <- function(index, K1, K2, rho1, rho2,
                 bp=100*(1+(pl2-pl1+cl2-cl1)/dK)))
 }
 
-BCtranche.delta <- function(index, Z, w, N=length(portolio)+1, tradedate = Sys.Date()){
+BCtranche.delta <- function(index, Z, w, N=length(index$portolio)+1, tradedate = Sys.Date()){
     ## computes the tranche delta (on current notional) by doing a proportional
     ## blip of all the curves
     ## if K2==1, then computes the delta using the lower attachment only
@@ -688,28 +688,52 @@ BCtranche.delta <- function(index, Z, w, N=length(portolio)+1, tradedate = Sys.D
 
     Portf.plus <- tweakportfolio(index$portfolio, eps)
     Portf.minus <- tweakportfolio(index$portfolio, -eps)
+    Portf.2plus <- tweakportfolio(index$portfolio, 2*eps)
     indexplus <- index
     indexplus$portfolio <- Portf.plus
     indexplus$defaultprob <- 1 - SPmatrix(Portf.plus, length(index$cs$dates))
+
+    index2plus <- index
+    index2plus$portfolio <- Portf.2plus
+    index2plus$defaultprob <- 1 - SPmatrix(Portf.2plus, length(index$cs$dates))
+
     indexminus <- index
     indexminus$portfolio <- Portf.minus
     indexminus$defaultprob <- 1 - SPmatrix(Portf.minus, length(index$cs$dates))
     dPVindex <- indexpv(indexplus, tradedate = tradedate, clean = FALSE)$bp -
         indexpv(indexminus, tradedate = tradedate, clean = FALSE)$bp
-    deltas <- rep(0, length(index$K)-1)
-    bpplus <- c(0)
-    bpminus <- c(0)
+    dPVindex2 <- indexpv(index2plus, tradedate = tradedate, clean = FALSE)$bp -
+        indexpv(index, tradedate = tradedate, clean = FALSE)$bp
+    dPVindex1 <- indexpv(indexplus, tradedate = tradedate, clean = FALSE)$bp-
+        indexpv(index, tradedate = tradedate, clean = FALSE)$bp
+    bpplus <- 0
+    bpminus <- 0
+    bp <- 0
+    bp2plus <- 0
     for(i in 1:(length(index$rho)-2)){
+        dist <- BClossdistC(index$defaultprob, index$issuerweights, index$recov, index$rho[i+1], Z, w, N)
         distplus <- BClossdistC(indexplus$defaultprob, index$issuerweights, index$recov, index$rho[i+1], Z, w, N)
         distminus <- BClossdistC(indexminus$defaultprob, index$issuerweights, index$recov, index$rho[i+1], Z, w, N)
+        dist2plus <- BClossdistC(index2plus$defaultprob, index$issuerweights, index$recov, index$rho[i+1], Z, w, N)
+        bp <- c(bp, tranche.pv(dist$L, dist$R, index$cs, 0, index$K[i+1]))
         bpplus <- c(bpplus, tranche.pv(distplus$L, distplus$R, index$cs, 0, index$K[i+1]))
         bpminus <- c(bpminus, tranche.pv(distminus$L, distminus$R, index$cs, 0, index$K[i+1]))
+        bp2plus <- c(bp2plus, tranche.pv(dist2plus$L, dist2plus$R, index$cs, 0, index$K[i+1]))
     }
     dPVtranche <- 1/diff(index$K[-length(index$K)])*(diff(bpplus) - diff(bpminus))
+    dPVtranche2 <-  1/diff(index$K[-length(index$K)])*(diff(bp2plus) - diff(bp))
+
     deltas <- dPVtranche/dPVindex
+    ##might need to check that in case of non 1 factor
     deltas <- c(deltas,
                 (1-(bpplus[length(bpplus)]-bpminus[length(bpminus)])/dPVindex)/(1-index$K[length(index$K)-1]))
-    return( deltas )
+    deltasplus <- dPVtranche2/dPVindex2
+    ##might need to check that in case of non 1 factor
+    deltasplus <- c(deltasplus,
+                (1-(bp2plus[length(bpplus)]-bp[length(bpminus)])/dPVindex2)/(1-index$K[length(index$K)-1]))
+    gammas <- (deltasplus-deltas)/dPVindex1/100
+
+    return( list(deltas=deltas, gammas=gammas) )
 }
 
 EL <- function(index){
-- 
cgit v1.2.3-70-g09d2