add fater version

4 files changed, 200 insertions, 6 deletions

diff --git a/c_layer/cdsbootstrap.c b/c_layer/cdsbootstrap.c
new file mode 100644
index 0000000..f3ae975
--- /dev/null
+++ b/c_layer/cdsbootstrap.c
@@ -0,0 +1,50 @@
+#include "cdsbootstrap.h"
+
+int cdsBootstrapPointFunction
+(double   hazardRate,
+ void    *data,
+ double  *pv)
+{
+    int         status    = FAILURE;
+
+    cds_bootstrap_ctx *context = (cds_bootstrap_ctx*)data;
+
+    TCurve         *discountCurve = context->discountCurve;
+    TCurve         *cdsCurve      = context->cdsCurve;
+    TDate           cdsBaseDate   = cdsCurve->fBaseDate;
+    TDate           stepinDate    = context->stepinDate;
+    TDate           cashSettleDate = context->cashSettleDate;
+
+    double          pvC; /* PV of contingent leg */
+    double          pvF; /* PV of fee leg */
+    cdsCurve->fArray[0].fRate = hazardRate;
+
+    if (JpmcdsContingentLegPV (context->cl,
+                               cdsBaseDate,
+                               cashSettleDate,
+                               stepinDate,
+                               discountCurve,
+                               cdsCurve,
+                               context->recoveryRate,
+                               &pvC) != SUCCESS)
+        goto done;
+
+    if (JpmcdsFeeLegPV(context->fl,
+                       cdsBaseDate,
+                       stepinDate,
+                       cashSettleDate,
+                       discountCurve,
+                       cdsCurve,
+                       1,
+                       &pvF) != SUCCESS)
+        goto done;
+    /* Note: price is discounted to cdsBaseDate */
+    *pv = pvC - context->spread * pvF;
+    status = SUCCESS;
+
+ done:
+
+    if (status != SUCCESS)
+        printf("Something went wrong");
+    return status;
+}
diff --git a/c_layer/cdsbootstrap.h b/c_layer/cdsbootstrap.h
new file mode 100644
index 0000000..2e19189
--- /dev/null
+++ b/c_layer/cdsbootstrap.h
@@ -0,0 +1,26 @@
+#include "isda/ldate.h"
+#include "isda/cxzerocurve.h"
+#include "isda/feeleg.h"
+#include "isda/contingentleg.h"
+#include "stdio.h"
+
+
+#define SUCCESS 0
+#define FAILURE -1
+
+typedef struct
+{
+    TDate           stepinDate;
+    TDate           cashSettleDate;
+    TCurve         *discountCurve;
+    TCurve         *cdsCurve;
+    double          recoveryRate;
+    double          spread;
+    TContingentLeg *cl;
+    TFeeLeg        *fl;
+} cds_bootstrap_ctx;
+
+
+int cdsBootstrapPointFunction(double   hazardRate,
+                              void    *data,
+                              double  *pv);
diff --git a/pyisda/flat_hazard.pyx b/pyisda/flat_hazard.pyx
index 7549a38..2257098 100644
--- a/pyisda/flat_hazard.pyx
+++ b/pyisda/flat_hazard.pyx
@@ -8,8 +8,8 @@ from curve cimport (TCurve, YieldCurve, SpreadCurve, JpmcdsCleanSpreadCurve,
                     JpmcdsFreeTCurve, JpmcdsMakeTCurve, CONTINUOUS)
 from date cimport pydate_to_TDate, ACT_360, ACT_365F
 from cdsone cimport JpmcdsStringToStubMethod, TStubMethod
-from cpython.mem cimport PyMem_Malloc, PyMem_Free
-from libc.stdlib cimport free
+from libc.stdlib cimport free, malloc
+from libc.stdio cimport printf
 cimport cython
 import numpy as np
 
