add dispersion method

1 files changed, 82 insertions, 1 deletions

diff --git a/pyisda/credit_index.pyx b/pyisda/credit_index.pyx
index e2bd58b..dbb9a6d 100644
--- a/pyisda/credit_index.pyx
+++ b/pyisda/credit_index.pyx
@@ -1,6 +1,6 @@
 #cython: cdivision=True, boundscheck=False, c_string_type=unicode, c_string_encoding=ascii
 from libc.stdlib cimport malloc, free
-from libc.math cimport nan, isnan
+from libc.math cimport nan, isnan, log
 from libc.string cimport memcpy, memset
 from libcpp.pair cimport pair
 from libcpp.memory cimport make_shared
@@ -667,6 +667,87 @@ cdef class CreditIndex(CurveList):
                     &spreads_view[p.second, 0])
         return spreads
 
+    def dispersion(self, YieldCurve yc not None, const TDate[::1] schedule=None, const bint use_log=True):
+        cdef:
+            shared_ptr[TCurve] sc
+            pair[CurveName, size_t] p
+            np.npy_intp n
+            const TDate* schedule_ptr
+            size_t i
+
+        if schedule is None:
+            schedule_ptr = self._maturities.const_data()
+            n = self._maturities.size()
+        else:
+            schedule_ptr = &schedule[0]
+            n = schedule.shape[0]
+
+        cdef:
+            np.ndarray std_spreads = np.PyArray_ZEROS(1, &n, np.NPY_DOUBLE, 0)
+            double[::1] std_spreads_view = std_spreads
+            TStubMethod stub_type
+            double Mk, w
+            double* spreads = <double*>malloc(n * sizeof(double))
+            double* M = <double*>malloc(n * sizeof(double))
+            double w_sum, w_sum2
+            bint first_time = True
+
+        JpmcdsStringToStubMethod(b"f/s", &stub_type)
+
+        with nogil:
+
+            for p in self.names:
+                if self.defaulted[p.second] != -1:
+                    continue
+                sc = self._curves[p.second]
+                w = self._weights[p.second]
+                JpmcdsCdsParSpreads(
+                    self.base_date,
+                    self.base_date + 1,
+                    self._start_date,
+                    n,
+                    schedule_ptr,
+                    True,
+                    NULL,
+                    &stub_type,
+                    ACT_360,
+                    MODIFIED,
+                    b'NONE',
+                    get_TCurve(yc),
+                    sc.get(),
+                    self.recovery_rates[p.second].get(),
+                    spreads)
+
+                if first_time:
+                    if use_log:
+                        for i in range(n):
+                            M[i] = log(spreads[i])
+                    else:
+                        memcpy(M, spreads, n * sizeof(double))
+                    w_sum = w
+                    w_sum2 = w * w
+                    first_time = False
+                    continue
+                w_sum += w
+                w_sum2 += w * w
+                if use_log:
+                    for i in range(n):
+                        S = log(spreads[i])
+                        Mk = M[i] + (w / w_sum) * (S - M[i])
+                        std_spreads_view[i] += w * (S - M[i]) * (S - Mk)
+                        M[i] = Mk
+                else:
+                    for i in range(n):
+                        S = spreads[i]
+                        Mk = M[i] + (w / w_sum) * (S - M[i])
+                        std_spreads_view[i] += w * (S - M[i]) * (S - Mk)
+                        M[i] = Mk
+            for i in range(n):
+                std_spreads_view[i] /= (w_sum - w_sum2 / w_sum)
+            free(M)
+            free(spreads)
+        return std_spreads
+
 cdef unsigned long fill_mask(const TDate maturity, const vector[TDate]& maturities,
                              const shared_ptr[TCurve]& sc) nogil:
     cdef: