diff options
Diffstat (limited to 'python')
| -rw-r--r-- | python/analytics/black.pxd | 2 | ||||
| -rw-r--r-- | python/analytics/black.pyx | 2 | ||||
| -rw-r--r-- | python/analytics/cms_spread.py | 96 | ||||
| -rw-r--r-- | python/analytics/cms_spread_utils.pyx | 68 | ||||
| -rw-r--r-- | python/tests/test_cms_spread.py | 16 |
5 files changed, 86 insertions, 98 deletions
diff --git a/python/analytics/black.pxd b/python/analytics/black.pxd new file mode 100644 index 00000000..e3b9fafa --- /dev/null +++ b/python/analytics/black.pxd @@ -0,0 +1,2 @@ +#cython: language_level=3 +cdef double cnd_erf(double d) nogil diff --git a/python/analytics/black.pyx b/python/analytics/black.pyx index 80d13a1a..9f733282 100644 --- a/python/analytics/black.pyx +++ b/python/analytics/black.pyx @@ -3,7 +3,7 @@ from libc.math cimport log, sqrt, erf from scipy.stats import norm import cython -cpdef double cnd_erf(double d): +cdef double cnd_erf(double d) nogil: """ 2 * Phi where Phi is the cdf of a Normal """ cdef double RSQRT2 = 0.7071067811865475 return 1 + erf(RSQRT2 * d) diff --git a/python/analytics/cms_spread.py b/python/analytics/cms_spread.py index 24c45e73..fa61c7cf 100644 --- a/python/analytics/cms_spread.py +++ b/python/analytics/cms_spread.py @@ -1,12 +1,12 @@ -import atexit +from . import cms_spread_utils +from .cms_spread_utils import h_call, h_put import datetime import matplotlib.pyplot as plt import numpy as np import pandas as pd import re from math import exp, sqrt, log, pi -from .black import bachelier, cnd_erf -from numba import cfunc, types, float64, vectorize +from .black import bachelier from quantlib.time.api import ( Date, Period, @@ -45,95 +45,7 @@ from utils.db import dawn_engine, serenitas_pool __all__ = ["CmsSpread"] -@vectorize( - [ - float64( - float64, - float64, - float64, - float64, - float64, - float64, - float64, - float64, - float64, - ) - ], - cache=True, - nopython=True, -) -def h_call(z, K, S1, S2, mu_x, mu_y, sigma_x, sigma_y, rho): - # conditionned on S2, integral wrt S1 - # z = (y - mu_y) / sigma_y - u1 = mu_x + rho * sigma_x * z - Ktilde = K + S2 * exp(mu_y + sigma_y * z) - u2 = log(Ktilde / S1) - - v = sigma_x * sqrt(1 - rho * rho) - v2 = sigma_x * sigma_x * (1 - rho * rho) - x = (u1 - u2) / v - return 0.5 * (S1 * exp(u1 + 0.5 * v2) * cnd_erf(x + v) - Ktilde * cnd_erf(x)) - - -@vectorize( - [ - float64( - float64, - float64, - float64, - float64, - float64, - float64, - float64, - float64, - float64, - ) - ], - cache=True, - nopython=True, -) -def h_put(z, K, S1, S2, mu_x, mu_y, sigma_x, sigma_y, rho): - # z = (y - mu_y) / sigma_y - u1 = mu_x + rho * sigma_x * z - Ktilde = K + S2 * exp(mu_y + sigma_y * z) - u2 = log(Ktilde / S1) - - v = sigma_x * sqrt(1 - rho * rho) - v2 = sigma_x * sigma_x * (1 - rho * rho) - x = (u2 - u1) / v - return 0.5 * (Ktilde * cnd_erf(x) - S1 * exp(u1 + 0.5 * v2) * cnd_erf(x - v)) - - -_sig = types.double(types.intc, types.CPointer(types.double)) - - -@cfunc(_sig, cache=True, nopython=True) -def _h1(n, args): - # z = (y - mu_y) / sigma_y - z = args[0] - K = args[1] - S1 = args[2] - S2 = args[3] - mu_x = args[4] - mu_y = args[5] - sigma_x = args[6] - sigma_y = args[7] - rho = args[8] - u1 = mu_x + rho * sigma_x * z - Ktilde = K + S2 * exp(mu_y + sigma_y * z) - u2 = log(Ktilde / S1) - - v = sigma_x * sqrt(1 - rho * rho) - v2 = sigma_x * sigma_x * (1 - rho * rho) - x = (u1 - u2) / v - return ( - 0.5 - * (S1 * exp(u1 + 0.5 * v2) * cnd_erf(x + v) - Ktilde * cnd_erf(x)) - * exp(-0.5 * z * z) - ) - - -_call_integrand = LowLevelCallable(_h1.ctypes) +_call_integrand = LowLevelCallable.from_cython(cms_spread_utils, "h1") def get_fixings(conn, tenor1, tenor2, fixing_date=None): diff --git a/python/analytics/cms_spread_utils.pyx b/python/analytics/cms_spread_utils.pyx new file mode 100644 index 00000000..17396054 --- /dev/null +++ b/python/analytics/cms_spread_utils.pyx @@ -0,0 +1,68 @@ +# cython: language_level=3, cdivision=True +from libc.math cimport log, exp, sqrt +from .black cimport cnd_erf +cimport cython +import numpy as np + +cdef api double h1(int n, double* data) nogil: + # z = (y - mu_y) / sigma_y + cdef: + double z = data[0] + double K = data[1] + double S1 = data[2] + double S2 = data[3] + double mu_x = data[4] + double mu_y = data[5] + double sigma_x = data[6] + double sigma_y = data[7] + double rho = data[8] + double u1, u2, Ktilde, v, v2, x + + u1 = mu_x + rho * sigma_x * z + Ktilde = K + S2 * exp(mu_y + sigma_y * z) + u2 = log(Ktilde / S1) + + v = sigma_x * sqrt(1 - rho * rho) + v2 = sigma_x * sigma_x * (1 - rho * rho) + x = (u1 - u2) / v + return ( + 0.5 + * (S1 * exp(u1 + 0.5 * v2) * cnd_erf(x + v) - Ktilde * cnd_erf(x)) + * exp(-0.5 * z * z) + ) + + +cpdef double[::1] h_call(double[::1] z, double K, double S1, double S2, double mu_x, double mu_y, double sigma_x, double sigma_y, double rho): + # conditionned on S2, integral wrt S1 + # z = (y - mu_y) / sigma_y + cdef double[::1] r = np.empty_like(z) + cdef: + double u1, u2, Ktilde, x + double v = sigma_x * sqrt(1 - rho * rho) + double v2 = sigma_x * sigma_x * (1 - rho * rho) + int i + + for i in range(z.shape[0]): + u1 = mu_x + rho * sigma_x * z[i] + Ktilde = K + S2 * exp(mu_y + sigma_y * z[i]) + u2 = log(Ktilde / S1) + x = (u1 - u2) / v + r[i] = 0.5 * (S1 * exp(u1 + 0.5 * v2) * cnd_erf(x + v) - Ktilde * cnd_erf(x)) + return r + +cpdef double[::1] h_put(double[::1] z, double K, double S1, double S2, double mu_x, double mu_y, double sigma_x, double sigma_y, double rho): + # z = (y - mu_y) / sigma_y + cdef: + double[::1] r = np.empty_like(z) + double u1, u2, Ktilde, x + double v = sigma_x * sqrt(1 - rho * rho) + double v2 = sigma_x * sigma_x * (1 - rho * rho) + int i + + for i in range(z.shape[0]): + u1 = mu_x + rho * sigma_x * z[i] + Ktilde = K + S2 * exp(mu_y + sigma_y * z[i]) + u2 = log(Ktilde / S1) + x = (u2 - u1) / v + r[i] = 0.5 * (Ktilde * cnd_erf(x) - S1 * exp(u1 + 0.5 * v2) * cnd_erf(x - v)) + return r diff --git a/python/tests/test_cms_spread.py b/python/tests/test_cms_spread.py index 9ee89bbb..dbf3f218 100644 --- a/python/tests/test_cms_spread.py +++ b/python/tests/test_cms_spread.py @@ -16,7 +16,6 @@ from quantlib.experimental.coupons.cms_spread_coupon import ( CappedFlooredCmsSpreadCoupon) from quantlib.experimental.coupons.lognormal_cmsspread_pricer import ( LognormalCmsSpreadPricer) -from scipy import LowLevelCallable from scipy.special import roots_hermitenorm from scipy.integrate import quad import ctypes @@ -90,13 +89,20 @@ class TestCmsSpread(unittest.TestCase): def test_h1_hcall(self): args = (self.cap, *self.params, self.ρ.value) - h1_fun = _call_integrand.function - + h1_capsule = _call_integrand.function + ctypes.pythonapi.PyCapsule_GetPointer.restype = ctypes.c_void_p + ctypes.pythonapi.PyCapsule_GetPointer.argtypes = [ctypes.py_object, ctypes.c_char_p] + ptr_type = ctypes.CFUNCTYPE(ctypes.c_double, ctypes.c_int, ctypes.POINTER(ctypes.c_double)) + # Get the function pointer from the capsule "capsule_mul" + h1_fun = ptr_type(ctypes.pythonapi.PyCapsule_GetPointer(h1_capsule, b"double (int, double *)")) + x = np.linspace(-5, 5, 11) + b = h_call(x, *args) * np.exp(-0.5* x * x) + i = 0 for x in np.linspace(-5, 5, 11): full_args = np.array((x, *args)) a = h1_fun(9, full_args.ctypes.data_as(ctypes.POINTER(ctypes.c_double))) - b = h_call(x, *args) * math.exp(-0.5 * x * x) - self.assertAlmostEqual(a, b) + self.assertAlmostEqual(a, b[i]) + i+= 1 def test_scipy_integrate(self): x, w = roots_hermitenorm(20) |