@@ -78,3 +78,113 @@ def strike_vec(double[:] spreads, YieldCurve yc, trade_date, value_date,
     if sc is not NULL:
         JpmcdsFreeTCurve(sc)
     return default_leg_pv, coupon_leg_pv
+
+cdef extern from "cdsbootstrap.h":
+    ctypedef struct cds_bootstrap_ctx:
+        TDate           stepinDate
+        TDate           cashSettleDate
+        TCurve         *discountCurve
+        TCurve         *cdsCurve;
+        double          recoveryRate
+        double          spread
+        TContingentLeg *cl
+        TFeeLeg        *fl
+    int cdsBootstrapPointFunction(...) nogil
+
+cdef extern from "isda/rtbrent.h" nogil:
+    ctypedef int (*TObjectFunc) (double x, void * para, double *f)
+
+    int JpmcdsRootFindBrent(
+        TObjectFunc funcd,              # (I) function to be solved */
+        void        *data,              # (I) data to pass into funcd */
+        double      boundLo,            # (I) lower bound on legal X */
+        double      boundHi,            # (I) upper bound on legal X */
+        int         numIterations,      # (I) Maximum number of iterations */
+        double      guess,              # (I) Initial guess */
+        double      initalXStep,        # (I) Size of step in x */
+        double      initialFDeriv,      # (I) Derivative, defaults to 0 */
+        double      xacc,               # (I) X accuracy tolerance */
+        double      facc,               # (I) function accuracy tolerance */
+        double      *solution)          # (O) root found */
+
+@cython.boundscheck(False)
+@cython.cdivision(True)
+def faster_strike_vec(double[:] spreads, YieldCurve yc, trade_date, value_date,
+                      start_date, end_date, double recovery_rate):
+    """Computes coupon and default leg for a vector of spreads"""
+
+    cdef TDate trade_date_c = pydate_to_TDate(trade_date)
+    cdef TDate step_in_date_c = trade_date_c + 1
+    cdef TDate value_date_c = pydate_to_TDate(value_date)
+    cdef TDate start_date_c = pydate_to_TDate(start_date)
+    cdef TDate end_date_c = pydate_to_TDate(end_date)
+    cdef np.ndarray[np.float_t] coupon_leg_pv = np.empty_like(spreads)
+    cdef np.ndarray[np.float_t] default_leg_pv = np.empty_like(spreads)
+    cdef int nspreads = spreads.shape[0]
+    cdef double h, guess
+    cdef TStubMethod stub_type
+    cdef TCurve* sc
+    cdef TContingentLeg* default_leg
+    cdef TFeeLeg* fee_leg
+    cdef cds_bootstrap_ctx* params
+    cdef size_t i
+    JpmcdsStringToStubMethod(b"f/s", &stub_type)
+
+    with nogil:
+        h = spreads[0] / (1 - recovery_rate)
+        sc = JpmcdsMakeTCurve(trade_date_c,
+                              &end_date_c,
+                              &h,
+                              1,
+                              <double>CONTINUOUS,
+                              ACT_365F) # JPMCDS_ACT_365F
+
+        default_leg = JpmcdsCdsContingentLegMake(start_date_c,
+                                                 end_date_c,
+                                                 1., # notional
+                                                 1) # protect_start = True
+        fee_leg = JpmcdsCdsFeeLegMake(start_date_c,
+                                      end_date_c,
+                                      1, # payAccOnDefault = True
+                                      NULL,
+                                      &stub_type,
+                                      1., # notional
+                                      1., # coupon_rate
+                                      ACT_360,  # JPMCDS_ACT_360
+                                      <long>'M',
+                                      b'NONE',
+                                      1) # protect_start = True
+        params = <cds_bootstrap_ctx*>malloc(sizeof(cds_bootstrap_ctx))
+        params.stepinDate = step_in_date_c
+        params.cashSettleDate = value_date_c
+        params.discountCurve = yc._thisptr
+        params.cdsCurve = sc
+        params.recoveryRate = recovery_rate
+        params.cl = default_leg
+        params.fl = fee_leg
+
+        for i in range(nspreads):
+            params.spread = spreads[i]
+            guess = spreads[i] / ( 1 - recovery_rate) * 365/360
+            if JpmcdsRootFindBrent(<TObjectFunc>cdsBootstrapPointFunction,
+                                   <void*> params,
+                                   0.0,   # boundLo */
+                                   1e10,   # boundHi */
+                                   100,    # numIterations */
+                                   guess,
+                                   0.0005, # initialXstep */
+                                   0,      # initialFDeriv */
+                                   1e-10,  # xacc */
+                                   1e-10,  # facc */
+                                   &h) != 0:
+                printf("failed to find the root")
+            sc.fArray[0].fRate = h
+            if JpmcdsFeeLegPV(fee_leg, trade_date_c, step_in_date_c, value_date_c,
+                              yc._thisptr, sc, True, &coupon_leg_pv[i]) != 0:
+                printf("Something went wrong")
+            default_leg_pv[i] = coupon_leg_pv[i] * spreads[i]
+        free(params)
+        JpmcdsFeeLegFree(fee_leg)
+        free(default_leg)
+        JpmcdsFreeTCurve(sc)
+    return default_leg_pv, coupon_leg_pv
diff --git a/setup.py b/setup.py
index 7c0b4f6..5958dd8 100644
--- a/setup.py
+++ b/setup.py
@@ -2,12 +2,20 @@ from distutils.core import setup
 from distutils.extension import Extension
 from Cython.Build import cythonize
 
-extensions = Extension("*", ["pyisda/*.pyx"],
-        libraries = ["cds"])
+all_extensions = Extension("*", ["pyisda/*.pyx"],
+                           include_dirs = ['c_layer'],
+                           libraries = ["cds"])
+
+c_extension = Extension("pyisda.flat_hazard",
+                        include_dirs = ['c_layer'],
+                        sources = ['pyisda/flat_hazard.pyx', 'c_layer/cdsbootstrap.c'],
+                        libraries = ['cds'])
+all_extensions = cythonize([c_extension, all_extensions], nthreads = 4,
+                           compiler_directives={'embedsignature':True})
+
 setup(
     name = "pyisda",
     version = '0.1',
     author = 'Guillaume Horel',
-    ext_modules = cythonize(extensions, nthreads = 4,
-        compiler_directives={'embedsignature':True}),
+    ext_modules = all_extensions,
     packages = ['pyisda'])