From 8966b9753f876250182ef75a8c4294d523f2126a Mon Sep 17 00:00:00 2001 From: Guillaume Horel Date: Thu, 8 May 2014 13:40:58 -0400 Subject: implement skew mapping by doing tranche loss proportion mapping --- R/tranche_functions.R | 49 ++++++++++++++++++++++++++++++++++++------ src/lossdistrib.c | 59 +++++++++++++++++++++++++++++++++++++++++++++++---- 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); +} -- cgit v1.2.3-70-g09d2