aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--lossdistrib.c195
-rw-r--r--tranche_functions.R184
2 files changed, 249 insertions, 130 deletions
diff --git a/lossdistrib.c b/lossdistrib.c
index 7e5a8908..e6eef029 100644
--- a/lossdistrib.c
+++ b/lossdistrib.c
@@ -1,25 +1,27 @@
#include <R.h>
#include <Rmath.h>
-void lossdistrib(double *p, int *np, double *w, double *S, double *lu, double *q) {
+#define MIN(x, y) (((x) < (y)) ? (x) : (y))
+
+void lossdistrib(double *p, int *np, double *w, double *S, int *N, double *q) {
/* recursive algorithm with first order correction for computing
the loss distribution.
p vector of default probabilities
np length of p
S vector of severities (should be same length as p)
- lu loss unit
+ N number of ticks in the grid
q the loss distribution */
- int i, j, N, d1, d2;
- double d, p1, p2, sum;
+ int i, j, d1, d2;
+ double lu, d, p1, p2, sum;
double *q1, *q2;
- N = ceil(1./(*lu) + 1);
- q1 = malloc( N * sizeof(double));
- q2 = malloc( N * sizeof(double));
+ lu = 1./(*N-1);
+ q1 = malloc( *N * sizeof(double));
+ q2 = malloc( *N * sizeof(double));
q[0] = 1;
for(i=0; i<(*np); i++){
- d = S[i] * w[i]/ (*lu);
+ d = S[i] * w[i]/ lu;
d1 = floor(d);
d2 = ceil(d);
p1 = p[i] * (d2-d);
@@ -28,25 +30,67 @@ void lossdistrib(double *p, int *np, double *w, double *S, double *lu, double *q
for(j=0; j<d1; j++){
q1[j]=0;
};
- for(j=0; j<N-d1; j++){
+ for(j=0; j<*N-d1; j++){
q1[d1+j] = p1*q[j];
};
for(j=0; j<d2; j++){
q2[j]=0;
};
- for(j=0; j<N-d2; j++){
+ for(j=0; j<*N-d2; j++){
q2[d2+j] = p2*q[j];
};
- for(j=0; j<N; j++){
+ for(j=0; j<*N; j++){
q[j] = q1[j] + q2[j] + (1-p[i]) * q[j];
};
}
/* correction for weight loss */
- sum = 0;
- for(j=0; j<N; j++){
- sum += q[j];
+ /* sum = 0; */
+ /* for(j=0; j<*N; j++){ */
+ /* sum += q[j]; */
+ /* } */
+ /* q[*N-1] += 1-sum; */
+ free(q1);
+ free(q2);
+}
+
+void lossdistrib_fast(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
+ np length of p
+ S vector of severities (should be same length as p)
+ N number of ticks in the grid
+ T where to truncate the distribution
+ q the loss distribution */
+
+ int i, j, d1, d2, M;
+ double lu, d, p1, p2, sum;
+ double *q1, *q2;
+
+ lu = 1./(*N-1);
+ q1 = calloc( *T, sizeof(double));
+ q2 = calloc( *T, sizeof(double));
+ q[0] = 1;
+ M = 1;
+ for(i=0; i<(*np); i++){
+ d = S[i] * w[i]/ lu;
+ d1 = floor(d);
+ d2 = ceil(d);
+ p1 = p[i] * (d2-d);
+ p2 = p[i] - p1;
+ for(j=0; j < MIN(M, *T); j++){
+ q1[j] = p1 * q[j];
+ q2[j] = p2 * q[j];
+ q[j] = (1-p[i]) * q[j];
+ }
+ for(j=0; j < MIN(M, *T-d1); j++){
+ q[d1+j] += q1[j];
+ };
+ for(j=0; j < MIN(M, *T-d2); j++){
+ q[d2+j] += q2[j];
+ };
+ M += d2;
}
- q[N-1] += 1-sum;
free(q1);
free(q2);
}
@@ -60,28 +104,28 @@ inline void aux2(double *q, double *qtemp, double p, int N, int i, int j){
}
}
-void lossdistrib_joint( double *p, int *np, double *w, double *S, double *lu, double *q) {
+void lossdistrib_joint( double *p, int *np, double *w, double *S, int *N, double *q) {
/* recursive algorithm with first order correction
computes jointly the loss and recovery distribution
p vector of default probabilities
np length of p
w vector of issuer weights (length np)
S vector of severities (should be same length as p)
- lu loss unit
+ N number of ticks in the grid
q the joint probability distribution */
- int i, j, k, l, N;
- double x, y, sum;
+ int i, j, k, l;
+ double lu, x, y, sum;
double alpha1, alpha2, beta1, beta2;
double w1, w2, w3, w4;
double *qtemp;
- N = ceil(1./(*lu) + 1);
- qtemp = malloc( N * N * sizeof(double));
+ lu = 1/(*N-1);
+ qtemp = malloc( (*N) * (*N) * sizeof(double));
q[0] = 1;
for(k=0; k<(*np); k++){
- x = S[k] * w[k]/(*lu);
- y = (1-S[k]) * w[k] /(*lu);
+ x = S[k] * w[k]/ lu;
+ y = (1-S[k]) * w[k] / lu;
i = floor(x);
j = floor(y);
alpha1 = i + 1 - x;
@@ -93,23 +137,23 @@ void lossdistrib_joint( double *p, int *np, double *w, double *S, double *lu, do
w3 = alpha2 * beta2;
w4 = alpha2 * beta1;
/* init qtemp to 0 */
- for(l=0; l < N*N; l++){
+ for(l=0; l < (*N)*(*N); l++){
qtemp[l] = 0;
}
- aux2(q, qtemp, w1 * p[k], N, i, j);
- aux2(q, qtemp, w2 * p[k], N, i, j+1);
- aux2(q, qtemp, w3 * p[k], N, i+1, j+1);
- aux2(q, qtemp, w4 * p[k], N, i+1, j);
- for(l=0; l < N * N; l++){
+ aux2(q, qtemp, w1 * p[k], *N, i, j);
+ aux2(q, qtemp, w2 * p[k], *N, i, j+1);
+ aux2(q, qtemp, w3 * p[k], *N, i+1, j+1);
+ aux2(q, qtemp, w4 * p[k], *N, i+1, j);
+ for(l=0; l < (*N) * (*N); l++){
q[l] = qtemp[l] + (1-p[k]) * q[l];
}
}
/* correction for weight loss */
sum = 0;
- for(k=0; k<N*N; k++){
+ for(k=0; k < (*N)*(*N); k++){
sum+=q[k];
}
- q[N*N-1]+= 1 - sum;
+ q[(*N)*(*N)-1]+= 1 - sum;
free(qtemp);
}
@@ -121,7 +165,7 @@ inline void aux(double *q, double *qtemp, double p, int N, int d) {
};
}
-void recovdist(double *dp, double *pp, int *n, double *w, double *S, double *lu, double *q) {
+void recovdist(double *dp, double *pp, int *n, double *w, double *S, int *N, double *q) {
/* recursive algorithm with first order correction for computing
the loss distribution.
dp vector of default probabilities
@@ -129,19 +173,19 @@ void recovdist(double *dp, double *pp, int *n, double *w, double *S, double *lu,
n length of p
S vector of severities (should be same length as p)
w vector of weights
- lu loss unit
+ N number of ticks in the grid
q the loss distribution */
- int i, j, N, d1l, d1u, d2l, d2u;
- double d1, d2, dp1, dp2, pp1, pp2, sum;
+ int i, j, d1l, d1u, d2l, d2u;
+ double lu, d1, d2, dp1, dp2, pp1, pp2, sum;
double *qtemp;
- N = ceil(1./(*lu) + 1);
- qtemp = malloc( N * sizeof(double));
+ lu = 1./(*N - 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);
+ d1 = w[i] * (1-S[i]) /lu;
+ d2 = w[i]/lu;
d1l = floor(d1);
d1u = d1l + 1;
d2l = floor(d2);
@@ -150,40 +194,40 @@ void recovdist(double *dp, double *pp, int *n, double *w, double *S, double *lu,
dp2 = dp[i] - dp1;
pp1 = pp[i] * (d2u -d2);
pp2 = pp[i] - pp1;
- for(j=0; j<N; j++){
+ 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++){
+ 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++){
+ for(j=0; j<(*N); j++){
sum += q[j];
}
- q[N-1] += 1-sum;
+ q[*N-1] += 1-sum;
free(qtemp);
}
-void lossdistrib_prepay_joint(double *dp, double *pp, int *n, double *w, double *S, double *lu, double *q) {
- int i, j1, j2, k, l, N;
- double x, y1, y2, sum;
+void lossdistrib_prepay_joint(double *dp, double *pp, int *n, double *w, double *S, int *N, double *q) {
+ int i, j1, j2, k, l;
+ double lu, x, y1, y2, sum;
double alpha1, alpha2, beta1, beta2;
double w1, w2, w3, w4;
double *qtemp;
- N = ceil(1./(*lu) + 1);
- qtemp = malloc( N * N * sizeof(double));
+ lu = 1./(*N-1);
+ qtemp = malloc( (*N) * (*N) * sizeof(double));
q[0] = 1;
for(k=0; k<(*n); k++){
- x = S[k] * w[k]/(*lu);
- y1 = (1-S[k]) * w[k]/(*lu);
- y2 = w[k]/(*lu);
+ 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);
@@ -196,25 +240,25 @@ void lossdistrib_prepay_joint(double *dp, double *pp, int *n, double *w, double
w3 = alpha2 * beta2;
w4 = alpha2 * beta1;
/* init qtemp to 0 */
- for(l=0; l < N*N; l++){
+ for(l=0; l < (*N)*(*N); l++){
qtemp[l] = 0;
}
- aux2(q, qtemp, w1 * dp[k], N, i, j1);
- aux2(q, qtemp, w2 * dp[k], N, i, j1+1);
- aux2(q, qtemp, w3 * dp[k], N, i+1, j1+1);
- aux2(q, qtemp, w4 * dp[k], N, i+1, j1);
- aux2(q, qtemp, pp[k]*(j2+1-y2), N, 0, j2);
- aux2(q, qtemp, pp[k]*(y2-j2), N, 0, j2+1);
- for(l=0; l < N * N; l++){
+ aux2(q, qtemp, w1 * dp[k], *N, i, j1);
+ aux2(q, qtemp, w2 * dp[k], *N, i, j1+1);
+ aux2(q, qtemp, w3 * dp[k], *N, i+1, j1+1);
+ aux2(q, qtemp, w4 * dp[k], *N, i+1, j1);
+ aux2(q, qtemp, pp[k]*(j2+1-y2), *N, 0, j2);
+ aux2(q, qtemp, pp[k]*(y2-j2), *N, 0, j2+1);
+ for(l=0; l < (*N)*(*N); l++){
q[l] = qtemp[l] + (1-dp[k]-pp[k]) * q[l];
}
}
/* correction for weight loss */
sum = 0;
- for(k=0; k<N*N; k++){
+ for(k=0; k<(*N)*(*N); k++){
sum+=q[k];
}
- q[N*N-1]+= 1 - sum;
+ q[(*N)*(*N)-1]+= 1 - sum;
free(qtemp);
}
@@ -234,7 +278,7 @@ 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, double *lu, 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, double *L, double *R) {
/*
computes the loss and recovery distribution over time with a flat gaussiancorrelation
inputs:
@@ -244,21 +288,20 @@ void BClossdist(double *SurvProb, int *dim1, int *dim2, double *issuerweights, d
Z: vector of factor values (length n)
w: vector of factor weights (length n)
rho: correlation beta
- lu: loss unit
+ N: number of ticks in the grid
outputs:
L: matrix of size (1/lu+1, dim2)
R: matrix of size (1/lu+1, dim2)
*/
- int t, i, j, N;
+ int t, i, j;
double g;
double *gshocked, *Rshocked, *Sshocked, *Lw, *Rw;
gshocked = calloc((*dim1), sizeof(double));
Rshocked = calloc((*dim1), sizeof(double));
Sshocked = calloc((*dim1), sizeof(double));
- N = ceil(1. / (*lu) + 1);
- Lw = malloc(N * sizeof(double));
- Rw = malloc(N * sizeof(double));
+ Lw = malloc((*N) * sizeof(double));
+ Rw = malloc((*N) * sizeof(double));
for(t=0; t < (*dim2); t++) {
for(i=0; i < *n; i++){
@@ -269,14 +312,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<*N; j++){
Lw[j]=0;
Rw[j]=0;
}
- lossdistrib(gshocked, dim1, issuerweights, Sshocked, lu, Lw);
- lossdistrib(gshocked, dim1, issuerweights, Rshocked, lu, Rw);
- addandmultiply(Lw, w[i], L + t * N, N);
- addandmultiply(Rw, w[i], R + t * N, N);
+ 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);
}
}
free(gshocked);
diff --git a/tranche_functions.R b/tranche_functions.R
index a7262ea0..69b14511 100644
--- a/tranche_functions.R
+++ b/tranche_functions.R
@@ -6,7 +6,7 @@ lossdistrib <- function(p){
q <- rep(0, (n+1))
q[1] <- 1
for(i in 1:n){
- q[-(n+1)] <- p[i]*q[-(n+1)]+(1-p[i])*q[-1]
+ q[-1] <- p[i]*q[-(n+1)]+(1-p[i])*q[-1]
q[1] <- (1-p[i])*q[1]
}
return(q)
@@ -23,16 +23,15 @@ lossdistrib.fft <- function(p){
return(1/(n+1)*Re(fft(v)))
}
-lossdistrib2 <- function(p, w, S, lu){
+lossdistrib2 <- function(p, w, S, N){
#recursive algorithm with first order correction
#p vector of default probabilities
#w vector of weigths
#S vector of severities
- #lu loss unit
+ #N number of ticks in the grid
n <- length(p)
- N <- 1/lu + 1
- eps <- 1e-10
- q <- rep(0, N+eps)
+ lu <- 1/(N-1)
+ q <- rep(0, N)
q[1] <- 1
for(i in 1:n){
d <- S[i] * w[i] / lu
@@ -44,11 +43,75 @@ lossdistrib2 <- function(p, w, S, lu){
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)
+ ## 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)
return(q)
}
-recovdist <- function(dp, pp, w, S, lu){
+lossdistrib2.fast.trunc <- function(p, w, S, N, truncated=1){
+ ## 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
+ ## (useful for tranche computations)
+
+ ## this is actually slower than lossdistrib2. But in C this is
+ ## twice as fast.
+ n <- length(p)
+ lu <- 1/(N-1)
+ q <- rep(0, truncated)
+ 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: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
+ 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
## where X_i = 0 w.p 1-dp[i]-pp[i]
@@ -58,9 +121,9 @@ recovdist <- function(dp, pp, w, S, lu){
## the pair of values floor(v/lu) and ceiling(v/lu) so that
## X_i has four non zero values
n <- length(dp)
- N <- 1/lu + 1
- q <- rep(0, ceiling(N))
+ q <- rep(0, N)
q[1] <- 1
+ lu <- 1/(N-1)
for(i in 1:n){
d1 <- w[i]*(1-S[i])/lu
d1l <- floor(d1)
@@ -82,20 +145,20 @@ recovdist <- function(dp, pp, w, S, lu){
return(q)
}
-lossdistrib.joint <- function(p, w, S, lu){
+lossdistrib.joint <- function(p, w, S, N){
## recursive algorithm with first order correction
## to compute the joint probability distribution of the loss and recovery
## inputs:
## p: vector of default probabilities
## w: vector of issuer weights
## S: vector of severities
- ## lu: loss unit
+ ## N: number of tick sizes on the grid
## output:
## q: matrix of joint loss, recovery probability
## colSums(q) is the recovery distribution marginal
## rowSums(q) is the loss distribution marginal
n <- length(p)
- N <- ceiling(1/lu + 1)
+ lu <- 1/(N-1)
q <- matrix(0, N, N)
q[1,1] <- 1
for(k in 1:n){
@@ -116,7 +179,7 @@ lossdistrib.joint <- function(p, w, S, lu){
return(q)
}
-lossdistribprepay.joint <- function(dp, pp, w, S, lu){
+lossdistribprepay.joint <- function(dp, pp, w, S, N){
## recursive algorithm with first order correction
## to compute the joint probability distribition of the loss and recovery
## inputs:
@@ -124,13 +187,13 @@ lossdistribprepay.joint <- function(dp, pp, w, S, lu){
## pp: vector of prepay probabilities
## w: vector of issuer weights
## S: vector of severities
- ## lu: loss unit
+ ## N: number of tick sizes on the grid
## outputs
## q: matrix of joint loss and recovery probability
## colSums(q) is the recovery distribution marginal
## rowSums(q) is the loss distribution marginal
n <- length(dp)
- N <- ceiling(1/lu + 1)
+ lu <- 1/(N-1)
q <- matrix(0, N, N)
q[1,1] <- 1
for(k in 1:n){
@@ -156,60 +219,64 @@ lossdistribprepay.joint <- function(dp, pp, w, S, lu){
return(q)
}
-
-lossdistribC <- function(p, w, S, lu){
+lossdistribC <- function(p, w, S, N){
## C version of lossdistrib2, roughly 50 times faster
dyn.load("lossdistrib.dll")
.C("lossdistrib", as.double(p), as.integer(length(p)),
- as.double(w), as.double(S), as.double(lu), q = double(1/lu+1))$q
+ as.double(w), as.double(S), as.integer(N), q = double(N))$q
}
-recovdistC <- function(dp, pp, w, S, lu){
+lossdistrib.fastC <- 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)),
+ as.double(w), as.double(S), as.integer(N), as.integer(T), q = double(T))$q
+}
+
+recovdistC <- function(dp, pp, w, S, N){
## C version of recovdist
dyn.load("lossdistrib.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
+ as.double(w), as.double(S), as.integer(N), q = double(N))$q
}
-lossdistribC.joint <- function(p, w, S, lu){
+lossdistribC.joint <- function(p, w, S, N){
## C version of lossdistrib.joint, roughly 20 times faster
- N <- ceiling(1/lu+1)
dyn.load("lossdistrib.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
+ as.double(S), as.integer(N), q = matrix(0, N, N))$q
}
-lossdistribprepayC.joint <- function(dp, pp, w, S, lu){
+lossdistribprepayC.joint <- function(dp, pp, w, S, N){
## C version of lossdistribprepay.joint
dyn.load("lossdistrib.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
+ as.double(w), as.double(S), as.integer(N), q=matrix(0, N, N))$q
}
-lossrecovdist <- function(defaultprob, prepayprob, w, S, lu, useC=TRUE){
+lossrecovdist <- function(defaultprob, prepayprob, w, S, N, useC=TRUE){
if(all(!prepayprob)){
if(useC){
- L <- lossdistribC(defaultprob, w, S, lu)
- R <- lossdistribC(defaultprob, w, 1-S, lu)
+ L <- lossdistribC(defaultprob, w, S, N)
+ R <- lossdistribC(defaultprob, w, 1-S, N)
}else{
- L <- lossdistrib2(defaultprob, w, S, lu)
- R <- lossdistrib2(defaultprob, w, 1-S, lu)
+ L <- lossdistrib2(defaultprob, w, S, N)
+ R <- lossdistrib2(defaultprob, w, 1-S, N)
}
}else{
- L <- lossdistribC(defaultprob, w, S, lu)
- R <- recovdistC(defaultprob, prepayprob, w, S, lu)
+ L <- lossdistribC(defaultprob, w, S, N)
+ R <- recovdistC(defaultprob, prepayprob, w, S, N)
}
return(list(L=L, R=R))
}
-lossrecovdist.term <- function(defaultprob, prepayprob, w, S, lu, useC=TRUE){
+lossrecovdist.term <- function(defaultprob, prepayprob, w, S, N, useC=TRUE){
## computes the loss and recovery distribution over time
- L <- array(0, dim=c(1/lu+1, ncol(defaultprob)))
- R <- array(0, dim=c(1/lu+1, ncol(defaultprob)))
+ L <- array(0, dim=c(N, ncol(defaultprob)))
+ R <- array(0, dim=c(N, ncol(defaultprob)))
if(all(!prepayprob)){
for(t in 1:ncol(defaultprob)){
- temp <- lossrecovdist(defaultprob[,t], 0, w, S[,t], lu, useC)
+ temp <- lossrecovdist(defaultprob[,t], 0, w, S[,t], N, useC)
L[,t] <- temp$L
R[,t] <- temp$R
}
@@ -324,21 +391,20 @@ tranche.bpvec <- function(K, L, R, cs){
return( r )
}
-BClossdist <- function(SurvProb, issuerweights, recov, rho, lu, n.int=100){
+BClossdist <- function(SurvProb, issuerweights, recov, rho, N, n.int=100){
quadrature <- gauss.quad.prob(n.int, "normal")
Z <- quadrature$nodes
w <- quadrature$weights
- n <- as.integer(1/lu+1)
- LZ <- matrix(0, n, n.int)
- RZ <- matrix(0, n, n.int)
- L <- matrix(0, n, ncol(SurvProb))
- R <- matrix(0, n, ncol(SurvProb))
+ LZ <- matrix(0, N, n.int)
+ RZ <- matrix(0, N, n.int)
+ L <- matrix(0, N, ncol(SurvProb))
+ R <- matrix(0, N, ncol(SurvProb))
for(t in 1:ncol(SurvProb)){
g <- 1 - SurvProb[, t]
for(i in 1:length(Z)){
g.shocked <- shockprob(g, rho, Z[i])
S.shocked <- shockseverity(1-recov, 1, Z[i], rho, g)
- temp <- lossrecovdist(g.shocked, 0, issuerweights, S.shocked, lu)
+ temp <- lossrecovdist(g.shocked, 0, issuerweights, S.shocked, N)
LZ[,i] <- temp$L
RZ[,i] <- temp$R
}
@@ -348,19 +414,18 @@ BClossdist <- function(SurvProb, issuerweights, recov, rho, lu, n.int=100){
list(L=L, R=R)
}
-BClossdistC <- function(SurvProb, issuerweights, recov, rho, lu, n.int=100){
+BClossdistC <- function(SurvProb, issuerweights, recov, rho, N, n.int=100){
dyn.load("lossdistrib.dll")
quadrature <- gauss.quad.prob(n.int, "normal")
Z <- quadrature$nodes
w <- quadrature$weights
- N <- as.integer(1/lu+1)
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.double(lu), L=L, R=R)
+ 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)
return(list(L=r$L,R=r$R))
}
-BCtranche.pv <- function(portfolio, index, coupon, K1, K2, rho1, rho2, lu=0.01){
+BCtranche.pv <- function(portfolio, index, coupon, K1, K2, rho1, rho2, N=101){
## computes the protection leg, couponleg, and bond price of a tranche
## in the base correlation setting
if(K1==0){
@@ -374,9 +439,9 @@ BCtranche.pv <- function(portfolio, index, coupon, K1, K2, rho1, rho2, lu=0.01){
issuerweights <- rep(1/length(portfolio), length(portfolio))
K <- adjust.attachments(c(K1,K2), index$loss, index$factor)
dK <- K[2] - K[1]
- dist2 <- BClossdistC(SurvProb, issuerweights, recov, rho2, lu)
+ dist2 <- BClossdistC(SurvProb, issuerweights, recov, rho2, N)
if(rho1!=0){
- dist1 <- BClossdistC(SurvProb, issuerweights, recov, rho1, lu)
+ dist1 <- BClossdistC(SurvProb, issuerweights, recov, rho1, N)
}
cl2 <- tranche.cl(dist2$L, dist2$R, cs, 0, K[2])
cl1 <- tranche.cl(dist1$L, dist1$R, cs, 0, K[1])
@@ -386,7 +451,7 @@ BCtranche.pv <- function(portfolio, index, coupon, K1, K2, rho1, rho2, lu=0.01){
bp=100*(1+(pl2-pl1+cl2-cl1)/dK)))
}
-BCtranche.delta <- function(portfolio, index, coupon, K1, K2, rho1, rho2, lu=0.01){
+BCtranche.delta <- function(portfolio, index, coupon, K1, K2, rho1, rho2, N=101){
## 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
@@ -414,11 +479,11 @@ BCtranche.delta <- function(portfolio, index, coupon, K1, K2, rho1, rho2, lu=0.0
return( delta )
}
-BCstrikes <- function(portfolio, index, coupon, K, rho) {
+BCstrikes <- function(portfolio, index, coupon, K, rho, N=101) {
## computes the strikes as a percentage of expected loss
EL <- c()
for(i in 2:length(K)){
- EL <- c(EL, -BCtranche.pv(portfolio, index, coupon, K[i-1], K[i], rho[i-1], rho[i], lu)$pl)
+ EL <- c(EL, -BCtranche.pv(portfolio, index, coupon, K[i-1], K[i], rho[i-1], rho[i], N)$pl)
}
Kmodified <- adjust.attachments(K, index$loss, index$factor)
return(cumsum(EL*diff(Kmodified))/sum(EL*diff(Kmodified)))
@@ -431,3 +496,14 @@ 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))
generated by cgit v1.2.3-70-g09d2 (git 2.52.0) at 2026-01-26 08:36:34 +0000