增加美式期权

Author: sfl666

american_option.py

@@ -0,0 +1,176 @@
+"""
+Author: shifulin
+Email: [email protected]
+"""
+from math import sqrt, exp, inf
+import numpy as np
+
+
+def _call_price(s, k, sigma, r, t, steps=100):
+    r_ = exp(r * (t / steps))
+    r_reciprocal = 1.0 / r_
+    u = exp(sigma * sqrt(t / steps))
+    d = 1.0 / u
+    u_square = u ** 2
+    p_u = (r_ - d) / (u - d)
+    p_d = 1.0 - p_u
+    prices = np.zeros(steps + 1)
+    prices[0] = s * d ** steps
+    for i in range(1, steps + 1):
+        prices[i] = prices[i - 1] * u_square
+    values = np.zeros(steps + 1)
+    for i in range(steps + 1):
+        values[i] = max(0.0, prices[i] - k)
+    for j in range(steps, 0, -1):
+        for i in range(j):
+            values[i] = (p_u * values[i + 1] + p_d * values[i]) * r_reciprocal
+            prices[i] = d * prices[i + 1]
+            values[i] = max(values[i], prices[i] - k)
+    # print(values)
+    return values[0]
+
+
+def _put_price(s, k, sigma, r, t, steps=100):
+    r_ = exp(r * (t / steps))
+    r_reciprocal = 1.0 / r_
+    u = exp(sigma * sqrt(t / steps))
+    d = 1.0 / u
+    u_square = u ** 2
+    p_u = (r_ - d) / (u - d)
+    p_d = 1.0 - p_u
+    prices = np.zeros(steps + 1)
+    prices[0] = s * d ** steps
+    for i in range(1, steps + 1):
+        prices[i] = prices[i - 1] * u_square
+    values = np.zeros(steps + 1)
+    for i in range(steps + 1):
+        values[i] = max(0, k - prices[i])
+    for j in range(steps, 0, -1):
+        for i in range(0, j):
+            values[i] = (p_u * values[i + 1] + p_d * values[i]) * r_reciprocal
+            prices[i] = d * prices[i + 1]
+            values[i] = max(values[i], k - prices[i])
+    return values[0]
+
+
+def call_price(s, k, sigma, r, t, steps=100):
+    return (_call_price(s, k, sigma, r, t, steps) + _call_price(s, k, sigma, r, t, steps + 1)) / 2.0
+
+
+def put_price(s, k, sigma, r, t, steps=100):
+    return (_put_price(s, k, sigma, r, t, steps) + _put_price(s, k, sigma, r, t, steps + 1)) / 2.0
+
+
+def delta(s, k, sigma, r, t, option_type, steps=100):
+    if t == 0.0:
+        if s == k:
+            return {'Call': 0.5, 'Put': -0.5}[option_type]
+        elif s > k:
+            return {'Call': 1.0, 'Put': 0.0}[option_type]
+        else:
+            return {'Call': 0.0, 'Put': -1.0}[option_type]
+    else:
+        price_func = {'Call': call_price, 'Put': put_price}[option_type]
+        return (price_func(s + 0.01, k, sigma, r, t, steps=steps) -
+                price_func(s - 0.01, k, sigma, r, t, steps=steps)) * 50.0
+
+
+def gamma(s, k, sigma, r, t, option_type, steps=100):
+    if t == 0.0:
+        return inf if s == k else 0.0
+    price_func = {'Call': call_price, 'Put': put_price}[option_type]
+    return (price_func(s + 0.01, k, sigma, r, t, steps=steps) +
+            price_func(s + 0.01, k, sigma, r, t, steps=steps) -
+            price_func(s, k, sigma, r, t, steps=steps) * 2.0) * 10000.0
+
+
+def theta(s, k, sigma, r, t, option_type, steps=100):
+    price_func = {'Call': call_price, 'Put': put_price}[option_type]
+    t_unit = 1.0 / 365.0
+    if t <= t_unit:
+        return price_func(s, k, sigma, r, 0.0001, steps=steps) - \
+               price_func(s, k, sigma, r, t, steps=steps)
+    else:
+        return price_func(s, k, sigma, r, t - t_unit, steps=steps) - \
+               price_func(s, k, sigma, r, t, steps=steps)
+
+
+def vega(s, k, sigma, r, t, option_type, steps=100):
+    price_func = {'Call': call_price, 'Put': put_price}[option_type]
+    if sigma < 0.02:
+        return 0.0
+    else:
+        return (price_func(s, k, sigma + 0.01, r, t, steps=steps) -
+                price_func(s, k, sigma - 0.01, r, t, steps=steps)) * 50.0
+
+
+def rho(s, k, sigma, r, t, option_type, steps=100):
+    price_func = {'Call': call_price, 'Put': put_price}[option_type]
+    return (price_func(s, k, sigma, r + 0.001, t, steps=steps) -
+            price_func(s, k, sigma, r - 0.001, t, steps=steps)) * 500.0
+
+
+def call_iv(c, s, k, t, r=0.03, sigma_min=0.01, sigma_max=1.0, e=0.00001, steps=100):
+    sigma_mid = (sigma_min + sigma_max) / 2.0
+    call_min = call_price(s, k, sigma_min, r, t, steps)
+    call_max = call_price(s, k, sigma_max, r, t, steps)
+    call_mid = call_price(s, k, sigma_mid, r, t, steps)
+    diff = c - call_mid
+    if c <= call_min:
+        return sigma_min
+    elif c >= call_max:
+        return sigma_max
+    while abs(diff) > e:
+        if c > call_mid:
+            sigma_min = sigma_mid
+        else:
+            sigma_max = sigma_mid
+        sigma_mid = (sigma_min + sigma_max) / 2.0
+        call_mid = call_price(s, k, sigma_mid, r, t, steps)
+        diff = c - call_mid
+    # print(sigma_mid)
+    return sigma_mid
+
+
+def put_iv(c, s, k, t, r=0.03, sigma_min=0.01, sigma_max=1.0, e=0.00001, steps=100):
+    sigma_mid = (sigma_min + sigma_max) / 2.0
+    put_min = put_price(s, k, sigma_min, r, t, steps)
+    put_max = put_price(s, k, sigma_max, r, t, steps)
+    put_mid = put_price(s, k, sigma_mid, r, t, steps)
+    diff = c - put_mid
+    if c <= put_min:
+        return sigma_min
+    elif c >= put_max:
+        return sigma_max
+    while abs(diff) > e:
+        if c > put_mid:
+            sigma_min = sigma_mid
+        else:
+            sigma_max = sigma_mid
+        sigma_mid = (sigma_min + sigma_max) / 2.0
+        put_mid = put_price(s, k, sigma_mid, r, t, steps)
+        diff = c - put_mid
+    return sigma_mid
+
+
+def my_test():
+    import matplotlib.pyplot as plt
+    a = np.linspace(1.0 / 365.0, 2, 100)
+    yc, yp = [], []
+    for i in a:
+        yc.append(vega(6.0, 5.0, 0.25, 0.03, i, option_type='Call', steps=100))
+        yp.append(vega(6.0, 5.0, 0.25, 0.03, i, option_type='Put', steps=100))
+    plt.plot(yc)
+    plt.plot(yp)
+    plt.show()
+
+
+def my_test2():
+    # print(call_price(5.0, 5.0, 0.1, 0.03, 0.4))
+    # call_price(5.0, 5.0, 0.25, 0.03, 0.4, 99)
+    print(call_iv(0.138, 3.046, 3.1, 0.5, r=0.03, sigma_min=0.01, sigma_max=1.0, e=0.00001, steps=100))
+
+
+if __name__ == '__main__':
+    my_test2()
+

