From 9a5c376a9adccabf3ad3d8f2a261a063699ab9a9 Mon Sep 17 00:00:00 2001
From: simaki <shota.imaki.0801@gmail.com>
Date: Wed, 23 Mar 2022 17:50:59 +0900
Subject: [PATCH] ENH: Analytical BS European binary formulas (#437) (#553)

---
 pfhedge/nn/functional.py | 72 +++++++++++++++++++++++++---------------
 1 file changed, 46 insertions(+), 26 deletions(-)

diff --git a/pfhedge/nn/functional.py b/pfhedge/nn/functional.py
index b3643157..aa5008da 100644
--- a/pfhedge/nn/functional.py
+++ b/pfhedge/nn/functional.py
@@ -703,6 +703,21 @@ def bilerp(
     return torch.lerp(lerp1, lerp2, weight2)
 
 
+def _bs_theta_gamma_relation(gamma: Tensor, spot: Tensor, volatility: Tensor) -> Tensor:
+    # theta = (1/2) * vola^2 * spot^2 * gamma
+    # by Black-Scholes formula
+    return gamma * volatility.square() * spot.square() / 2
+
+
+def _bs_vega_gamma_relation(
+    gamma: Tensor, spot: Tensor, time_to_maturity: Tensor, volatility: Tensor
+) -> Tensor:
+    # vega = vola * spot^2 * time * gamma
+    # in Black-Scholes model
+    # See Chapter 5 Appendix A, Bergomi "Stochastic volatility modeling"
+    return gamma * volatility * spot.square() * time_to_maturity
+
+
 def bs_european_price(
     log_moneyness: Tensor,
     time_to_maturity: Tensor,
@@ -849,14 +864,13 @@ def bs_european_binary_gamma(
 
     See :func:`pfhedge.nn.BSEuropeanBinaryOption.gamma` for details.
     """
-    # TODO(simaki): Directly compute gamma.
-    return autogreek.gamma_from_delta(
-        bs_european_binary_delta,
-        log_moneyness=log_moneyness,
-        time_to_maturity=time_to_maturity,
-        volatility=volatility,
-        strike=strike,
-    )
+    s, t, v = broadcast_all(log_moneyness, time_to_maturity, volatility)
+    spot = s.exp() * strike
+
+    d2_tensor = d2(s, t, v)
+    w = volatility * time_to_maturity.square()
+
+    return -npdf(d2_tensor).div(w * spot.square()) * (1 + d2_tensor.div(w))
 
 
 def bs_european_binary_vega(
@@ -869,14 +883,16 @@ def bs_european_binary_vega(
 
     See :func:`pfhedge.nn.BSEuropeanBinaryOption.vega` for details.
     """
-    # TODO(simaki): Directly compute gamma.
-    return autogreek.vega(
-        bs_european_binary_price,
+    gamma = bs_european_binary_gamma(
         log_moneyness=log_moneyness,
         time_to_maturity=time_to_maturity,
         volatility=volatility,
         strike=strike,
     )
+    spot = log_moneyness.exp() * strike
+    return _bs_vega_gamma_relation(
+        gamma, spot=spot, time_to_maturity=time_to_maturity, volatility=volatility
+    )
 
 
 def bs_european_binary_theta(
@@ -889,14 +905,14 @@ def bs_european_binary_theta(
 
     See :func:`pfhedge.nn.BSEuropeanBinaryOption.theta` for details.
     """
-    # TODO(simaki): Directly compute theta.
-    return autogreek.theta(
-        bs_european_binary_price,
+    gamma = bs_european_binary_gamma(
         log_moneyness=log_moneyness,
         time_to_maturity=time_to_maturity,
         volatility=volatility,
         strike=strike,
     )
+    spot = log_moneyness.exp() * strike
+    return _bs_theta_gamma_relation(gamma, spot=spot, volatility=volatility)
 
 
 def bs_american_binary_price(
@@ -975,15 +991,17 @@ def bs_american_binary_vega(
 
     See :func:`pfhedge.nn.BSAmericanBinaryOption.vega` for details.
     """
-    # TODO(simaki): Compute analytically
-    return autogreek.vega(
-        bs_american_binary_price,
+    gamma = bs_american_binary_gamma(
         log_moneyness=log_moneyness,
         max_log_moneyness=max_log_moneyness,
         time_to_maturity=time_to_maturity,
         volatility=volatility,
         strike=strike,
     )
+    spot = log_moneyness.exp() * strike
+    return _bs_vega_gamma_relation(
+        gamma, spot=spot, time_to_maturity=time_to_maturity, volatility=volatility
+    )
 
 
 def bs_american_binary_theta(
@@ -997,15 +1015,15 @@ def bs_american_binary_theta(
 
     See :func:`pfhedge.nn.BSAmericanBinaryOption.theta` for details.
     """
-    # TODO(simaki): Compute analytically
-    return autogreek.theta(
-        bs_american_binary_price,
+    gamma = bs_american_binary_gamma(
         log_moneyness=log_moneyness,
         max_log_moneyness=max_log_moneyness,
         time_to_maturity=time_to_maturity,
         volatility=volatility,
         strike=strike,
     )
+    spot = log_moneyness.exp() * strike
+    return _bs_theta_gamma_relation(gamma, spot=spot, volatility=volatility)
 
 
 def bs_lookback_price(
@@ -1100,15 +1118,17 @@ def bs_lookback_vega(
 
     See :func:`pfhedge.nn.BSLookbackOption.vega` for details.
     """
-    # TODO(simaki): Calculate analytically
-    return autogreek.vega(
-        bs_lookback_price,
+    gamma = bs_lookback_gamma(
         log_moneyness=log_moneyness,
         max_log_moneyness=max_log_moneyness,
         time_to_maturity=time_to_maturity,
         volatility=volatility,
         strike=strike,
     )
+    spot = log_moneyness.exp() * strike
+    return _bs_vega_gamma_relation(
+        gamma, spot=spot, time_to_maturity=time_to_maturity, volatility=volatility
+    )
 
 
 def bs_lookback_theta(
@@ -1122,15 +1142,15 @@ def bs_lookback_theta(
 
     See :func:`pfhedge.nn.BSLookbackOption.theta` for details.
     """
-    # TODO(simaki): Calculate analytically
-    return autogreek.theta(
-        bs_lookback_price,
+    gamma = bs_lookback_gamma(
         log_moneyness=log_moneyness,
         max_log_moneyness=max_log_moneyness,
         time_to_maturity=time_to_maturity,
         volatility=volatility,
         strike=strike,
     )
+    spot = log_moneyness.exp() * strike
+    return _bs_theta_gamma_relation(gamma, spot=spot, volatility=volatility)
 
 
 def box_muller(