add GHquad function, get rid of statmod dependency

4 files changed, 44 insertions, 20 deletions

diff --git a/R/tranche_functions.R b/R/tranche_functions.R
index 2e7f0f3..5fa954a 100644
--- a/R/tranche_functions.R
+++ b/R/tranche_functions.R
@@ -1,5 +1,3 @@
-library(statmod)

-

 ## todo:

 ## -investigate other ways to interpolate the random severities on the grid

 ## I'm thinking that at eah severity that we add to the distribution, round it down

@@ -11,22 +9,16 @@ library(statmod)
 ## recovery with 0.1-1 support, so it's not clear that there is a big gain.

 ## - do the correlation adjustments when computing the deltas since it seems to be

 ## the market standard

-hostname <- system("hostname", intern=TRUE)

 

-checkSymbol <- function(name){

-    if(!is.loaded(name)){

-        if(.Platform$OS.type == "unix"){

-            root.dir <- "/home/share/CorpCDOs"

-            dyn.load(file.path(root.dir, "code", "R", paste0("lossdistrib",

-                                                             hostname,

-                                                             .Platform$dynlib.ext)))

-        }else{

-            root.dir <- "//WDSENTINEL/share/CorpCDOs"

-            dyn.load(file.path(root.dir, "code", "R", paste0("lossdistrib",

-                                                             .Platform$dynlib.ext)))

-        }

-    }

+GHquad <- function(n){

+    n <- as.integer(n)

+    Z <- double(n)

+    w <- double(n)

+    result <- .C("GHquad", n, Z=Z, w=w)

+    result[[1]] <- NULL

+    return(result)

 }

+

 lossdistrib <- function(p){

     ## basic recursive algorithm of Andersen, Sidenius and Basu

     n <- length(p)

@@ -630,9 +622,9 @@ tranche.pvvec <- function(K, L, R, cs){
 BClossdist <- function(defaultprob, issuerweights, recov, rho, Z, w,

                        N=length(recov)+1, defaultflag=FALSE, n.int=500){

   if(missing(Z)){

-    quadrature <- gauss.quad.prob(n.int, "normal")

-    Z <- quadrature$nodes

-    w <- quadrature$weights

+    quadrature <- GHquad(n.int)

+    Z <- quadrature$Z

+    w <- quadrature$w

   }

   ## do not use if weights are not gaussian, results would be incorrect

   ## since shockseverity is invalid in that case (need to use stochasticrecov)

diff --git a/src/Makevars b/src/Makevars
index f6790ff..0ef4741 100644
--- a/src/Makevars
+++ b/src/Makevars
@@ -1,2 +1,2 @@
 PKG_CFLAGS=$(SHLIB_OPENMP_CFLAGS)

-PKG_LIBS=$(SHLIB_OPENMP_CFLAGS)

+PKG_LIBS=$(SHLIB_OPENMP_CFLAGS) -llapack

diff --git a/src/lossdistrib.c b/src/lossdistrib.c
index 77ea335..f12ab2f 100644
--- a/src/lossdistrib.c
+++ b/src/lossdistrib.c
@@ -11,8 +11,38 @@ extern int dgemv_(char* trans, int *m, int *n, double* alpha, double* A, int* ld
 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

+    int INFO;

+    int i;

+    // Initialize array for workspace

+    double * WORK   = malloc(sizeof(double)*(2*(*n)-2));

+

+    // Initialize array for eigenvectors

+    double * V      = malloc(sizeof(double)*(*n)*(*n));

+

+    for(i = 0; i<(*n)-1; i++){

+      w[i] = sqrt((i+1.)/2);

+    }

+

+    // Run eigen decomposition

+    dstev_(&JOBZ, n, Z, w, V, n, WORK, &INFO);

+

+    for (i=0; i<(*n); i++) {

+      w[i] = V[i*(*n)] * V[i*(*n)];

+      Z[i] *= sqrt(2);

+    }

+

+    // Deallocate temporary arrays

+    free(WORK);

+    free(V);

+}

+

 void lossdistrib(double *p, int *np, double *w, double *S, int *N, int *defaultflag,

 		 double *q) {

   /* recursive algorithm with first order correction for computing

diff --git a/src/lossdistrib.h b/src/lossdistrib.h
index d4bd974..785b4e8 100644
--- a/src/lossdistrib.h
+++ b/src/lossdistrib.h
@@ -44,3 +44,5 @@ void BClossdist(double *SurvProb, int *dim1, int *dim2, double *issuerweights,
                 int *defaultflag, double *L, double *R);

 

 double quantile(double* Z, double* w, int nZ, double p0);

+

+void GHquad(int *n, double* Z, double* w);