aboutsummaryrefslogtreecommitdiffstats
path: root/lossdistrib.c
diff options
context:
space:
mode:
Diffstat (limited to 'lossdistrib.c')
-rw-r--r--lossdistrib.c140
1 files changed, 96 insertions, 44 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);
}
generated by cgit v1.2.3-70-g09d2 (git 2.52.0) at 2026-01-22 18:21:14 +0000