diff options
| -rw-r--r-- | lossdistrib.c | 24 | ||||
| -rw-r--r-- | tranche_functions.R | 98 |
2 files changed, 43 insertions, 79 deletions
diff --git a/lossdistrib.c b/lossdistrib.c index e6eef029..d9af0f89 100644 --- a/lossdistrib.c +++ b/lossdistrib.c @@ -53,7 +53,7 @@ void lossdistrib(double *p, int *np, double *w, double *S, int *N, double *q) { free(q2);
}
-void lossdistrib_fast(double *p, int *np, double *w, double *S, int *N, int *T, double *q) {
+void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N, int *T, double *q) {
/* recursive algorithm with first order correction for computing
the loss distribution.
p vector of default probabilities
@@ -278,7 +278,8 @@ void addandmultiply(double *X, double alpha, double *Y, int n) { }
}
-void BClossdist(double *SurvProb, int *dim1, int *dim2, double *issuerweights, double *recov, double *Z, double *w, int *n, double *rho, int *N, double *L, double *R) {
+void BClossdist(double *SurvProb, int *dim1, int *dim2, double *issuerweights,
+ double *recov, double *Z, double *w, int *n, double *rho, int *N, int *T, double *L, double *R) {
/*
computes the loss and recovery distribution over time with a flat gaussiancorrelation
inputs:
@@ -289,9 +290,10 @@ void BClossdist(double *SurvProb, int *dim1, int *dim2, double *issuerweights, d w: vector of factor weights (length n)
rho: correlation beta
N: number of ticks in the grid
+ T: number of probabilities to compute
outputs:
- L: matrix of size (1/lu+1, dim2)
- R: matrix of size (1/lu+1, dim2)
+ L: matrix of size (T, dim2)
+ R: matrix of size (T, dim2)
*/
int t, i, j;
double g;
@@ -300,8 +302,8 @@ void BClossdist(double *SurvProb, int *dim1, int *dim2, double *issuerweights, d gshocked = calloc((*dim1), sizeof(double));
Rshocked = calloc((*dim1), sizeof(double));
Sshocked = calloc((*dim1), sizeof(double));
- Lw = malloc((*N) * sizeof(double));
- Rw = malloc((*N) * sizeof(double));
+ Lw = malloc((*T) * sizeof(double));
+ Rw = malloc((*T) * sizeof(double));
for(t=0; t < (*dim2); t++) {
for(i=0; i < *n; i++){
@@ -312,14 +314,14 @@ void BClossdist(double *SurvProb, int *dim1, int *dim2, double *issuerweights, d Rshocked[j] = 1 - Sshocked[j];
}
/* reset Lw and Rw to 0 */
- for(j=0; j<*N; j++){
+ for(j=0; j<*T; j++){
Lw[j]=0;
Rw[j]=0;
}
- lossdistrib(gshocked, dim1, issuerweights, Sshocked, N, Lw);
- lossdistrib(gshocked, dim1, issuerweights, Rshocked, N, Rw);
- addandmultiply(Lw, w[i], L + t * (*N), *N);
- addandmultiply(Rw, w[i], R + t * (*N), *N);
+ lossdistrib_truncated(gshocked, dim1, issuerweights, Sshocked, N, T, Lw);
+ lossdistrib_truncated(gshocked, dim1, issuerweights, Rshocked, N, T, Rw);
+ addandmultiply(Lw, w[i], L + t * (*T), *T);
+ addandmultiply(Rw, w[i], R + t * (*T), *T);
}
}
free(gshocked);
diff --git a/tranche_functions.R b/tranche_functions.R index 69b14511..3c41d45f 100644 --- a/tranche_functions.R +++ b/tranche_functions.R @@ -43,51 +43,24 @@ lossdistrib2 <- function(p, w, S, N){ q2 <- c(rep(0,d2), p2*q[1:(N-d2)])
q <- q1 + q2 + (1-p[i])*q
}
- ## q[length(q)] <- q[length(q)]+1-sum(q)
- return(q)
-}
-
-lossdistrib2.fast <- function(p, w, S, N){
- #recursive algorithm with first order correction
- #p vector of default probabilities
- #w vector of weigths
- #S vector of severities
- #N number of ticks in the grid
- ## this is actually slower than lossdistrib2. But in C this is
- ## twice as fast.
- n <- length(p)
- lu <- 1/(N-1)
- q <- rep(0, N)
- q[1] <- 1
- M <- 1
- for(i in 1:n){
- d <- S[i] * w[i] / lu
- d1 <- floor(d)
- d2 <- ceiling(d)
- p1 <- p[i]*(d2-d)
- p2 <- p[i] - p1
- q1 <- p1*q[1:M]
- q2 <- p2*q[1:M]
- q[1:M] <- (1-p[i])*q[1:M]
- q[(d1+1):(M+d1)] <- q[(d1+1):(M+d1)]+q1
- q[(d2+1):(M+d2)] <- q[(d2+1):(M+d2)]+q2
- M <- M+d2
- }
- ## q[length(q)] <- q[length(q)]+1-sum(q)
+ q[length(q)] <- q[length(q)]+1-sum(q)
return(q)
}
-lossdistrib2.fast.trunc <- function(p, w, S, N, truncated=1){
+lossdistrib2.truncated <- function(p, w, S, N, cutoff=N){
## recursive algorithm with first order correction
## p vector of default probabilities
## w vector of weigths
## S vector of severities
- ## N number of ticks in the grid
- ## truncated: where to stop computing the exact probabilities
+ ## N number of ticks in the grid (for best accuracy should
+ ## be a multiple of the number of issuers)
+ ## cutoff where to stop computing the exact probabilities
## (useful for tranche computations)
## this is actually slower than lossdistrib2. But in C this is
## twice as fast.
+ ## for high severities, M can become bigger than N, and there is
+ ## some probability mass escaping.
n <- length(p)
lu <- 1/(N-1)
q <- rep(0, truncated)
@@ -99,18 +72,16 @@ lossdistrib2.fast.trunc <- function(p, w, S, N, truncated=1){ d2 <- ceiling(d)
p1 <- p[i]*(d2-d)
p2 <- p[i] - p1
- q1 <- p1*q[1:min(M, truncated-d1)]
- q2 <- p2*q[1:min(M, truncated-d2)]
- q[1:min(M, truncated)] <- (1-p[i])*q[1:min(M, truncated)]
- q[(d1+1):min(M+d1, truncated)] <- q[(d1+1):min(M+d1, truncated)]+q1
- q[(d2+1):min(M+d2, truncated)] <- q[(d2+1):min(M+d2, truncated)]+q2
+ q1 <- p1*q[1:min(M, cutoff-d1)]
+ q2 <- p2*q[1:min(M, cutoff-d2)]
+ q[1:min(M, cutoff)] <- (1-p[i])*q[1:min(M, cutoff)]
+ q[(d1+1):min(M+d1, cutoff)] <- q[(d1+1):min(M+d1, cutoff)]+q1
+ q[(d2+1):min(M+d2, cutoff)] <- q[(d2+1):min(M+d2, cutoff)]+q2
M <- M+d2
}
- ## q[length(q)] <- q[length(q)]+1-sum(q)
return(q)
}
-
recovdist <- function(dp, pp, w, S, N){
## computes the recovery distribution for a sum of independent variables
## R=\sum_{i=1}^n X_i
@@ -226,7 +197,7 @@ lossdistribC <- function(p, w, S, N){ as.double(w), as.double(S), as.integer(N), q = double(N))$q
}
-lossdistrib.fastC <- function(p, w, S, N, T=N){
+lossdistribC.truncated <- function(p, w, S, N, T=N){
## C version of lossdistrib2, roughly 50 times faster
dyn.load("lossdistrib.dll")
.C("lossdistrib_fast", as.double(p), as.integer(length(p)),
@@ -336,14 +307,13 @@ trancherecov <- function(R, K1, K2){ pos(R - 1 + K2) - pos(R - 1 +K1)
}
-tranche.cl <- function(L, R, cs, K1, K2, scaled=FALSE){
+tranche.cl <- function(L, R, cs, K1, K2, Ngrid=nrow(L), scaled=FALSE){
## computes the couponleg of a tranche
## if scaled is TRUE, scale it by the size of the tranche (K2-K1)
if(K1==K2){
return( 0 )
}else{
- n <- nrow(L)
- support <- seq(0, 1, length=n)
+ support <- seq(0, 1, length=Ngrid)[1:nrow(L)]
size <- K2 - K1 - crossprod(trancheloss(support, K1, K2), L) -
crossprod(trancherecov(support, K1, K2), R)
sizeadj <- as.numeric(0.5 * (size + c(K2-K1, size[-length(size)])))
@@ -355,14 +325,13 @@ tranche.cl <- function(L, R, cs, K1, K2, scaled=FALSE){ }
}
-tranche.pl <- function(L, R, cs, K1, K2, scaled=FALSE){
+tranche.pl <- function(L, R, cs, K1, K2, Ngrid=nrow(L), scaled=FALSE){
## computes the protection leg of a tranche
## if scaled
if(K1==K2){
return(0)
}else{
- n <- nrow(L)
- support <- seq(0, 1, length=n)
+ support <- seq(0, 1, length=Ngrid)[1:nrow(L)]
cf <- K2 - K1 - crossprod(trancheloss(support, K1, K2), L)
cf <- c(K2 - K1, cf)
if(scaled){
@@ -373,8 +342,8 @@ tranche.pl <- function(L, R, cs, K1, K2, scaled=FALSE){ }
}
-tranche.pv <- function(L, R, cs, K1, K2){
- return( tranche.pl(L, R, cs, K1, K2) + tranche.cl(L, R, cs, K1, K2))
+tranche.pv <- function(L, R, cs, K1, K2, Ngrid=nrow(L)){
+ return( tranche.pl(L, R, cs, K1, K2, Ngrid) + tranche.cl(L, R, cs, K1, K2, Ngrid))
}
adjust.attachments <- function(K, losstodate, factor){
@@ -383,15 +352,16 @@ adjust.attachments <- function(K, losstodate, factor){ return( pmin(pmax((K-losstodate)/factor, 0),1) )
}
-tranche.bpvec <- function(K, L, R, cs){
+tranche.pvvec <- function(K, L, R, cs){
r <- rep(0, length(K)-1)
for(i in 1:(length(K)-1)){
- r[i] <- tranche.bp(L, R, cs, K[i], K[i+1])
+ r[i] <- 1/(K[i+1]-K[i]) * tranche.pv(L, R, cs, K[i], K[i+1])
}
return( r )
}
-BClossdist <- function(SurvProb, issuerweights, recov, rho, N, n.int=100){
+BClossdist <- function(SurvProb, issuerweights, recov, rho, N=length(recov)+1,
+ n.int=100){
quadrature <- gauss.quad.prob(n.int, "normal")
Z <- quadrature$nodes
w <- quadrature$weights
@@ -414,14 +384,17 @@ BClossdist <- function(SurvProb, issuerweights, recov, rho, N, n.int=100){ list(L=L, R=R)
}
-BClossdistC <- function(SurvProb, issuerweights, recov, rho, N, n.int=100){
+BClossdistC <- function(SurvProb, issuerweights, recov, rho,
+ N=length(issuerweights)+1, T=N, n.int=100){
dyn.load("lossdistrib.dll")
quadrature <- gauss.quad.prob(n.int, "normal")
Z <- quadrature$nodes
w <- quadrature$weights
- L <- matrix(0, N, dim(SurvProb)[2])
- R <- matrix(0, N, dim(SurvProb)[2])
- r <- .C("BClossdist", SurvProb, as.integer(dim(SurvProb)[1]), as.integer(dim(SurvProb)[2]), as.double(issuerweights), as.double(recov), as.double(Z), as.double(w), as.integer(n.int), as.double(rho), as.integer(N), L=L, R=R)
+ L <- matrix(0, T, dim(SurvProb)[2])
+ R <- matrix(0, T, dim(SurvProb)[2])
+ r <- .C("BClossdist", SurvProb, as.integer(dim(SurvProb)[1]), as.integer(dim(SurvProb)[2]),
+ as.double(issuerweights), as.double(recov), as.double(Z), as.double(w),
+ as.integer(n.int), as.double(rho), as.integer(N), as.integer(T), L=L, R=R)
return(list(L=r$L,R=r$R))
}
@@ -496,14 +469,3 @@ delta.factor <- function(K1, K2, index){ -adjust.attachments(K1, index$loss, index$factor))/(K2-K1)
return( factor )
}
-
-rho <- 0.999
-N <- 101
-L <- matrix(0, N, 100)
-for(i in 1:100){
- ptilde <- shockprob(p, rho, Z[i])
- Sshocked <- shockseverity(S, 1, Z[i], rho, p)
- L[,i] <- lossdistrib2(ptilde, w, Sshocked, N)
-}
-Lw <- L%*%quadrature$weights
-crossprod(Lw, seq(0, 1, length=N))
|