function to compute the average recovery distribution

1 files changed, 41 insertions, 0 deletions

diff --git a/src/lossdistrib.c b/src/lossdistrib.c
index b1b78c4..d8257bd 100644
--- a/src/lossdistrib.c
+++ b/src/lossdistrib.c
@@ -260,6 +260,47 @@ void lossdistrib_joint(const double *p, const double* pp, const int *np,
     free(qtemp);
 }
 
+void joint_default_averagerecov_distrib(const double *p,const int *np,
+                                        const double *S, const 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
+       S vector of severities (should be same length as p)
+       N number of ticks in the grid
+       q the joint probability distribution */
+
+    double lu = 1./(*N-1);
+    q[0] = 1;
+
+    double* newrecov = malloc(*N * sizeof(double));
+    unsigned int* index = malloc(*N * sizeof(unsigned int));
+    double* weights = malloc(*N * sizeof(double));
+    q[0] = 1;
+    double temp;
+    for(int i = 0; i < (*np); i++) {
+        for(int j = i + 1; j > 0; j--) {
+            for(int n = 0; n < *N; n++) {
+                newrecov[n] = ( (j-1) * n * lu + (1-S[i])) / (double)j;
+                weights[n] = modf(newrecov[n] * (*N-1), &temp);
+                index[n] = (unsigned int) temp;
+            }
+            if(j < i+1) {
+                cblas_dscal(*N, 1-p[i], q+(*N)*j, 1);
+            }
+            for(int k = 0; k < *N; k++) {
+                q[j*(*N)+index[k]+1] += weights[k] * p[i] * q[(j-1)*(*N)+k];
+                q[j*(*N)+index[k]] += (1-weights[k]) * p[i] * q[(j-1)*(*N)+k];
+            }
+        }
+        cblas_dscal(*N, 1-p[i], q, 1);
+    }
+    free(newrecov);
+    free(index);
+}
+
+
 void recovdist(const double *dp, const double *pp, const int *n, const double *w,
                const double *S, const int *N, double *q) {
     /* recursive algorithm with first order correction for computing