implement skew mapping by doing tranche loss proportion mapping

2 files changed, 98 insertions, 10 deletions

diff --git a/R/tranche_functions.R b/R/tranche_functions.R
index b575c33..a5b99d8 100644
--- a/R/tranche_functions.R
+++ b/R/tranche_functions.R
@@ -651,7 +651,7 @@ BClossdistC <- function(defaultprob, issuerweights, recov, rho, Z, w,
     L <- matrix(0, N, dim(defaultprob)[2])

     R <- matrix(0, N, dim(defaultprob)[2])

     rho <- rep(rho, length(issuerweights))

-    r <- .C("BClossdist", defaultprob, dim(defaultprob)[1], dim(defaultprob)[2],

+    r <- .C("BCloss_recov_dist", defaultprob, dim(defaultprob)[1], dim(defaultprob)[2],

             as.double(issuerweights), as.double(recov), as.double(Z), as.double(w),

             as.integer(length(Z)), as.double(rho), as.integer(N), as.logical(defaultflag), L=L, R=R)

     return(list(L=r$L,R=r$R))

@@ -718,17 +718,54 @@ BCtranche.delta <- function(portfolio, index, coupon, K1, K2, rho1, rho2, Z, w,
     return( delta )

 }

 

-BCstrikes <- function(defaultprob, issuerweights, recov, cs, Kmodified, rho, Z, w, N=101) {

+EL <- function(defaultprob, issuerweights, recov, cs){

+    ## computes the expected loss of a portfolio (time discounted)

+    ## given the default curves and recovery

+    ## should be very close to the protection leg of the portfolio of cds

+    ELvec <- crossprod( issuerweights * (1-recov), defaultprob)

+    return( crossprod(cs$df, diff(c(0, as.numeric(ELvec)))) )

+}

+

+ELtrunc <- function(defaultprob, issuerweights, recov, cs, K, rho, Z, w, Ngrid){

+    ## computes the expected loss of a portfolio below strike K

+    ## could be written faster by using a truncated version of lossdistrib

+    rho <- rep(rho, length(issuerweights))

+    L <- matrix(0, Ngrid, dim(defaultprob)[2])

+    L <- .C("BCloss_dist", defaultprob, dim(defaultprob)[1], dim(defaultprob)[2],

+            as.double(issuerweights), as.double(recov), as.double(Z), as.double(w),

+            as.integer(length(Z)), as.double(rho), as.integer(Ngrid), FALSE, L=L)$L

+    support <- seq(0, 1, length=Ngrid)

+    trancheloss <- pmin(support, K)

+    ELvec <- crossprod(trancheloss, L)

+    return( crossprod(cs$df, diff(c(0, as.numeric(ELvec)))) )

+}

+

+BCstrikes <- function(defaultprob, issuerweights, recov, cs, K, rho, Z, w, N=101) {

     ## computes the strikes as a percentage of expected loss

     ## Kmodified is the current attachment points (adjusted for losses)

-    EL <- c()

+    el <- EL(defaultprob, issuerweights, recov, cs)

+    ELvec <- c()

     for(i in 2:length(K)){

-        EL <- c(EL, -BCtranche.pv(defaultprob, issuerweights, recov, cs,

-                                  K[i-1], K[i], rho[i-1], rho[i], Z, w, N)$pl)

+        ELvec <- c(ELvec, ELtrunc(defaultprob, issuerweights, recov, cs, K[i], rho[i], Z, w, N))

     }

-    return(cumsum(EL*diff(Kmodified))/sum(EL*diff(Kmodified)))

+    return( ELvec/el )

 }

 

+skewmapping <- function(defaultprob1, issuerweights1, recov1, cs1, rhofun,

+                        defaultprob2, issuerweights2, recov2, cs2, K2, Z, w, N=101) {

+

+    EL1 <- EL(defaultprob1, issuerweights1, recov1, cs1)

+    EL2 <- EL(defaultprob2, issuerweights2, recov2, cs2)

+    f <- function(x, ...){

+        return(abs(ELtrunc(defaultprob1, issuerweights1, recov1, cs1,

+                           0, x, 0, rhofun(x), Z, w, N)/EL1-

+                   ELtrunc(defaultprob2, issuerweights2, recov2, cs2,

+                           0, K2, 0, rhofun(x), Z, w, N)/EL2))

+    }

+    return( optimize(f, interval=c(0,1)) )

+}

+

+

 delta.factor <- function(K1, K2, index){

     ## compute the factor to convert from delta on current notional to delta on original notional

     ## K1 and K2 original strikes

diff --git a/src/lossdistrib.c b/src/lossdistrib.c
index f12ab2f..e5f2c53 100644
--- a/src/lossdistrib.c
+++ b/src/lossdistrib.c
@@ -728,10 +728,10 @@ void lossdistrib_joint_Z(double *dp, int *ndp, double *w,
   free(qmat);

 }

 

-void BClossdist(double *defaultprob, int *dim1, int *dim2,

-                double *issuerweights, double *recov, double *Z, double *w,

-                int *n, double *rho, int *N, int *defaultflag,

-                double *L, double *R) {

+void BCloss_recov_dist(double *defaultprob, int *dim1, int *dim2,

+		       double *issuerweights, double *recov, double *Z, double *w,

+		       int *n, double *rho, int *N, int *defaultflag,

+		       double *L, double *R) {

   /*

     computes the loss and recovery distribution over time with a flat gaussian

     correlation

@@ -788,3 +788,54 @@ void BClossdist(double *defaultprob, int *dim1, int *dim2,
   free(Lw);

   free(Rw);

 }

+

+

+void BCloss_dist(double *defaultprob, int *dim1, int *dim2,

+      double *issuerweights, double *recov, double *Z, double *w,

+      int *n, double *rho, int *N, int *defaultflag,

+      double *L) {

+  /*

+    computes the loss distribution over time with a flat gaussian

+    correlation

+    inputs:

+      defaultprob: matrix of size dim1 x dim2. dim1 is the number of issuers

+      and dim2 number of time steps

+      issuerweights: vector of issuer weights (length dim1)

+      recov: vector of recoveries (length dim1)

+      Z: vector of factor values (length n)

+      w: vector of factor weights (length n)

+      rho: correlation beta vector (length dim1)

+      N: number of ticks in the grid

+      defaultflag: if true, computes the default distribution

+    outputs:

+      L: matrix of size (N, dim2)

+   */

+  int t, i, j;

+  double g;

+  double *gshocked, *Sshocked, *Lw;

+  int one = 1;

+  double alpha = 1;

+  double beta = 0;

+

+  gshocked = malloc((*dim1) * (*n) * sizeof(double));

+  Sshocked = malloc((*dim1) * (*n) * sizeof(double));

+  Lw = malloc((*N) * (*n) * sizeof(double));

+

+  for(t=0; t < (*dim2); t++) {

+    memset(Lw, 0, (*N) * (*n) * sizeof(double));

+    #pragma omp parallel for private(j, g)

+    for(i=0; i < *n; i++){

+      for(j=0; j < (*dim1); j++){

+        g = defaultprob[j + (*dim1) * t];

+        gshocked[j+(*dim1)*i] = shockprob(g, rho[j], Z[i], 0);

+        Sshocked[j+(*dim1)*i] = shockseverity(1-recov[j], Z[i], rho[j], g);

+      }

+      lossdistrib_blas(gshocked + (*dim1) * i, dim1, issuerweights, Sshocked + (*dim1)*i, N, defaultflag,

+		       Lw + i * (*N));

+    }

+    dgemv_("n", N, n, &alpha, Lw, N, w, &one, &beta, L + t * (*N), &one);

+  }

+  free(gshocked);

+  free(Sshocked);

+  free(Lw);

+}