try to compute shocks between two yieldcurves

2 files changed, 123 insertions, 0 deletions

diff --git a/pyisda/curve.pxd b/pyisda/curve.pxd
index 48790ff..ea01e08 100644
--- a/pyisda/curve.pxd
+++ b/pyisda/curve.pxd
@@ -237,3 +237,9 @@ cpdef enum Fmt:
     Plain = 0
     LZ4 = 1
     Packed = 2
+
+cdef void forward_rates(const TRatePt* arr, const TDate base_date, const int n, double* f, int *t) noexcept nogil
+
+cdef class CurveShock:
+    cdef vector[int] t
+    cdef vector[double] s
diff --git a/pyisda/curve.pyx b/pyisda/curve.pyx
index b0c0e29..3f2fda5 100644
--- a/pyisda/curve.pyx
+++ b/pyisda/curve.pyx
@@ -65,6 +65,16 @@ cdef double survival_prob(const TCurve* curve, TDate start_date, TDate maturity_
             u *= (1 + eps)
         return u if log else exp(u)
 
+cdef void forward_rates(const TRatePt* arr, const TDate base_date, const int n, double* f, int *t) noexcept nogil:
+    cdef:
+        size_t i
+
+    t[0] = arr[0].fDate - base_date
+    f[0] = arr[0].fRate
+    for i in range(1, n):
+        t[i] = arr[i].fDate - base_date
+        f[i] = (arr[i].fRate * t[i] - arr[i-1].fRate * t[i-1]) / (t[i] - t[i-1])
+
 cdef class Curve(object):
 
     cdef inline const TCurve* get_TCurve(self) noexcept nogil:
@@ -584,6 +594,113 @@ cdef class YieldCurve(Curve):
         yc.buf.reset(<char*>forward_curve, char_free)
         return yc
 
+    @cython.cdivision(True)
+    def __sub__(self, YieldCurve yc not None):
+        cdef:
+            const TCurve* curve1 = self.get_TCurve()
+            const TCurve* curve2 = yc.get_TCurve()
+            size_t i, j
+            double* f1
+            int* t1
+            double* f2
+            int* t2
+            int n1, n2
+            CurveShock cs = CurveShock()
+
+        n1 = curve1.fNumItems
+        n2 = curve2.fNumItems
+        t1 = <int*>malloc(n1 * sizeof(int))
+        f1 = <double*>malloc(n1 * sizeof(double))
+        t2 = <int*>malloc(n2 * sizeof(int))
+        f2 = <double*>malloc(n2 * sizeof(double))
+        if curve2.fBaseDate >= curve1.fBaseDate:
+            raise ValueError()
+        forward_rates(curve1.fArray, curve1.fBaseDate, n1, f1, t1)
+        forward_rates(curve2.fArray, curve2.fBaseDate, n2, f2, t2)
+
+        for i in range(n1):
+            t1[i] -= curve2.fBaseDate - curve1.fBaseDate
+
+        i = 0
+        j = 0
+        while i < n1 and j < n2:
+            cs.s.push_back(f1[j] / f2[i] - 1)
+            if t1[i] < t2[j]:
+                cs.t.push_back(t1[i])
+                i += 1
+            elif t2[j] < t1[i]:
+                cs.t.push_back(t2[j])
+                j += 1
+            else:
+                cs.t.push_back(t1[i])
+                i += 1
+                j += 1
+        free(t1)
+        free(f1)
+        free(t2)
+        free(f2)
+        return cs
+
+    @cython.cdivision(True)
+    def __add__(self, CurveShock shock not None):
+        cdef:
+            YieldCurve shocked_yc = YieldCurve.__new__(YieldCurve)
+            const TCurve* curve = self.get_TCurve()
+            TCurve* shocked_curve
+            size_t i, j
+            char* buf
+            int n1, n2
+            int t1, t2
+            double *f
+            vector[double] f_shocked
+            vector[int] t_shocked
+            double acc
+
+        n1 = curve.fNumItems
+        f = <double*>malloc(sizeof(double) * n1)
+        t = <int*>malloc(sizeof(int) * n1)
+        forward_rates(curve.fArray, curve.fBaseDate, n1, f, t)
+
+        i = 0
+        j = 0
+        while i < n1 and j < shock.t.size():
+            f_shocked.push_back(f[i] * (1 + shock.s[j]))
+            if t[i] < shock.t[j]:
+                i += 1
+                t_shocked.push_back(t[i])
+            elif shock.t[j] < t[i]:
+                j += 1
+                t_shocked.push_back(shock.t[j])
+            else:
+                i += 1
+                j += 1
+                t_shocked.push_back(t[i])
+        free(t)
+        free(f)
+
+        shocked_yc.buf_size = sizeof(TCurve) + f_shocked.size() * sizeof(TRatePt)
+        buf = <char*>malloc(self.buf_size)
+        shocked_yc.buf.reset(buf, char_free)
+        shocked_curve = <TCurve*>buf
+        shocked_curve.fNumItems = f_shocked.size()
+        shocked_curve.fDayCountConv = curve.fDayCountConv
+        shocked_curve.fBasis = curve.fBasis
+        shocked_curve.fBaseDate = curve.fBaseDate
+
+        t1 = 0
+        t2 = 0
+        acc = 0
+        for i in range(shocked_curve.fNumItems):
+            t2 = t_shocked[i]
+            acc += f_shocked[i] * (t2 - t1)
+            shocked_curve.fArray[i] = TRatePt(curve.fBaseDate + t2, acc / t2)
+            t1 = t2
+        return shocked_yc
+
+cdef class CurveShock:
+     def inspect(self):
+         return self.t, self.s
+
 @cython.cdivision(True)
 cdef void tweak_curve(const TCurve* sc, TCurve* sc_tweaked, double epsilon,
                       unsigned long mask) noexcept nogil: