add an function for computing truncated expected loss

2 files changed, 36 insertions, 1 deletions

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);