new function using truncated dists

1 files changed, 103 insertions, 16 deletions

diff --git a/src/lossdistrib.c b/src/lossdistrib.c
index 51e273b..9330cfb 100644
--- a/src/lossdistrib.c
+++ b/src/lossdistrib.c
@@ -6,15 +6,6 @@
 
 #define MIN(x, y) (((x) < (y)) ? (x) : (y))
 
-extern int dgemv_(char* trans, int *m, int *n, double* alpha, double* A, int* lda,
-                  double* x, int* incx, double* beta, double* y, int* incy);
-extern double ddot_(int* n, double* dx, int* incx, double* dy, int* incy);
-extern int dscal_(int* n, double* da, double* dx, int* incx);
-extern int daxpy_(int* n, double* da, double* dx, int* incx, double* dy, int* incy);
-extern int dstev_(char* JOBZ, int* n, double* D, double* E, double* Z, int* ldz, double* WORK, int* INFO);
-
-extern void openblas_set_num_threads(int);
-
 void GHquad(int* n, double* Z, double* w) {
     // Setup for eigenvalue computations
     char JOBZ   = 'V'; // Compute eigenvalues & vectors
@@ -170,7 +161,7 @@ void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N,
        output:
        q the loss distribution */
 
-    int i, j, d1, d2, M;
+    int i, d1, d2, M;
     double lu, d, p1, p2;
     double *qtemp;
     int one = 1;
@@ -200,20 +191,27 @@ void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N,
     free(qtemp);
 }
 
+static inline void posK(int T, double K, double lu, double* val){
+    int i = 0;
+    for(i = 0; i < T; i++){
+        val[i] = K-lu*i;
+    }
+}
+
 void exp_trunc(double *p, int *np, double *w, double *S, int *N, double *K,
                double *r) {
     double lu;
     double *qtemp;
-    double lambda;
+
     lu = 1./(*N+1);
     int T = (int) floor((*K) * (*N))+1;
     int zero = 0;
+    int one = 1;
     qtemp = calloc( T, sizeof(double));
-    int i;
     lossdistrib_truncated(p, np, w, S, N, &T, &zero, qtemp);
-    for(i = 0; i < T; i++){
-        *r += (*K - lu*i) * qtemp[i];
-    }
+    double* val = malloc(T * sizeof(double));
+    posK(T, *K, lu, val);
+    *r = ddot_(&T, val, &one, qtemp, &one);
     free(qtemp);
 }
 
@@ -678,7 +676,7 @@ void stochasticrecov(double* R, double* Rtilde, double* Z, double* w, int* nZ,
     }else{
         ptemp = (1 - *R) / (1 - *Rtilde) * *porig;
         fitprob(Z, w, nZ, rho, &ptemp, &ptilde);
-#pragma omp parallel for
+        #pragma omp parallel for
         for(i = 0; i < *nZ; i++){
             q[i] = fabs(1 - (1 - *Rtilde) * exp( shockprob(ptilde, *rho, Z[i], 1) -
                                                  shockprob(*pmod, *rho, Z[i], 1)));
@@ -804,6 +802,95 @@ void BCloss_recov_dist(double *defaultprob, int *dim1, int *dim2,
     free(Rw);
 }
 
+void BCloss_recov_trunc(double *defaultprob, int *dim1, int *dim2,
+                        double *issuerweights, double *recov, double *Z, double *w,
+                        int *n, double *rho, int *N, double * K, int *defaultflag,
+                        double *ELt, double *ERt) {
+    /*
+      computes EL_t = E[(K-L_t)^+] and ER_t = E[(K-(1-R_t))^+]
+       the the put options on loss and recovery over time
+       with a flat gaussian correlation.
+      inputs:
+      defaultprob: matrix of size dim1 x dim2. dim1 is the number of issuers
+      and dim2 number of time steps
+      issuerweights: vector of issuer weights (length dim1)
+      recov: vector of recoveries (length dim1)
+      Z: vector of factor values (length n)
+      w: vector of factor weights (length n)
+      rho: correlation beta vector (length dim1)
+      N: number of ticks in the grid
+      K: the strike
+      defaultflag: if true, computes the default distribution
+      outputs:
+      ELt: vector of length dim2
+      ERt: vector of length dim2
+    */
+    int t, i, j;
+    double g, Ktilde;
+    int one = 1;
+    int T = (int) floor((*K) * (*N))+1;
+    int Ttilde;
+    double lu = 1./(*N+1);
+    double* valL = malloc( T * sizeof(double));
+    posK(T, *K, lu, valL);
+    double* EL = malloc( (*n) * sizeof(double));
+    double* ER = malloc( (*n) * sizeof(double));
+    double* Lw = malloc(T * (*n) * sizeof(double));
+    double alpha = 1;
+    double beta = 0;
+
+    for(t=0; t < (*dim2); t++) {
+        memset(Lw, 0, T * (*n) * sizeof(double));
+        #pragma omp parallel for private(j, g, Ktilde, Ttilde)
+        for(i=0; i < *n; i++){
+            double* Rw = NULL;
+            double* Rshocked = malloc((*dim1) * sizeof(double));
+            double* Sshocked = malloc((*dim1) * sizeof(double));
+            double* gshocked = malloc((*dim1) * sizeof(double));
+            double* gshockedbar = malloc((*dim1) * sizeof(double));
+            double* valR = NULL;
+
+            for(j=0; j < (*dim1); j++){
+                g = defaultprob[j + (*dim1) * t];
+                gshocked[j] = shockprob(g, rho[j], Z[i], 0);
+                Sshocked[j] = shockseverity(1-recov[j], Z[i], rho[j], g);
+                gshockedbar[j] = 1 - gshocked[j];
+                Rshocked[j] = 1 - Sshocked[j];
+            }
+
+            lossdistrib_truncated(gshocked, dim1, issuerweights, Sshocked,
+                                  N, &T, defaultflag, Lw + i * T);
+            ER[i] = 0;
+            Ktilde = *K - ddot_(dim1, issuerweights, &one, Sshocked, &one);
+            if(Ktilde > 0){
+                Ttilde = (int) floor(Ktilde * (*N))+1;
+                Rw = calloc(Ttilde, sizeof(double));
+                lossdistrib_truncated(gshockedbar, dim1, issuerweights, Rshocked,
+                                      N, &Ttilde, defaultflag, Rw);
+                valR = malloc(Ttilde * sizeof(double));
+                posK(Ttilde, Ktilde, lu, valR);
+                ER[i] = ddot_(&Ttilde, Rw, &one, valR, &one);
+            }
+            if(Rw != NULL){
+                free(Rw);
+            }
+            if(valR != NULL){
+                free(valR);
+            }
+            free(Rshocked);
+            free(Sshocked);
+            free(gshocked);
+            free(gshockedbar);
+        }
+        dgemv_("t", &T, n, &alpha, Lw, &T, valL, &one, &beta, EL, &one);
+        ELt[t] = ddot_(n, EL, &one, w, &one);
+        ERt[t] = ddot_(n, ER, &one, w, &one);
+    }
+    free(Lw);
+    free(EL);
+    free(ER);
+    free(valL);
+}
 
 void BCloss_dist(double *defaultprob, int *dim1, int *dim2,
                  double *issuerweights, double *recov, double *Z, double *w,