1 files changed, 58 insertions, 0 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));
 }