1import math 2 3import numpy as np 4 5import unittest 6from numba.core.compiler import compile_isolated, compile_extra, Flags 7from numba.core import types, typing 8from numba.tests.support import TestCase 9 10 11RISKFREE = 0.02 12VOLATILITY = 0.30 13 14 15A1 = 0.31938153 16A2 = -0.356563782 17A3 = 1.781477937 18A4 = -1.821255978 19A5 = 1.330274429 20RSQRT2PI = 0.39894228040143267793994605993438 21 22 23def cnd_array(d): 24 K = 1.0 / (1.0 + 0.2316419 * np.abs(d)) 25 ret_val = (RSQRT2PI * np.exp(-0.5 * d * d) * 26 (K * (A1 + K * (A2 + K * (A3 + K * (A4 + K * A5)))))) 27 return np.where(d > 0, 1.0 - ret_val, ret_val) 28 29 30def cnd(d): 31 K = 1.0 / (1.0 + 0.2316419 * math.fabs(d)) 32 ret_val = (RSQRT2PI * math.exp(-0.5 * d * d) * 33 (K * (A1 + K * (A2 + K * (A3 + K * (A4 + K * A5)))))) 34 if d > 0: 35 ret_val = 1.0 - ret_val 36 return ret_val 37 38 39def blackscholes_arrayexpr(stockPrice, optionStrike, optionYears, Riskfree, 40 Volatility): 41 S = stockPrice 42 X = optionStrike 43 T = optionYears 44 R = Riskfree 45 V = Volatility 46 sqrtT = np.sqrt(T) 47 d1 = (np.log(S / X) + (R + 0.5 * V * V) * T) / (V * sqrtT) 48 d2 = d1 - V * sqrtT 49 cndd1 = cnd_array(d1) 50 cndd2 = cnd_array(d2) 51 52 expRT = np.exp(- R * T) 53 54 callResult = (S * cndd1 - X * expRT * cndd2) 55 putResult = (X * expRT * (1.0 - cndd2) - S * (1.0 - cndd1)) 56 return callResult, putResult 57 58 59def blackscholes_arrayexpr_jitted(stockPrice, optionStrike, optionYears, 60 Riskfree, Volatility): 61 S = stockPrice 62 X = optionStrike 63 T = optionYears 64 R = Riskfree 65 V = Volatility 66 sqrtT = np.sqrt(T) 67 d1 = (np.log(S / X) + (R + 0.5 * V * V) * T) / (V * sqrtT) 68 d2 = d1 - V * sqrtT 69 cndd1 = cnd_array_jitted(d1) 70 cndd2 = cnd_array_jitted(d2) 71 72 expRT = np.exp(- R * T) 73 74 callResult = (S * cndd1 - X * expRT * cndd2) 75 putResult = (X * expRT * (1.0 - cndd2) - S * (1.0 - cndd1)) 76 return callResult, putResult 77 78 79def blackscholes_scalar(callResult, putResult, stockPrice, optionStrike, 80 optionYears, Riskfree, Volatility): 81 S = stockPrice 82 X = optionStrike 83 T = optionYears 84 R = Riskfree 85 V = Volatility 86 for i in range(len(S)): 87 sqrtT = math.sqrt(T[i]) 88 d1 = (math.log(S[i] / X[i]) + (R + 0.5 * V * V) * T[i]) / (V * sqrtT) 89 d2 = d1 - V * sqrtT 90 cndd1 = cnd(d1) 91 cndd2 = cnd(d2) 92 93 expRT = math.exp((-1. * R) * T[i]) 94 callResult[i] = (S[i] * cndd1 - X[i] * expRT * cndd2) 95 putResult[i] = (X[i] * expRT * (1.0 - cndd2) - S[i] * (1.0 - cndd1)) 96 97 98def blackscholes_scalar_jitted(callResult, putResult, stockPrice, optionStrike, 99 optionYears, Riskfree, Volatility): 100 S = stockPrice 101 X = optionStrike 102 T = optionYears 103 R = Riskfree 104 V = Volatility 105 for i in range(len(S)): 106 sqrtT = math.sqrt(T[i]) 107 d1 = (math.log(S[i] / X[i]) + (R + 0.5 * V * V) * T[i]) / (V * sqrtT) 108 d2 = d1 - V * sqrtT 109 cndd1 = cnd_jitted(d1) 110 cndd2 = cnd_jitted(d2) 111 112 expRT = math.exp((-1. * R) * T[i]) 113 callResult[i] = (S[i] * cndd1 - X[i] * expRT * cndd2) 114 putResult[i] = (X[i] * expRT * (1.0 - cndd2) - S[i] * (1.0 - cndd1)) 115 116 117def randfloat(rand_var, low, high): 118 return (1.0 - rand_var) * low + rand_var * high 119 120 121class TestBlackScholes(TestCase): 122 def test_array_expr(self): 123 flags = Flags() 124 flags.set("enable_pyobject") 125 126 global cnd_array_jitted 127 scalty = types.float64 128 arrty = types.Array(scalty, 1, 'C') 129 cr1 = compile_isolated(cnd_array, args=(arrty,), flags=flags) 130 cnd_array_jitted = cr1.entry_point 131 cr2 = compile_isolated(blackscholes_arrayexpr_jitted, 132 args=(arrty, arrty, arrty, scalty, scalty), 133 flags=flags) 134 jitted_bs = cr2.entry_point 135 136 OPT_N = 400 137 iterations = 10 138 139 140 stockPrice = randfloat(self.random.random_sample(OPT_N), 5.0, 30.0) 141 optionStrike = randfloat(self.random.random_sample(OPT_N), 1.0, 100.0) 142 optionYears = randfloat(self.random.random_sample(OPT_N), 0.25, 10.0) 143 144 args = stockPrice, optionStrike, optionYears, RISKFREE, VOLATILITY 145 146 callResultGold, putResultGold = blackscholes_arrayexpr(*args) 147 callResultNumba, putResultNumba = jitted_bs(*args) 148 149 delta = np.abs(callResultGold - callResultNumba) 150 L1norm = delta.sum() / np.abs(callResultGold).sum() 151 print("L1 norm: %E" % L1norm) 152 print("Max absolute error: %E" % delta.max()) 153 self.assertEqual(delta.max(), 0) 154 155 def test_scalar(self): 156 flags = Flags() 157 158 # Compile the inner function 159 global cnd_jitted 160 cr1 = compile_isolated(cnd, (types.float64,)) 161 cnd_jitted = cr1.entry_point 162 # Manually type the compiled function for calling into 163 tyctx = cr1.typing_context 164 ctx = cr1.target_context 165 signature = typing.make_concrete_template("cnd_jitted", cnd_jitted, 166 [cr1.signature]) 167 tyctx.insert_user_function(cnd_jitted, signature) 168 169 # Compile the outer function 170 array = types.Array(types.float64, 1, 'C') 171 argtys = (array,) * 5 + (types.float64, types.float64) 172 cr2 = compile_extra(tyctx, ctx, blackscholes_scalar_jitted, 173 args=argtys, return_type=None, flags=flags, 174 locals={}) 175 jitted_bs = cr2.entry_point 176 177 OPT_N = 400 178 iterations = 10 179 180 callResultGold = np.zeros(OPT_N) 181 putResultGold = np.zeros(OPT_N) 182 183 callResultNumba = np.zeros(OPT_N) 184 putResultNumba = np.zeros(OPT_N) 185 186 stockPrice = randfloat(self.random.random_sample(OPT_N), 5.0, 30.0) 187 optionStrike = randfloat(self.random.random_sample(OPT_N), 1.0, 100.0) 188 optionYears = randfloat(self.random.random_sample(OPT_N), 0.25, 10.0) 189 190 args = stockPrice, optionStrike, optionYears, RISKFREE, VOLATILITY 191 192 blackscholes_scalar(callResultGold, putResultGold, *args) 193 jitted_bs(callResultNumba, putResultNumba, *args) 194 195 delta = np.abs(callResultGold - callResultNumba) 196 L1norm = delta.sum() / np.abs(callResultGold).sum() 197 print("L1 norm: %E" % L1norm) 198 print("Max absolute error: %E" % delta.max()) 199 self.assertAlmostEqual(delta.max(), 0) 200 201 202if __name__ == "__main__": 203 unittest.main() 204