closed-form solutions to Asian average strike and price options + tes…

…ts, moving Black-Sscholes formulae from AF example to `utils.financial_formulae`, support for dividend yield in B-S formulae + test (#497)

examples/PyMPDATA_examples/Arabas_and_Farhat_2020/setup1_european_corridor.py

@@ -1,5 +1,5 @@
 import numpy as np
-import PyMPDATA_examples.Arabas_and_Farhat_2020.Black_Scholes_1973 as BS73
+import PyMPDATA_examples.utils.financial_formulae.Black_Scholes_1973 as BS73
 from pystrict import strict
 
 

examples/PyMPDATA_examples/Arabas_and_Farhat_2020/setup2_american_put.py

@@ -1,6 +1,6 @@
 import numpy as np
-import PyMPDATA_examples.Arabas_and_Farhat_2020.Bjerksund_and_Stensland_1993 as BS93
-import PyMPDATA_examples.Arabas_and_Farhat_2020.Black_Scholes_1973 as BS73
+import PyMPDATA_examples.utils.financial_formulae.Bjerksund_and_Stensland_1993 as BS93
+import PyMPDATA_examples.utils.financial_formulae.Black_Scholes_1973 as BS73
 from pystrict import strict
 
 

examples/PyMPDATA_examples/utils/financial_formulae/Bjerksund_and_Stensland_1993.py

@@ -0,0 +1,62 @@
+import numpy as np
+import PyMPDATA_examples.utils.financial_formulae.Black_Scholes_1973 as BS
+
+
+def _phi(
+    S: [np.ndarray, float],
+    gamma: float,
+    H: float,
+    I: float,
+    r: float,
+    b: float,
+    var: float,
+    T: float,
+):
+    lmbd = (-r + gamma * b + 0.5 * gamma * (gamma - 1) * var) * T
+    d = -(np.log(S / H) + (b + (gamma - 0.5) * var) * T) / np.sqrt(var * T)
+    kappa = 2 * b / var + (2 * gamma - 1)
+    return (
+        np.exp(lmbd)
+        * np.power(S, gamma)
+        * (
+            BS.N(d)
+            - pow((I / S), kappa) * BS.N(d - 2 * np.log(I / S) / np.sqrt(var * T))
+        )
+    )
+
+
+def c_amer(
+    S: [np.ndarray, float],
+    K: [float, np.ndarray],
+    T: float,
+    r: float,
+    b: float,
+    sgma: float,
+):
+    if b >= r:
+        return BS.c_euro(S, K=K, T=T, r=r, b=b, sgma=sgma)
+
+    var = sgma * sgma
+    beta = (0.5 - b / var) + np.sqrt(pow((b / var - 0.5), 2) + 2 * r / var)
+    BInf = beta / (beta - 1) * K
+    B0 = np.maximum(K, r / (r - b) * K)
+    ht = -(b * T + 2 * sgma * np.sqrt(T)) * B0 / (BInf - B0)
+    I = B0 + (BInf - B0) * (1 - np.exp(ht))
+    alpha = (I - K) * pow(I, -beta)
+
+    return np.where(
+        S >= I,
+        S - K,
+        alpha * np.power(S, beta)
+        + (
+            -alpha * _phi(S, gamma=beta, H=I, I=I, r=r, b=b, var=var, T=T)
+            + _phi(S, gamma=1, H=I, I=I, r=r, b=b, var=var, T=T)
+            - _phi(S, gamma=1, H=K, I=I, r=r, b=b, var=var, T=T)
+            - K * _phi(S, gamma=0, H=I, I=I, r=r, b=b, var=var, T=T)
+            + K * _phi(S, gamma=0, H=K, I=I, r=r, b=b, var=var, T=T)
+        ),
+    )
+
+
+def p_amer(S: [np.ndarray, float], K: float, T: float, r: float, b: float, sgma: float):
+    return c_amer(K, K=S, T=T, r=r - b, b=-b, sgma=sgma)

examples/PyMPDATA_examples/utils/financial_formulae/Black_Scholes_1973.py

@@ -0,0 +1,32 @@
+import numpy as np
+from scipy.special import erf  # pylint: disable=no-name-in-module
+
+
+def N(x: float):
+    return (1 + erf(x / np.sqrt(2))) / 2
+
+
+def c_euro(S: np.ndarray, K: float, T: float, r: float, b: float, sgma: float):
+    d1 = (np.log(S / K) + (b + sgma * sgma / 2) * T) / sgma / np.sqrt(T)
+    d2 = d1 - sgma * np.sqrt(T)
+    return S * np.exp(b - r) * N(d1) - K * np.exp(-r * T) * N(d2)
+
+
+def p_euro(S: np.ndarray, K: float, T: float, r: float, b: float, sgma: float):
+    d1 = (np.log(S / K) + (b + sgma * sgma / 2) * T) / sgma / np.sqrt(T)
+    d2 = d1 - sgma * np.sqrt(T)
+    return K * np.exp(-r * T) * N(-d2) - S * np.exp((b - r) * T) * N(-d1)
+
+
+def c_euro_with_dividend(
+    S: np.ndarray, K: float, T: float, r: float, sgma: float, dividend_yield: float
+):
+    b = r - dividend_yield
+    return c_euro(S, K, T, r, b, sgma)
+
+
+def p_euro_with_dividend(
+    S: np.ndarray, K: float, T: float, r: float, sgma: float, dividend_yield: float
+):
+    b = r - dividend_yield
+    return p_euro(S, K, T, r, b, sgma)

examples/PyMPDATA_examples/utils/financial_formulae/asian_option.py

@@ -0,0 +1,60 @@
+# pylint: disable=line-too-long
+"""
+Closed-forms for geometric Asian options are taken from:
+[Derivatives Markets Appendix 19A](https://media.pearsoncmg.com/ph/bp/bridgepages/teamsite/mcdonald/McDonald-web-19-A.pdf)
+"""
+
+import numpy as np
+
+from .Black_Scholes_1973 import c_euro_with_dividend, p_euro_with_dividend
+
+# for fun in (c_euro_with_dividend, p_euro_with_dividend):
+#     locals()['geometric_asian_average_price_'+fun.__name__[0]] = lambda **args: fun(
+#         **{key:value for key, value in args.items() if key not in ('sgma', 'dividend_yield')},
+#         sgma=args['sgma']/np.sqrt(3),
+#         dividend_yield=0.5*(args['r'] + args['dividend_yield'] + args['sgma']**2/6),
+#     )
+
+
+def geometric_asian_average_price_c(S, K, T, r, sgma, dividend_yield):
+    return c_euro_with_dividend(
+        S=S,
+        K=K,
+        T=T,
+        r=r,
+        sgma=sgma / np.sqrt(3),
+        dividend_yield=0.5 * (r + dividend_yield + sgma**2 / 6),
+    )
+
+
+def geometric_asian_average_price_p(S, K, T, r, sgma, dividend_yield):
+    return p_euro_with_dividend(
+        S=S,
+        K=K,
+        T=T,
+        r=r,
+        sgma=sgma / np.sqrt(3),
+        dividend_yield=0.5 * (r + dividend_yield + sgma**2 / 6),
+    )
+
+
+def geometric_asian_average_strike_c(S, K, T, r, sgma, dividend_yield):
+    return c_euro_with_dividend(
+        S=S,
+        K=K,
+        T=T,
+        dividend_yield=dividend_yield,
+        sgma=sgma * np.sqrt(T / 3),
+        r=0.5 * (r + dividend_yield + sgma**2 / 6),
+    )
+
+
+def geometric_asian_average_strike_p(S, K, T, r, sgma, dividend_yield):
+    return p_euro_with_dividend(
+        S=S,
+        K=K,
+        T=T,
+        dividend_yield=dividend_yield,
+        sgma=sgma * np.sqrt(T / 3),
+        r=0.5 * (r + dividend_yield + sgma**2 / 6),
+    )

tests/smoke_tests/utils/test_financial_formulae.py

@@ -0,0 +1,62 @@
+# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring
+# pylint: disable=too-many-arguments,invalid-name
+import numpy as np
+import pytest
+from PyMPDATA_examples.utils.financial_formulae import Black_Scholes_1973 as BS73
+from PyMPDATA_examples.utils.financial_formulae import asian_option as AO
+
+
+class TestFinancialFormulae:
+    @staticmethod
+    @pytest.mark.parametrize(
+        "funs",
+        (
+            {"normal": BS73.c_euro, "with_dividend": BS73.c_euro_with_dividend},
+            {"normal": BS73.p_euro, "with_dividend": BS73.p_euro_with_dividend},
+        ),
+    )
+    @pytest.mark.parametrize("price", (np.array([95, 100, 105]),))
+    @pytest.mark.parametrize("strike", (100, 10))
+    @pytest.mark.parametrize("time_to_maturity", (1, 0.5))
+    @pytest.mark.parametrize("risk_free_rate", (0.05, 0.001))
+    @pytest.mark.parametrize("volatility", (0.2, 0.5))
+    @pytest.mark.parametrize("dividend_yield", (0.02, 0))
+    def test_black_scholes_with_dividend(
+        funs: dict,
+        price,
+        strike,
+        time_to_maturity,
+        risk_free_rate,
+        volatility,
+        dividend_yield,
+    ):
+        common_args = {
+            "S": price,
+            "K": strike,
+            "T": time_to_maturity,
+            "sgma": volatility,
+            "r": risk_free_rate,
+        }
+        price_dividend = funs["with_dividend"](
+            dividend_yield=dividend_yield, **common_args
+        )
+        price_normal = funs["normal"](b=risk_free_rate - dividend_yield, **common_args)
+        assert np.allclose(price_dividend, price_normal)
+
+    @staticmethod
+    @pytest.mark.parametrize(
+        "fun, expected_value",
+        (
+            (AO.geometric_asian_average_price_c, 3.246),
+            (AO.geometric_asian_average_price_p, 2.026),
+            (AO.geometric_asian_average_strike_c, 3.725),
+            (AO.geometric_asian_average_strike_p, 1.869),
+        ),
+    )
+    def test_asian_geometric_average(fun: callable, expected_value):
+        """
+        Analytic results are taken from [Derivatives Markets](
+        https://faculty.ksu.edu.sa/sites/default/files/derivatives_markets_3e_0.pdf) page 413
+        """
+        price = fun(S=40, K=40, T=1, r=0.08, sgma=0.3, dividend_yield=0)
+        assert np.allclose(price, expected_value, rtol=1e-3)