Added Numba parallelisation to pricing via monte carlo

Author: domokane

MANIFEST,in renamed to MANIFEST.in

@@ -1,4 +1,4 @@
-## Quickstart Guide to Bonds
+## Quickstart Guide to Bond Analysis using financepy
 To analyse a bond first load up the Date and other utility classes such as Frequency and DayCounts. The easiest way to do this is to use a wildcard import. The downside is that it imports a lot of unnecessary classes and will take a second or so. But it is probably the simplest way to start until you become familiar with the library structure.
 
 ```

__init__.py

@@ -1 +1 @@
-__version__ = "1.0.0"
+__version__ = "1.0.1"

docs/QUICKSTART_BONDS.md

@@ -2,6 +2,8 @@
 
 from enum import Enum
 from numba import njit, float64, int64
+import numba as nb
+
 import numpy as np
 from ..utils.helpers import label_to_string
 
@@ -122,7 +124,12 @@ def draw(rt, a, b, sigma, dt):
 ########################################################################################
 
 
-@njit(float64[:](float64, float64, float64, float64, float64, float64, int64, int64))
+@njit(
+    float64[:](float64, float64, float64, float64, float64, float64, int64, int64),
+    parallel=True,
+    fastmath=True,
+    cache=True,
+)
 def rate_path_mc(r0, a, b, sigma, t, dt, seed, scheme):
     """Generate a path of CIR rates using a number of numerical schemes."""
 
@@ -136,7 +143,7 @@ def rate_path_mc(r0, a, b, sigma, t, dt, seed, scheme):
 
         sigmasqrt_dt = sigma * np.sqrt(dt)
 
-        for _ in range(0, num_paths):
+        for _ in nb.prange(num_paths):
 
             r = r0
             z = np.random.normal(0.0, 1.0, size=num_steps - 1)
@@ -154,7 +161,7 @@ def rate_path_mc(r0, a, b, sigma, t, dt, seed, scheme):
         x = np.exp(-a * dt)
         y = 1.0 - x
 
-        for _ in range(0, num_paths):
+        for _ in nb.prange(num_paths):
 
             r = r0
             z = np.random.normal(0.0, 1.0, size=num_steps - 1)
@@ -171,7 +178,7 @@ def rate_path_mc(r0, a, b, sigma, t, dt, seed, scheme):
         sigmasqrt_dt = sigma * np.sqrt(dt)
         sigma2dt = sigma * sigma * dt / 4.0
 
-        for _ in range(0, num_paths):
+        for _ in nb.prange(num_paths):
 
             r = r0
             z = np.random.normal(0.0, 1.0, size=num_steps - 1)
@@ -190,7 +197,7 @@ def rate_path_mc(r0, a, b, sigma, t, dt, seed, scheme):
         bhat = b - sigma * sigma / 4.0 / a
         sqrt_dt = np.sqrt(dt)
 
-        for _ in range(0, num_paths):
+        for _ in nb.prange(num_paths):
 
             r = r0
             z = np.random.normal(0.0, 1.0, size=num_steps - 1)
@@ -204,7 +211,7 @@ def rate_path_mc(r0, a, b, sigma, t, dt, seed, scheme):
 
     elif scheme == CIRNumericalScheme.EXACT.value:
 
-        for _ in range(0, num_paths):
+        for _ in nb.prange(num_paths):
 
             r = r0
 
@@ -229,7 +236,10 @@ def rate_path_mc(r0, a, b, sigma, t, dt, seed, scheme):
         int64,
         int64,
         int64,
-    )
+    ),
+    parallel=True,
+    fastmath=True,
+    cache=True,
 )
 def zero_price_mc(r0, a, b, sigma, t, dt, num_paths, seed, scheme):
     """Determine the CIR zero price using Monte Carlo."""
@@ -246,7 +256,7 @@ def zero_price_mc(r0, a, b, sigma, t, dt, num_paths, seed, scheme):
 
         sigmasqrt_dt = sigma * np.sqrt(dt)
 
-        for _ in range(0, num_paths):
+        for _ in nb.prange(num_paths):
 
             r = r0
             rsum = r
@@ -268,7 +278,7 @@ def zero_price_mc(r0, a, b, sigma, t, dt, num_paths, seed, scheme):
         x = np.exp(-a * dt)
         y = 1.0 - x
 
-        for _ in range(0, num_paths):
+        for _ in nb.prange(num_paths):
 
             r = r0
             rsum = r0
@@ -290,7 +300,7 @@ def zero_price_mc(r0, a, b, sigma, t, dt, num_paths, seed, scheme):
         sigmasqrt_dt = sigma * np.sqrt(dt)
         sigma2dt = sigma * sigma * dt / 4.0
 
-        for _ in range(0, num_paths):
+        for _ in nb.prange(num_paths):
 
             r = r0
             rsum = r
@@ -313,7 +323,7 @@ def zero_price_mc(r0, a, b, sigma, t, dt, num_paths, seed, scheme):
         bhat = b - sigma * sigma / 4.0 / a
         sqrt_dt = np.sqrt(dt)
 
-        for _ in range(0, num_paths):
+        for _ in nb.prange(num_paths):
 
             r = r0
             rsum = r
@@ -331,7 +341,7 @@ def zero_price_mc(r0, a, b, sigma, t, dt, num_paths, seed, scheme):
 
     elif scheme == CIRNumericalScheme.EXACT.value:
 
-        for _ in range(0, num_paths):
+        for _ in nb.prange(num_paths):
 
             r = r0
             rsum = r

financepy/__init__.py

@@ -2,6 +2,8 @@
 
 import numpy as np
 from numba import njit  # , float64, int64
+import numba as nb
+
 from ..utils.math import cholesky
 from ..utils.error import FinError
 
@@ -14,7 +16,7 @@
 ########################################################################################
 
 
-@njit(fastmath=True, cache=True)
+@njit(fastmath=True, cache=True, parallel=True)
 def get_paths_times(num_paths, num_time_steps, t, mu, stock_price, volatility, seed):
     """Get the simulated GBM process for a single asset with even num paths and
     time steps. Inputs include the number of time steps, paths, the drift mu,
@@ -46,7 +48,7 @@ def get_paths_times(num_paths, num_time_steps, t, mu, stock_price, volatility, s
 
     for it in range(1, num_time_steps + 1):
         g_1d = np.random.standard_normal((num_paths_even))
-        for ip in range(0, num_paths_even):
+        for ip in nb.prange(num_paths_even):
             w = np.exp(g_1d[ip] * vsqrt_dt)
             ip_start = ip * 2
             s_all[ip_start, it] = s_all[ip_start, it - 1] * m * w
@@ -121,7 +123,7 @@ def get_assets_paths_times(
     # or use g_corr = np.einsum('ijk,kl->ijl', g, c)
 
     # Calculate the Cholesky dot product
-    for ip in range(0, num_paths_even):
+    for ip in nb.prange(num_paths_even):
         for it in range(0, num_time_steps + 1):
             for ia in range(0, num_assets):
                 g_corr[ip][it][ia] = 0.0
@@ -131,7 +133,7 @@ def get_assets_paths_times(
     for ia in range(0, num_assets):
         s_all[ia, :, 0] = stock_prices[ia]
 
-    for ip in range(0, num_paths_even):
+    for ip in nb.prange(num_paths_even):
         ip_start = ip * 2
         for it in range(1, num_time_steps + 1):
             for ia in range(0, num_assets):
@@ -205,13 +207,13 @@ def get_assets_paths(
 
     # Calculate the dot product
     for ia in range(0, num_assets):
-        for ip in range(0, num_paths_even):
+        for ip in nb.prange(num_paths_even):
             g_corr[ip][ia] = 0.0
             for ib in range(0, num_assets):
                 g_corr[ip][ia] += g[ip][ib] * c[ia][ib]
 
     for ia in range(0, num_assets):
-        for ip in range(0, num_paths_even):
+        for ip in nb.prange(num_paths_even):
             z = g_corr[ip, ia]
             w = np.exp(z * vsqrt_dts[ia])
             ip_start = ip * 2

financepy/models/cir_montecarlo.py

@@ -56,6 +56,7 @@ class HestonNumericalScheme(Enum):
     ),
     cache=True,
     fastmath=True,
+    parallel=True,
 )
 def get_paths(s0, r, q, v0, kappa, theta, sigma, rho, t, dt, num_paths, seed, scheme):
 

financepy/models/gbm_process_simulator.py

@@ -4,6 +4,8 @@
 from enum import Enum
 
 from numba import njit, float64, int64
+import numba as nb
+
 import numpy as np
 
 from ..utils.error import FinError
@@ -137,6 +139,7 @@ class FinHestonNumericalScheme(Enum):
     ),
     cache=True,
     fastmath=True,
+    parallel=True,
 )
 def get_heston_paths(
     num_paths,
@@ -164,7 +167,7 @@ def get_heston_paths(
 
     if scheme == FinHestonNumericalScheme.EULER_SCHEME.value:
         # Basic scheme to first order with truncation on variance
-        for i_path in range(0, num_paths):
+        for i_path in nb.prange(num_paths):
             s = s0
             v = v0
             for i_step in range(1, num_steps + 1):
@@ -188,7 +191,7 @@ def get_heston_paths(
 
     elif scheme == FinHestonNumericalScheme.EULERLOG_SCHEME.value:
         # Basic scheme to first order with truncation on variance
-        for i_path in range(0, num_paths):
+        for i_path in nb.prange(num_paths):
             x = log(s0)
             v = v0
             for i_step in range(1, num_steps + 1):
@@ -220,7 +223,7 @@ def get_heston_paths(
         c1 = sigma2 * q * (1.0 - q) / kappa
         c2 = theta * sigma2 * ((1.0 - q) ** 2) / 2.0 / kappa
 
-        for i_path in range(0, num_paths):
+        for i_path in nb.prange(num_paths):
             x = log(s0)
             vn = v0
             for i_step in range(1, num_steps + 1):
@@ -283,6 +286,7 @@ class FinGBMNumericalScheme(Enum):
     float64[:, :](int64, int64, float64, float64, float64, float64, int64, int64),
     cache=True,
     fastmath=True,
+    parallel=True,
 )
 def get_gbm_paths(num_paths, num_annual_steps, t, mu, stock_price, sigma, scheme, seed):
 
@@ -298,7 +302,7 @@ def get_gbm_paths(num_paths, num_annual_steps, t, mu, stock_price, sigma, scheme
         s_all[:, 0] = stock_price
         for it in range(1, num_time_steps + 1):
             g1_d = np.random.standard_normal((num_paths))
-            for ip in range(0, num_paths):
+            for ip in nb.prange(num_paths):
                 w = np.exp(g1_d[ip] * vsqrt_dt)
                 s_all[ip, it] = s_all[ip, it - 1] * m * w
 
@@ -308,7 +312,7 @@ def get_gbm_paths(num_paths, num_annual_steps, t, mu, stock_price, sigma, scheme
         s_all[:, 0] = stock_price
         for it in range(1, num_time_steps + 1):
             g1_d = np.random.standard_normal((num_paths))
-            for ip in range(0, num_paths):
+            for ip in nb.prange(num_paths):
                 w = np.exp(g1_d[ip] * vsqrt_dt)
                 s_all[ip, it] = s_all[ip, it - 1] * m * w
                 s_all[ip + num_paths, it] = s_all[ip + num_paths, it - 1] * m / w
@@ -355,7 +359,7 @@ def get_vasicek_paths(
     if scheme == FinVasicekNumericalScheme.NORMAL.value:
         rate_path = np.empty((num_paths, num_steps + 1))
         rate_path[:, 0] = r0
-        for i_path in range(0, num_paths):
+        for i_path in nb.prange(num_paths):
             r = r0
             z = np.random.normal(0.0, 1.0, size=num_steps)
             for i_step in range(1, num_steps + 1):
@@ -364,7 +368,7 @@ def get_vasicek_paths(
     elif scheme == FinVasicekNumericalScheme.ANTITHETIC.value:
         rate_path = np.empty((2 * num_paths, num_steps + 1))
         rate_path[:, 0] = r0
-        for i_path in range(0, num_paths):
+        for i_path in nb.prange(num_paths):
             r1 = r0
             r2 = r0
             z = np.random.normal(0.0, 1.0, size=num_steps)
@@ -410,7 +414,7 @@ def get_cir_paths(
 
     if scheme == CIRNumericalScheme.EULER_SCHEME.value:
         sigma_sqrt_dt = sigma * sqrt(dt)
-        for i_path in range(0, num_paths):
+        for i_path in nb.prange(num_paths):
             r = r0
             z = np.random.normal(0.0, 1.0, size=num_steps)
             for i_step in range(1, num_steps + 1):
@@ -426,7 +430,7 @@ def get_cir_paths(
     elif scheme == CIRNumericalScheme.LOGNORMAL_SCHEME.value:
         x = exp(-kappa * dt)
         y = 1.0 - x
-        for i_path in range(0, num_paths):
+        for i_path in nb.prange(num_paths):
             r = r0
             z = np.random.normal(0.0, 1.0, size=num_steps)
             for i_step in range(1, num_steps + 1):
@@ -439,7 +443,7 @@ def get_cir_paths(
     elif scheme == CIRNumericalScheme.MILSTEIN_SCHEME.value:
         sigma_sqrt_dt = sigma * sqrt(dt)
         sigma2dt = sigma * sigma * dt / 4.0
-        for i_path in range(0, num_paths):
+        for i_path in nb.prange(num_paths):
             r = r0
             z = np.random.normal(0.0, 1.0, size=num_steps)
             for i_step in range(1, num_steps + 1):
@@ -455,7 +459,7 @@ def get_cir_paths(
     elif scheme == CIRNumericalScheme.KAHLJACKEL_SCHEME.value:
         bhat = theta - sigma * sigma / 4.0 / kappa
         sqrt_dt = sqrt(dt)
-        for i_path in range(0, num_paths):
+        for i_path in nb.prange(num_paths):
             r = r0
             z = np.random.normal(0.0, 1.0, size=num_steps)
             for i_step in range(1, num_steps + 1):

financepy/models/heston.py

@@ -2,6 +2,7 @@
 
 from math import sqrt, exp
 from numba import njit, float64, int64
+import numba as nb
 import numpy as np
 
 from ..utils.helpers import label_to_string
@@ -71,7 +72,12 @@ def zero_price(r0, a, b, sigma, t):
 ########################################################################################
 
 
-@njit(float64[:](float64, float64, float64, float64, float64, float64, int64))
+@njit(
+    float64[:](float64, float64, float64, float64, float64, float64, int64),
+    parallel=True,
+    fastmath=True,
+    cache=True,
+)
 def rate_path_mc(r0, a, b, sigma, t, dt, seed):
     """Generate a path of short rates using Vasicek model"""
 
@@ -83,7 +89,7 @@ def rate_path_mc(r0, a, b, sigma, t, dt, seed):
 
     sigmasqrt_dt = sigma * sqrt(dt)
 
-    for _ in range(0, num_paths):
+    for _ in nb.prange(num_paths):
 
         r = r0
         z = np.random.normal(0.0, 1.0, size=num_steps - 1)
@@ -100,16 +106,17 @@ def rate_path_mc(r0, a, b, sigma, t, dt, seed):
 
 @njit(
     float64(float64, float64, float64, float64, float64, float64, int64, int64),
-    fastmath=True,
     cache=True,
+    fastmath=True,
+    parallel=True,
 )
 def zero_price_mc(r0, a, b, sigma, t, dt, num_paths, seed):
     """Generate zero price by Monte Carlo using Vasicek model"""
     np.random.seed(seed)
     num_steps = int(t / dt)
     sigmasqrt_dt = sigma * sqrt(dt)
     zcb = 0.0
-    for _ in range(0, num_paths):
+    for _ in nb.prange(num_paths):
         z = np.random.normal(0.0, 1.0, size=num_steps)
         rsum = 0.0
         r = r0