diff options
| -rw-r--r-- | lossdistrib.c | 58 | ||||
| -rw-r--r-- | tranche_functions.R | 127 |
2 files changed, 148 insertions, 37 deletions
diff --git a/lossdistrib.c b/lossdistrib.c index 62016e22..511d2169 100644 --- a/lossdistrib.c +++ b/lossdistrib.c @@ -121,6 +121,64 @@ void lossdistrib2( double *p, int *np, double *S, double *lu, double *q) { free(qtemp);
}
+void aux(double *q, double *qtemp, double p, int N, int d) {
+ /* helper function for the lossdistrib functions */
+ int j;
+ for(j=d; j<N; j++){
+ qtemp[j] += p * q[j];
+ };
+}
+
+void lossdistrib3(double *dp, double *pp, int *n, double *S, double *w, double *lu, double *q) {
+ /* recursive algorithm with first order correction for computing
+ the loss distribution.
+ dp vector of default probabilities
+ pp vector of prepay probabilities
+ n length of p
+ S vector of severities (should be same length as p)
+ w vector of weights
+ lu loss unit
+ q the loss distribution */
+
+ int i, j, N, d1l, d1u, d2l, d2u;
+ double d1, d2, dp1, dp2, pp1, pp2, sum;
+ double *qtemp;
+
+ N = ceil(1./(*lu) + 1);
+ qtemp = malloc( N * sizeof(double));
+ q[0] = 1;
+ for(i=0; i<(*n); i++){
+ d1 = w[i] * (1-S[i]) / *lu;
+ d2 = w[i]/ *lu;
+ d1l = floor(d1);
+ d1u = d1l + 1;
+ d2l = floor(d2);
+ d2u = d1u + 1;
+ dp1 = dp[i] * (d1u-d1);
+ dp2 = dp[i] - dp1;
+ pp1 = pp[i] * (d2u -d2);
+ pp2 = pp[i] - pp1;
+ for(j=0; j<N; j++){
+ qtemp[j]=0;
+ }
+ aux(q, qtemp, dp1, N, d1l);
+ aux(q, qtemp, dp2, N, d1u);
+ aux(q, qtemp, pp1, N, d2l);
+ aux(q, qtemp, pp2, N, d2u);
+ for(j=0; j<N; j++){
+ q[j] = qtemp[j] + (1-dp[i]-pp[i]) * q[j];
+ };
+ }
+ /* correction for weight loss */
+ sum = 0;
+ for(j=0; j<N; j++){
+ sum += q[j];
+ }
+ q[N-1] += 1-sum;
+ free(q);
+ free(qtemp);
+}
+
double shockprob(double p, double rho, double Z, int give_log){
return( pnorm( (qnorm(p, 0, 1, 1, 0) -sqrt(rho) * Z)/sqrt(1 - rho), 0, 1, 1, give_log));
}
diff --git a/tranche_functions.R b/tranche_functions.R index 9fdd3dcb..1b20fdc7 100644 --- a/tranche_functions.R +++ b/tranche_functions.R @@ -1,7 +1,7 @@ library(statmod)
lossdistrib <- function(p){
- #basic recursive algorithm of Andersen, Sidenius and Basu
+ ## basic recursive algorithm of Andersen, Sidenius and Basu
n <- length(p)
q <- rep(0, (n+1))
q[1] <- 1
@@ -13,19 +13,20 @@ lossdistrib <- function(p){ }
lossdistrib.fft <- function(p){
- #computes loss distribution using the fft
- #complexity is of order O(n*m)+O(m*log(m))
- #where m is the size of the grid and n the number of probabilities.
- #this is slower than the recursive algorithm
- theta <- 2*pi*1i*(0:n)/(n+1)
- Phi <- 1 - p + p%o%exp(theta)
- v <- apply(Phi, 2, prod)
- return(1/(n+1)*Re(fft(v)))
+ ## computes loss distribution using the fft
+ ## complexity is of order O(n*m)+O(m*log(m))
+ ## where m is the size of the grid and n the number of probabilities.
+ ## this is slower than the recursive algorithm
+ theta <- 2*pi*1i*(0:n)/(n+1)
+ Phi <- 1 - p + p%o%exp(theta)
+ v <- apply(Phi, 2, prod)
+ return(1/(n+1)*Re(fft(v)))
}
-lossdistrib2 <- function(p, S, lu){
+lossdistrib2 <- function(p, w, S, lu){
#recursive algorithm with first order correction
#p vector of default probabilities
+ #w vector of weigths
#S vector of severities
#lu loss unit
n <- length(p)
@@ -34,7 +35,7 @@ lossdistrib2 <- function(p, S, lu){ q <- rep(0, N+eps)
q[1] <- 1
for(i in 1:n){
- d <- S[i]/(n*lu)
+ d <- S[i] * w[i] / lu
d1 <- floor(d)
d2 <- ceiling(d)
p1 <- p[i]*(d2-d)
@@ -47,7 +48,7 @@ lossdistrib2 <- function(p, S, lu){ return(q)
}
-lossdistrib3 <- function(dp, pp, w, S, lu){
+recovdist <- function(dp, pp, w, S, lu){
## computes the recovery distribution for a sum of independent variables
## R=\sum_{i=1}^n X_i
## where X_i = 0 w.p 1-dp[i]-pp[i]
@@ -64,13 +65,14 @@ lossdistrib3 <- function(dp, pp, w, S, lu){ d1 <- w[i]*(1-S[i])/lu
d1l <- floor(d1)
d1u <- ceiling(d1)
- d2 <- w[i]/lu
+ d2 <- w[i] / lu
d2l <- floor(d2)
d2u <- ceiling(d2)
dp1 <- dp[i] * (d1u-d1)
dp2 <- dp[i] - dp1
pp1 <- pp[i] * (d2u - d2)
pp2 <- pp[i] - pp1
+ cat(dp1, dp2, pp1, pp2, "\n")
q1 <- c(rep(0, d1l), dp1 * q[1:(N-d1l)])
q2 <- c(rep(0, d1u), dp2 * q[1:(N-d1u)])
q3 <- c(rep(0, d2l), pp1 * q[1:(N-d2l)])
@@ -80,19 +82,20 @@ lossdistrib3 <- function(dp, pp, w, S, lu){ return(q)
}
-lossdistrib.joint <- function(p, S, lu){
- #recursive algorithm with first order correction
- #to compute the joint probability distribition of the loss and recovery
- #p vector of default probabilities
- #S vector of severities
- #lu loss unit
+lossdistrib.joint <- function(p, w, S, lu){
+ ## recursive algorithm with first order correction
+ ## to compute the joint probability distribution of the loss and recovery
+ ## p vector of default probabilities
+ ## w vector of issuer weights
+ ## S vector of severities
+ ## lu loss unit
n <- length(p)
N <- ceiling(1/lu + 1)
q <- matrix(0, N, N)
q[1,1] <- 1
for(k in 1:n){
- x <- S[k]/(n*lu)
- y <- (1-S[k])/(n*lu)
+ x <- S[k] * w[k]/lu
+ y <- (1-S[k]) * w[k]/lu
i <- floor(x)
j <- floor(y)
weights <- c((i+1-x)*(j+1-y), (i+1-x)*(y-j), (x-i)*(y-j), (j+1-y)*(x-i))
@@ -108,34 +111,84 @@ lossdistrib.joint <- function(p, S, lu){ return(q)
}
-lossdistribC <- function(p, S, lu){
- #C version of lossdistrib2, roughly 50 times faster
+lossdistribprepay.joint <- function(dp, pp, w, S, lu){
+ ## recursive algorithm with first order correction
+ ## to compute the joint probability distribition of the loss and recovery
+ ## dp vector of default probabilities
+ ## pp vector of prepay probabilities
+ ## w vector of issuer weights
+ ## S vector of severities
+ ## lu loss unit
+ n <- length(dp)
+ N <- ceiling(1/lu + 1)
+ q <- matrix(0, N, N)
+ q[1,1] <- 1
+ for(k in 1:n){
+ x <- S[k] * w[k]/lu
+ y1 <- (1-S[k]) * w[k]/lu
+ y2 <- w[k]/lu
+ i <- floor(x)
+ j1 <- floor(y1)
+ j2 <- floor(y2)
+ weights <- c((i+1-x)*(j1+1-y1), (i+1-x)*(y1-j1), (x-i)*(y1-j1), (j1+1-y1)*(x-i))
+ dpsplit <- dp[k] * weights
+ ppsplit <- pp[k] * c(j2+1-y2, y2-j2)
+ qtemp <- matrix(0, N, N)
+ qtemp[(i+1):N,(j1+1):N] <- qtemp[(i+1):N,(j1+1):N] + dpsplit[1] * q[1:(N-i),1:(N-j1)]
+ qtemp[(i+1):N,(j1+2):N] <- qtemp[(i+1):N,(j1+2):N] + dpsplit[2] * q[1:(N-i), 1:(N-j1-1)]
+ qtemp[(i+2):N,(j1+2):N] <- qtemp[(i+2):N,(j1+2):N] + dpsplit[3] * q[1:(N-i-1), 1:(N-j1-1)]
+ qtemp[(i+2):N,(j1+1):N] <- qtemp[(i+2):N, (j1+1):N] + dpsplit[4] * q[1:(N-i-1), 1:(N-j1)]
+ qtemp[, (j2+1):N] <- qtemp[,(j2+1):N]+ppsplit[1]*q[,1:(N-j2)]
+ qtemp[, (j2+2):N] <- qtemp[,(j2+2):N]+ppsplit[2]*q[,1:(N-j2-1)]
+ q <- qtemp + (1-pp[k]-dp[k]) * q
+ }
+ q[length(q)] <- q[length(q)] + 1 - sum(q)
+ return(q)
+}
+
+
+lossdistribC <- function(p, w, S, lu){
+ ## C version of lossdistrib2, roughly 50 times faster
dyn.load("lossdistrib.dll")
.C("lossdistrib", as.double(p), as.integer(length(p)),
- as.double(S), as.double(lu), q = double(1/lu+1))$q
+ as.double(w), as.double(S), as.double(lu), q = double(1/lu+1))$q
}
-lossdistribC.joint <- function(p, S, lu){
- #C version of lossdistrib.joint, roughly 20 times faster
+recovdistC <- function(dp, pp, w, S, lu){
+ ## C version of recovdist
+ dyn.load("lossdistrib2.dll")
+ .C("recovdist", as.double(dp), as.double(pp), as.integer(length(dp)),
+ as.double(w), as.double(S), as.double(lu), q = double(ceiling(1/lu+1)))$q
+}
+
+lossdistribC.joint <- function(p, w, S, lu){
+ ## C version of lossdistrib.joint, roughly 20 times faster
N <- ceiling(1/lu+1)
- dyn.load("lossdistrib.dll")
- .C("lossdistrib2", as.double(p), as.integer(length(p)), as.double(S),
- as.double(lu), q = matrix(0, N, N))$q
+ dyn.load("lossdistrib2.dll")
+ .C("lossdistrib_joint", as.double(p), as.integer(length(p)), as.double(w),
+ as.double(S), as.double(lu), q = matrix(0, N, N))$q
+}
+
+lossdistribprepayC.joint <- function(dp, pp, w, S, lu){
+ ## C version of lossdistribprepay.joint
+ dyn.load("lossdistrib2.dll")
+ N <- ceiling(1/lu+1)
+ .C("lossdistrib_prepay_joint", as.double(dp), as.double(pp), as.integer(length(dp)),
+ as.double(w), as.double(S), as.double(lu), q=matrix(0, N, N))$q
}
-lossrecovdist <- function(defaultprob, prepayprob, S, lu, useC=TRUE){
+lossrecovdist <- function(defaultprob, prepayprob, w, S, lu, useC=TRUE){
if(prepayprob==0){
if(useC){
- L <- lossdistribC(defaultprob, S, lu)
- R <- lossdistribC(defaultprob, 1-S, lu)
+ L <- lossdistribC(defaultprob, w, S, lu)
+ R <- lossdistribC(defaultprob, w, 1-S, lu)
}else{
- L <- lossdistrib2(defaultprob, S, lu)
- R <- lossdistrib2(defaultprob, 1-S, lu)
+ L <- lossdistrib2(defaultprob, w, S, lu)
+ R <- lossdistrib2(defaultprob, w, 1-S, lu)
}
}else{
- u <- prepayprob/defaultprob
- L <- lossdistribC(defaultprob+prepayprob, S/(1+u), lu)
- R <- lossdistribC(defaultprob+prepayprob, (u+1-S)/(1+u), lu)
+ L <- lossdistribC(defaultprob, w, S, lu)
+ R <- recovdistC(defaultprob, prepayprob, w, S, lu)
}
return(list(L=L, R=R))
}
|