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