rework flat_hazard into an optim module

3 files changed, 196 insertions, 197 deletions

diff --git a/pyisda/flat_hazard.pyx b/pyisda/flat_hazard.pyx
deleted file mode 100644
index 88cd82a..0000000
--- a/pyisda/flat_hazard.pyx
+++ /dev/null
@@ -1,195 +0,0 @@
-from date cimport dcc, TDate
-import numpy as np
-cimport numpy as np
-
-from legs cimport (JpmcdsCdsContingentLegMake, JpmcdsCdsFeeLegMake,
-                   TContingentLeg, TFeeLeg, JpmcdsContingentLegPV, JpmcdsFeeLegPV,
-                   JpmcdsFeeLegFree)
-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 libc.stdlib cimport free, malloc
-from libc.stdio cimport printf
-cimport cython
-
-np.import_array()
-
-@cython.boundscheck(False)
-def 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 TCurve* sc
-    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 TStubMethod stub_type
-
-    JpmcdsStringToStubMethod(b"f/s", &stub_type)
-
-    cdef TContingentLeg* default_leg = JpmcdsCdsContingentLegMake(
-        start_date_c,
-        end_date_c,
-        1, # notional
-        1) # protect_start = True
-    cdef TFeeLeg* fee_leg = JpmcdsCdsFeeLegMake(start_date_c,
-                                                end_date_c,
-                                                1,
-                                                NULL,
-                                                &stub_type,
-                                                1., # notional
-                                                1., # coupon_rate
-                                                ACT_360,  # JPMCDS_ACT_360
-                                                <long>'M',
-                                                b'NONE',
-                                                1) # protect_start = True
-    cdef size_t i
-    cdef double zero_upfront = 0
-    for i in range(spreads.shape[0]):
-        sc = JpmcdsCleanSpreadCurve(trade_date_c,
-                                    yc._thisptr.get(),
-                                    start_date_c,
-                                    step_in_date_c,
-                                    value_date_c,
-                                    1, # length of end_dates_c
-                                    &end_date_c,
-                                    &spreads[i],
-                                    &zero_upfront,
-                                    0,
-                                    &recovery_rate,
-                                    1, # pay_accrued_on_default True
-                                    NULL,
-                                    ACT_360, # JPMCDS_ACT_360
-                                    &stub_type,
-                                    <long>'M',
-                                    b'NONE')
-        if JpmcdsContingentLegPV(default_leg, trade_date_c, value_date_c,
-                                 step_in_date_c, yc._thisptr.get(), sc,
-                                 recovery_rate, &default_leg_pv[i]) != 0:
-            raise ValueError("Something went wrong")
-        if JpmcdsFeeLegPV(fee_leg, trade_date_c, step_in_date_c, value_date_c,
-                          yc._thisptr.get(), sc, True, &coupon_leg_pv[i]) != 0:
-            raise ValueError("Something went wrong")
-    JpmcdsFeeLegFree(fee_leg)
-    free(default_leg)
-    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 pv_vec(double[:] spreads, YieldCurve yc, trade_date, value_date,
-           start_date, end_date, double recovery_rate, double fixed_rate):
-    """Computes vector of cds clean pvs 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.npy_intp nspreads = spreads.shape[0]
-    cdef double h, guess, hi, lo, coupon_leg_pv
-    cdef np.ndarray[np.float_t] pv = np.empty_like(spreads)
-    cdef TStubMethod stub_type
-    cdef TCurve* sc
-    cdef TContingentLeg* default_leg
-    cdef TFeeLeg* fee_leg
-    cdef cds_bootstrap_ctx* params
-    cdef Py_ssize_t i
-    JpmcdsStringToStubMethod(b"f/s", &stub_type)
-    cdef double eta = 365./360.
-
-    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.get()
-        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) * eta
-            lo = guess * 0.9
-            hi = guess * 1.1
-            if JpmcdsRootFindBrent(<TObjectFunc>cdsBootstrapPointFunction,
-                                   <void*> params,
-                                   lo,   # boundLo */
-                                   hi,   # boundHi */
-                                   100,    # numIterations */
-                                   guess,
-                                   0.0005, # initialXstep */
-                                   0,      # initialFDeriv */
-                                   1e-10,  # xacc */
-                                   1e-10,  # facc */
-                                   &h) != 0:
-                printf("Failed to find hazard rate for: %f\n", spreads[i])
-                break
-            sc.fArray[0].fRate = h
-            if JpmcdsFeeLegPV(fee_leg, trade_date_c, step_in_date_c, value_date_c,
-                              yc._thisptr.get(), sc, True, &coupon_leg_pv) != 0:
-                printf("Something went wrong\n")
-            pv[i] = coupon_leg_pv * (spreads[i] - fixed_rate)
-        free(params)
-        JpmcdsFeeLegFree(fee_leg)
-        free(default_leg)
-        JpmcdsFreeTCurve(sc)
-    return pv
diff --git a/pyisda/optim.pyx b/pyisda/optim.pyx
new file mode 100644
index 0000000..540a10b
--- /dev/null
+++ b/pyisda/optim.pyx
@@ -0,0 +1,194 @@
+from libc.stdlib cimport abort
+from libc.stdlib cimport free, malloc
+from libc.stdio cimport printf
+from libc.math cimport exp, sqrt
+from cdsone cimport JpmcdsStringToStubMethod, TStubMethod
+from curve cimport (TCurve, YieldCurve, JpmcdsFreeTCurve, JpmcdsNewTCurve,
+                    CONTINUOUS)
+from date cimport TDate, pydate_to_TDate, ACT_360, ACT_365F
+from legs cimport (JpmcdsCdsContingentLegMake, JpmcdsCdsFeeLegMake,
+                   TContingentLeg, JpmcdsFeeLegPV, JpmcdsFeeLegFree, TFeeLeg)
+
+cimport cython
+from cpython.pycapsule cimport PyCapsule_New, PyCapsule_GetPointer, PyCapsule_Destructor
+
+cdef extern from "cdsbootstrap.h" nogil:
+    ctypedef struct cds_bootstrap_ctx:
+        TDate           stepinDate
+        TDate           cashSettleDate
+        TCurve         *discountCurve
+        TCurve         *cdsCurve;
+        double          recoveryRate
+        double          spread
+        TContingentLeg *cl
+        TFeeLeg        *fl
+    int cdsBootstrapPointFunction(...)
+
+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 */
+
+ctypedef struct my_ctx:
+    cds_bootstrap_ctx params
+    double fixed_rate
+    double G
+    double sigmaT
+    double S0
+
+cdef free_my_ctx(object cap):
+    cdef my_ctx* params_ext = <my_ctx*> \
+        PyCapsule_GetPointer(cap, "cds_bootstrap_context")
+    JpmcdsFeeLegFree(params_ext.params.fl)
+    JpmcdsFreeTCurve(params_ext.params.cdsCurve)
+    free(params_ext)
+
+
+def init_context(YieldCurve yc not None, trade_date, value_date, start_date,
+                 end_date, double recovery_rate, double fixed_rate, double G,
+                 double sigmaT, double S0):
+    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 TStubMethod stub_type
+    cdef TContingentLeg* default_leg
+    cdef TFeeLeg* fee_leg
+
+    cdef my_ctx* params_ext = <my_ctx*>malloc(sizeof(my_ctx))
+    cdef cds_bootstrap_ctx* params = &params_ext.params
+    params_ext.fixed_rate = fixed_rate
+    params_ext.G = G
+    params_ext.sigmaT = sigmaT
+    params_ext.S0 = S0
+
+    JpmcdsStringToStubMethod(b"f/s", &stub_type)
+    cdef double eta = 365./360.
+    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.stepinDate = step_in_date_c
+    params.cashSettleDate = value_date_c
+    params.discountCurve = yc._thisptr.get()
+    params.cdsCurve = JpmcdsNewTCurve(trade_date_c,
+                                      1,
+                                      <double>CONTINUOUS,
+                                      ACT_365F) # JPMCDS_ACT_365F
+    params.cdsCurve.fArray[0].fDate = end_date_c
+    params.recoveryRate = recovery_rate
+    params.cl = default_leg
+    params.fl = fee_leg
+    return PyCapsule_New(<void*>params_ext, "cds_bootstrap_context",
+                         <PyCapsule_Destructor>free_my_ctx)
+
+
+@cython.cdivision(True)
+cdef api double pv(double Z, void* ctx):
+    cdef:
+        my_ctx* params_ext = <my_ctx*>ctx
+        cds_bootstrap_ctx* params = &params_ext.params
+        double sigmaT = params_ext.sigmaT
+        double S = params_ext.S0 * exp(-0.5 * sigmaT * sigmaT + sigmaT * Z)
+        double guess = S / ( 1 - params.recoveryRate) * 365. / 360.
+        double lo = guess * 0.9
+        double hi = guess * 1.1
+        double h
+        double coupon_leg_pv
+    params.spread = S
+    if JpmcdsRootFindBrent(<TObjectFunc>cdsBootstrapPointFunction,
+                           <void*>params,
+                           lo,   # boundLo */
+                           hi,   # boundHi */
+                           100,    # numIterations */
+                           guess,
+                           0.0005, # initialXstep */
+                           0,      # initialFDeriv */
+                           1e-10,  # xacc */
+                           1e-10,  # facc */
+                           &h) != 0:
+        printf("Failed to find hazard rate for: %f\n", S)
+        abort()
+    params.cdsCurve.fArray[0].fRate = h
+    if JpmcdsFeeLegPV(params.fl, params.cdsCurve.fBaseDate, params.stepinDate,
+                      params.cashSettleDate,
+                      params.discountCurve, params.cdsCurve, True, &coupon_leg_pv) != 0:
+        printf("Something went wrong\n")
+        abort()
+    return 0.3989422804014327 * exp(-0.5 * Z * Z) * \
+        ((S - params_ext.fixed_rate) * coupon_leg_pv - params_ext.G)
+
+def update_context(object ctx, double S0):
+    cdef my_ctx* params_ext = <my_ctx*>PyCapsule_GetPointer(ctx, "cds_bootstrap_context")
+    params_ext.S0 = S0 * 1e-4
+
+@cython.boundscheck(False)
+@cython.cdivision(True)
+def expected_pv(double[:] tilt, double[:] w, double S0, object ctx):
+    """ computes :math:`\frac{1}{2\pi}\int_{-\inf}^{\inf} PV(S) \exp(-\frac{z^2}{2}) dz`
+    where :math:`S = S_0\exp(z\sigma\sqrt{T} - \frac{\sigma^2 T}{2}`
+    using Gauss-Hermite quadrature. PV(S) is the PV of a swap with spread S.
+
+    tilt : for efficiency reason, is a vector containing :math:`\exp(z\sigma\sqrt{T} - \frac{\sigma^2 T}{2}`
+    for the Gauss-Hermite abcisses.
+    w : Gauss-Hermite weights.
+    """
+
+    cdef:
+        my_ctx* params_ext = <my_ctx*>PyCapsule_GetPointer(ctx, "cds_bootstrap_context")
+        cds_bootstrap_ctx* params = &params_ext.params
+        double sigmaT = params_ext.sigmaT
+        double S, guess, lo, hi, h, coupon_leg_pv, z
+        size_t i
+        double r = 0.
+
+    for i in range(tilt.shape[0]):
+        S = S0 * tilt[i] * 1e-4
+        guess = S / ( 1 - params.recoveryRate) * 365. / 360.
+        lo = guess * 0.9
+        hi = guess * 1.1
+        params.spread = S
+        if JpmcdsRootFindBrent(<TObjectFunc>cdsBootstrapPointFunction,
+                               <void*>params,
+                               lo,   # boundLo */
+                               hi,   # boundHi */
+                               100,    # numIterations */
+                               guess,
+                               0.0005, # initialXstep */
+                               0,      # initialFDeriv */
+                               1e-10,  # xacc */
+                               1e-10,  # facc */
+                               &h) != 0:
+            printf("Failed to find hazard rate for: %f\n", S)
+            abort()
+        params.cdsCurve.fArray[0].fRate = h
+        if JpmcdsFeeLegPV(params.fl, params.cdsCurve.fBaseDate, params.stepinDate,
+                          params.cashSettleDate,
+                          params.discountCurve, params.cdsCurve, True, &coupon_leg_pv) != 0:
+            printf("Something went wrong\n")
+            abort()
+        r += w[i] * ((S - params_ext.fixed_rate) * coupon_leg_pv)
+    return r
diff --git a/setup.py b/setup.py
index 359540d..ff8a420 100644
--- a/setup.py
+++ b/setup.py
@@ -10,9 +10,9 @@ all_extensions = Extension("*", ["pyisda/*.pyx"],
                            extra_compile_args=['-fopenmp'],
                            extra_link_args=['-fopenmp', '-Wl,--strip-all'])
 
-c_extension = Extension("pyisda.flat_hazard",
+c_extension = Extension("pyisda.optim",
                         include_dirs=['c_layer', numpy.get_include()],
-                        sources=['pyisda/flat_hazard.pyx', 'c_layer/cdsbootstrap.c'],
+                        sources=['pyisda/optim.pyx', 'c_layer/cdsbootstrap.c'],
                         libraries=['cds'],
                         language='c++')