From 6eac91bf7bcf3515e661d74b49e8179938bd1488 Mon Sep 17 00:00:00 2001
From: Guillaume Horel <guillaume.horel@serenitascapital.com>
Date: Wed, 7 May 2014 10:50:13 -0400
Subject: add GHquad function, get rid of statmod dependency

---
 R/tranche_functions.R | 32 ++++++++++++--------------------
 src/Makevars          |  2 +-
 src/lossdistrib.c     | 30 ++++++++++++++++++++++++++++++
 src/lossdistrib.h     |  2 ++
 4 files changed, 45 insertions(+), 21 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);
-- 
cgit v1.2.3-70-g09d2