diff options
| -rw-r--r-- | lossdistrib.c | 140 | ||||
| -rw-r--r-- | tranche_functions.R | 6 |
2 files changed, 99 insertions, 47 deletions
diff --git a/lossdistrib.c b/lossdistrib.c index 511d2169..2d237e22 100644 --- a/lossdistrib.c +++ b/lossdistrib.c @@ -1,7 +1,7 @@ #include <R.h>
#include <Rmath.h>
-void lossdistrib(double *p, int *np, double *S, double *lu, double *q) {
+void lossdistrib(double *p, int *np, double *w, double *S, double *lu, double *q) {
/* recursive algorithm with first order correction for computing
the loss distribution.
p vector of default probabilities
@@ -19,7 +19,7 @@ void lossdistrib(double *p, int *np, double *S, double *lu, double *q) { q2 = malloc( N * sizeof(double));
q[0] = 1;
for(i=0; i<(*np); i++){
- d = S[i]/(*np * (*lu));
+ d = S[i] * w[i]/ (*lu);
d1 = floor(d);
d2 = ceil(d);
p1 = p[i] * (d2-d);
@@ -51,16 +51,26 @@ void lossdistrib(double *p, int *np, double *S, double *lu, double *q) { free(q2);
}
-void lossdistrib2( double *p, int *np, double *S, double *lu, double *q) {
+inline void aux2(double *q, double *qtemp, double p, int N, int i, int j){
+ int m, n;
+ for(m=i; m<N; m++){
+ for(n=j; n<N; n++){
+ qtemp[m+N*n]+= p * q[(m-i)+N*(n-j)];
+ }
+ }
+}
+
+void lossdistrib_joint( double *p, int *np, double *w, double *S, double *lu, 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
q the joint probability distribution */
- int i, j, k, l, m, n, N;
+ int i, j, k, l, N;
double x, y, sum;
double alpha1, alpha2, beta1, beta2;
double w1, w2, w3, w4;
@@ -70,8 +80,8 @@ void lossdistrib2( double *p, int *np, double *S, double *lu, double *q) { qtemp = malloc( N * N * sizeof(double));
q[0] = 1;
for(k=0; k<(*np); k++){
- x = S[k]/(*np * (*lu));
- y = (1-S[k])/(*np * (*lu));
+ x = S[k] * w[k]/(*lu);
+ y = (1-S[k]) * w[k] /(*lu);
i = floor(x);
j = floor(y);
alpha1 = i + 1 - x;
@@ -86,50 +96,32 @@ void lossdistrib2( double *p, int *np, double *S, double *lu, double *q) { for(l=0; l < N*N; l++){
qtemp[l] = 0;
}
- for(m=i; m<N; m++){
- for(n=j; n<N; n++){
- qtemp[m+N*n]+= w1 * p[k]* q[(m-i)+ N*(n-j)];
- }
- }
- for(m=i; m<N; m++){
- for(n=j+1; n<N; n++){
- qtemp[m+N*n]+= w2 * p[k]* q[(m-i)+ N*(n-j-1)];
- }
- }
- for(m=i+1; m<N; m++){
- for(n=j+1; n<N; n++){
- qtemp[m+N*n]+= w3 * p[k]* q[(m-i-1)+ N*(n-j-1)];
- }
- }
- for(m=i+1; m<N; m++){
- for(n=j; n<N; n++){
- qtemp[m+N*n]+= w4 * p[k]* q[(m-i-1)+ N*(n-j)];
- }
- }
- for(m=0; m<N; m++){
- for(n=0; n<N;n++){
- q[m+N*n] = (1-p[k]) * q[m+N*n]+ qtemp[m+N*n];
- }
+ 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(l=0; l<N*N; l++){
- sum+=q[l];
+ for(k=0; k<N*N; k++){
+ sum+=q[k];
}
- q[N * N -1]+=1-sum;
+ q[N*N-1]+= 1 - sum;
free(qtemp);
}
-void aux(double *q, double *qtemp, double p, int N, int d) {
+inline 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];
+ qtemp[j] += p * q[j-d];
};
}
-void lossdistrib3(double *dp, double *pp, int *n, double *S, double *w, double *lu, double *q) {
+void recovdist(double *dp, double *pp, int *n, double *w, double *S, double *lu, double *q) {
/* recursive algorithm with first order correction for computing
the loss distribution.
dp vector of default probabilities
@@ -148,12 +140,12 @@ void lossdistrib3(double *dp, double *pp, int *n, double *S, double *w, double * 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);
- d2u = d1u + 1;
+ d2u = d2l + 1;
dp1 = dp[i] * (d1u-d1);
dp2 = dp[i] - dp1;
pp1 = pp[i] * (d2u -d2);
@@ -175,7 +167,54 @@ void lossdistrib3(double *dp, double *pp, int *n, double *S, double *w, double * sum += q[j];
}
q[N-1] += 1-sum;
- free(q);
+ 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;
+ double alpha1, alpha2, beta1, beta2;
+ double w1, w2, w3, w4;
+ double *qtemp;
+
+ N = ceil(1./(*lu) + 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);
+ i = floor(x);
+ j1 = floor(y1);
+ j2 = floor(y2);
+ alpha1 = i + 1 - x;
+ alpha2 = 1 - alpha1;
+ beta1 = j1 + 1 - y1;
+ beta2 = 1 - beta1;
+ w1 = alpha1 * beta1;
+ w2 = alpha1 * beta2;
+ w3 = alpha2 * beta2;
+ w4 = alpha2 * beta1;
+ /* init qtemp to 0 */
+ 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++){
+ q[l] = qtemp[l] + (1-dp[k]-pp[k]) * q[l];
+ }
+ }
+ /* correction for weight loss */
+ sum = 0;
+ for(k=0; k<N*N; k++){
+ sum+=q[k];
+ }
+ q[N*N-1]+= 1 - sum;
free(qtemp);
}
@@ -195,8 +234,21 @@ void addandmultiply(double *X, double alpha, double *Y, int n) { }
}
-void BClossdist(double *SurvProb, int *dim1, int *dim2, 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 *recov, double *issuerweights, double *Z, double *w, int *n, double *rho, double *lu, double *L, double *R) {
+ /*
+ computes the loss and recovery distribution over time with a flat gaussiancorrelation
+ inputs:
+ Survprob: matrix of size dim1 x dim2. dim1 is the number of issuers and dim2 number of time steps
+ recov: vector of recoveries (length dim1)
+ issuerweights: vector of issuer weights (length dim2)
+ Z: vector of factor values (length n)
+ w: vector of factor weights (length n)
+ rho: correlation beta
+ lu: loss unit
+ outputs:
+ L: matrix of size (1/lu+1, dim2)
+ R: matrix of size (1/lu+1, dim2)
+ */
int t, i, j, N;
double g;
double *gshocked, *Rshocked, *Sshocked, *Lw, *Rw;
@@ -221,8 +273,8 @@ void BClossdist(double *SurvProb, int *dim1, int *dim2, double *recov, double *Z Lw[j]=0;
Rw[j]=0;
}
- lossdistrib(gshocked, dim1, Sshocked, lu, Lw);
- lossdistrib(gshocked, dim1, Rshocked, lu, Rw);
+ 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);
}
diff --git a/tranche_functions.R b/tranche_functions.R index 1b20fdc7..42597c68 100644 --- a/tranche_functions.R +++ b/tranche_functions.R @@ -156,7 +156,7 @@ lossdistribC <- function(p, w, S, lu){ recovdistC <- function(dp, pp, w, S, lu){
## C version of recovdist
- dyn.load("lossdistrib2.dll")
+ 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
}
@@ -164,14 +164,14 @@ recovdistC <- function(dp, pp, w, S, lu){ lossdistribC.joint <- function(p, w, S, lu){
## C version of lossdistrib.joint, roughly 20 times faster
N <- ceiling(1/lu+1)
- dyn.load("lossdistrib2.dll")
+ 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
}
lossdistribprepayC.joint <- function(dp, pp, w, S, lu){
## C version of lossdistribprepay.joint
- dyn.load("lossdistrib2.dll")
+ 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
|