From 74fa05ce92c4d23d14cc6672ffefe3379a804bcd Mon Sep 17 00:00:00 2001
From: Guillaume Horel <guillaume.horel@serenitascapital.com>
Date: Tue, 10 Mar 2015 18:25:03 -0400
Subject: add an function for computing truncated expected loss

---
 R/distrib.R       | 18 ++++++++++++++++++
 src/lossdistrib.c | 19 ++++++++++++++++++-
 2 files changed, 36 insertions(+), 1 deletion(-)

diff --git a/R/distrib.R b/R/distrib.R
index 277b324..b2c028a 100644
--- a/R/distrib.R
+++ b/R/distrib.R
@@ -302,6 +302,24 @@ lossdistC.truncated <- function(p, w, S, N, T=N, defaultflag=FALSE){
        q = double(N))$q
 }
 
+exp.trunc <- function(p, w, S, N, K){
+    ## computes E[(K-L)^+]
+    r <- 0
+    .C("exp_trunc", as.double(p), as.integer(length(p)),
+       as.double(w), as.double(S), as.integer(N), as.double(K), res = r)$res
+}
+
+rec.trunc <- function(p, w, S, N, K){
+    ## computes E[(K-(1-R))^+] = E[(\tilde K- \bar R)]
+    ## where \tilde K = K-sum_i w_i S_i and \bar R=\sum_i w_i R_i (1-X_i)
+    Ktilde <- K-crossprod(w, S)
+    if(Ktilde < 0){
+        return( 0 )
+    }else{
+        return( exp.trunc(1-p, w, 1-S, N, Ktilde) )
+    }
+}
+
 recovdistC <- function(dp, pp, w, S, N){
     ## C version of recovdist
     .C("recovdist", as.double(dp), as.double(pp), as.integer(length(dp)),
diff --git a/src/lossdistrib.c b/src/lossdistrib.c
index d28b1b2..51e273b 100644
--- a/src/lossdistrib.c
+++ b/src/lossdistrib.c
@@ -200,6 +200,23 @@ void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N,
     free(qtemp);
 }
 
+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;
+    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];
+    }
+    free(qtemp);
+}
+
 void lossdistrib_joint(double *p, int *np, double *w, double *S, int *N, int *defaultflag, double *q) {
     /* recursive algorithm with first order correction
        computes jointly the loss and recovery distribution
@@ -711,7 +728,7 @@ void lossdistrib_joint_Z(double *dp, int *ndp, double *w,
     double beta = 0;
     int one = 1;
 
-#pragma omp parallel for private(j)
+    #pragma omp parallel for private(j)
     for(i = 0; i < *nZ; i++){
         for(j = 0; j < *ndp; j++){
             dpshocked[j + (*ndp) * i] = shockprob(dp[j], rho[j], Z[i], 0);
-- 
cgit v1.2.3-70-g09d2