european_option.py

@@ -3,7 +3,7 @@
 Email: [email protected]
 """
 from math import log, sqrt, exp
-# import numpy as np
+import numpy as np
 from scipy.stats import norm
 
 
@@ -17,50 +17,44 @@
 #     d1 = (np.log(s / k) + (r + pow(sigma, 2) / 2.0) * t) / (sigma * np.sqrt(t))
 #     d2 = d1 - sigma * np.sqrt(t)
 #     return k * np.exp(-r * t) * norm.cdf(-d2) - s * norm.cdf(-d1)
-#
-#
-# def delta(s, k, sigma, r, t, option_type, position_type):
-#     d1 = (np.log(s / k) + (r + pow(sigma, 2) / 2.0) * t) / (sigma * np.sqrt(t))
-#     if option_type == 'Call':
-#         if position_type == 'long':
-#             return norm.cdf(d1)
-#         else:
-#             return -norm.cdf(d1)
-#     else:
-#         if position_type == 'long':
-#             return norm.cdf(d1) - 1.0
-#         else:
-#             return 1.0 - norm.cdf(d1)
-#
-#
-# def gamma(s, k, sigma, r, t):
-#     d1 = (np.log(s / k) + (r + pow(sigma, 2) / 2.0) * t) / (sigma * np.sqrt(t))
-#     return np.exp(-pow(d1, 2) / 2.0) / (s * sigma * np.sqrt(2.0 * np.pi * t))
-#
-#
-# def theta(s, k, sigma, r, t, option_type):
-#     d1 = (np.log(s / k) + (r + pow(sigma, 2) / 2.0) * t) / (sigma * np.sqrt(t))
-#     d2 = d1 - sigma * np.sqrt(t)
-#     theta_call = -(s * sigma * np.exp(-pow(d1, 2) / 2.0)) / (2.0 * np.sqrt(2.0 * np.pi * t)) - \
-#         r * k * np.exp(-r * t) * norm.cdf(d2)
-#     if option_type == 'Call':
-#         return theta_call
-#     else:
-#         return theta_call + r * k * np.exp(-r * t)
-#
-#
-# def vega(s, k, sigma, r, t):
-#     d1 = (np.log(s / k) + (r + pow(sigma, 2) / 2.0) * t) / (sigma * np.sqrt(t))
-#     return s * np.sqrt(t) * np.exp(-pow(d1, 2) / 2.0) / np.sqrt(2.0 * np.pi)
-#
-#
-# def rho(s, k, sigma, r, t, option_type):
-#     d1 = (np.log(s / k) + (r + pow(sigma, 2) / 2.0) * t) / (sigma * np.sqrt(t))
-#     d2 = d1 - sigma * np.sqrt(t)
-#     if option_type == 'Call':
-#         return k * t * np.exp(-r * t) * norm.cdf(d2)
-#     else:
-#         return -k * t * np.exp(-r * t) * norm.cdf(-d2)
+
+
+def delta(s, k, sigma, r, t, option_type):
+    d1 = (np.log(s / k) + (r + pow(sigma, 2) / 2.0) * t) / (sigma * np.sqrt(t))
+    if option_type == 'Call':
+        return norm.cdf(d1)
+    else:
+        return norm.cdf(d1) - 1.0
+
+
+def gamma(s, k, sigma, r, t):
+    d1 = (np.log(s / k) + (r + pow(sigma, 2) / 2.0) * t) / (sigma * np.sqrt(t))
+    return np.exp(-pow(d1, 2) / 2.0) / (s * sigma * np.sqrt(2.0 * np.pi * t))
+
+
+def theta(s, k, sigma, r, t, option_type):
+    d1 = (np.log(s / k) + (r + pow(sigma, 2) / 2.0) * t) / (sigma * np.sqrt(t))
+    d2 = d1 - sigma * np.sqrt(t)
+    theta_call = -(s * sigma * np.exp(-pow(d1, 2) / 2.0)) / (2.0 * np.sqrt(2.0 * np.pi * t)) - \
+        r * k * np.exp(-r * t) * norm.cdf(d2)
+    if option_type == 'Call':
+        return theta_call
+    else:
+        return theta_call + r * k * np.exp(-r * t)
+
+
+def vega(s, k, sigma, r, t):
+    d1 = (np.log(s / k) + (r + pow(sigma, 2) / 2.0) * t) / (sigma * np.sqrt(t))
+    return s * np.sqrt(t) * np.exp(-pow(d1, 2) / 2.0) / np.sqrt(2.0 * np.pi)
+
+
+def rho(s, k, sigma, r, t, option_type):
+    d1 = (np.log(s / k) + (r + pow(sigma, 2) / 2.0) * t) / (sigma * np.sqrt(t))
+    d2 = d1 - sigma * np.sqrt(t)
+    if option_type == 'Call':
+        return k * t * np.exp(-r * t) * norm.cdf(d2)
+    else:
+        return -k * t * np.exp(-r * t) * norm.cdf(-d2)
 
 
 def bs_call(s, k, sigma, r, t):
@@ -124,6 +118,18 @@ def my_test():
     call_iv(0.138, 3.046, 3.1, 0.5, r=0.03, sigma_min=0.01, sigma_max=1.0, e=0.000001)
 
 
+def my_test2():
+    import matplotlib.pyplot as plt
+    a = np.linspace(0, 0.8, 100)
+    yc, yp = [], []
+    for i in a:
+        yc.append(vega(6.0, 5.0, i, 0.03, 0.5))
+        yp.append(vega(6.0, 5.0, i, 0.03, 0.5))
+    plt.plot(yc)
+    plt.plot(yp)
+    plt.show()
+
+
 if __name__ == '__main__':
-    my_test()
+    my_test2()
 

historical_implied_volatility.py

@@ -14,7 +14,8 @@
 from sina_future_kline_api import get_future_day_kline
 from sina_commodity_option_api import get_option_kline as get_future_option_day_kline
 from sina_etf_option_api import get_option_day_kline as get_etf_option_day_kline
-from european_option import call_iv, put_iv
+import european_option
+import american_option
 
 
 ETF_SPOT_CODE = {
@@ -104,9 +105,9 @@ def align_kline(option_kline, spot_kline, window_size):
 
 def cal_historical_iv(option_kline, spot_kline, strike_price, expiry_date, r, option_type, exercise_type):
     if exercise_type == 'european':
-        iv_func = call_iv if option_type == 'call' else put_iv
+        iv_func = european_option.call_iv if option_type == 'call' else european_option.put_iv
     else:
-        pass
+        iv_func = american_option.call_iv if option_type == 'call' else american_option.put_iv
     x, y = [], []
     for option, spot in zip(option_kline, spot_kline):
         x.append(str(option[0]))
@@ -156,6 +157,8 @@ def main(option_code, spot_code, strike_price, expiry_date, option_type, exercis
 
 
 if __name__ == '__main__':
-    main('cu2003C51000', 'cu2003', 51000.0, '20200224', 'call', 'european', 5)
+    # main('cu2003C51000', 'cu2003', 51000.0, '20200224', 'call', 'european', 5)
     # main('au2004P340', 'au2004', 340.0, '20200325', 'put', 'european', 10)
     # main('10002062', '510050', 3.0, '20200122', 'put', 'european', 20)
+    # main('m2005C2800', 'm2005', 2800.0, '20200408', 'call', 'american', 5)
+    main('m2005P2700', 'm2005', 2700.0, '20200408', 'put', 'american', 10)

sina_commodity_option_api.py

@@ -17,6 +17,7 @@
     'ru': {'product': 'ru_o', 'exchange': 'shfe'},
     'cu': {'product': 'cu_o', 'exchange': 'shfe'},
     'au': {'product': 'au_o', 'exchange': 'shfe'},
+    'rm': {'product': 'rm', 'exchange': 'czce'},
 }
 URL_T_QUOTATION = "http://stock.finance.sina.com.cn/futures/api/openapi.php/OptionService.getOptionData?" \
                   "type=futures&product={product}&exchange={exchange}&pinzhong={code}"

volatility_surface.py

@@ -246,6 +246,7 @@ def main(cate, exchange, underlying, dividend=True, is_fit=True):
     ax_iv_sf_call = fig.add_subplot(gs[:2, :2], projection='3d')
     ax_iv_sf_call.view_init(ELEV, AZIM)
     # surf_call = ax_iv_sf_call.plot_surface(x, y, array(call_y2), rstride=1, cstride=1, cmap='rainbow')
+    # print(x.shape, y.shape, array(call_y2).shape)
     surf_call = ax_iv_sf_call.plot_wireframe(x, y, array(call_y2), rstride=1, cstride=1, cmap='rainbow')
     ax_iv_sf_call.set_yticklabels(dates_label)
     ax_iv_sf_call.set_xlabel('Strike Price')