From 5d2fd34b3995bb429598bd384b4e4b791f8ef412 Mon Sep 17 00:00:00 2001
From: Justin Lars Kirkby <33753270+jkirkby3@users.noreply.github.com>
Date: Sun, 21 Apr 2024 17:12:36 -0400
Subject: [PATCH] Update black_scholes.py

Add gamma, theta, rho formulas
---
 fypy/pricing/analytical/black_scholes.py | 70 ++++++++++++++++++++++++
 1 file changed, 70 insertions(+)

diff --git a/fypy/pricing/analytical/black_scholes.py b/fypy/pricing/analytical/black_scholes.py
index 7080ea9..ed80c6a 100644
--- a/fypy/pricing/analytical/black_scholes.py
+++ b/fypy/pricing/analytical/black_scholes.py
@@ -59,6 +59,76 @@ def black76_price_strikes(F: float,
     return prices
 
 
+def black76_gamma(F: float,
+                  K: Union[float, np.ndarray],
+                  vol: Union[float, np.ndarray],
+                  disc: float,
+                  T: float) -> Union[float, np.ndarray]:
+    """
+    Gamma(s) for strike(s)
+    :param F: float, forward price
+    :param K: float or array, the Strike(s)
+    :param vol: float or array, the Volatility(ies) ... if float, all strikes get same vol, else a vol smile
+    :param disc: float, the discount factor, e.g. 0.99
+    :param T: float, time to maturity of option
+    :return: float or np.ndarray, same shape as strikes
+    """
+    vol_st = vol * np.sqrt(T)
+    d_1 = (np.log(F / K) + 0.5 * vol_st ** 2) / vol_st
+    return disc * norm.pdf(d_1) / (F * vol_st)
+
+
+def black76_theta(F: float,
+                  K: Union[float, np.ndarray],
+                  is_call: bool,
+                  vol: Union[float, np.ndarray],
+                  disc: float,
+                  T: float) -> Union[float, np.ndarray]:
+    """
+    Theta(s) for strike(s)
+    :param F: float, forward price
+    :param K: float or array, the Strike(s)
+    :param is_call: bool, determines if ALL strikes are call or all are put
+    :param vol: float or array, the Volatility(ies) ... if float, all strikes get same vol, else a vol smile
+    :param disc: float, the discount factor, e.g. 0.99
+    :param T: float, time to maturity of option
+    :return: float or np.ndarray, same shape as strikes
+    """
+    sqt = np.sqrt(T)
+    vol_st = vol * sqt
+    d_1 = (np.log(F / K) + 0.5 * vol_st ** 2) / vol_st
+    d_2 = d_1 - vol_st
+
+    term1 = - F * disc * norm.pdf(d_1) * vol / (2 * sqt)
+    r = - np.log(disc) / T
+
+    phi = 1 if is_call else -1
+    return term1 - phi * r * K * disc * norm.cdf(phi * d_2) + phi * r * F * disc * norm.cdf(phi * d_1)
+
+
+def black76_rho(F: float,
+                K: Union[float, np.ndarray],
+                is_call: bool,
+                vol: Union[float, np.ndarray],
+                disc: float,
+                T: float) -> Union[float, np.ndarray]:
+    """
+    Rhos(s) for strike(s)
+    :param F: float, forward price
+    :param K: float or array, the Strike(s)
+    :param is_call: bool, determines if ALL strikes are call or all are put
+    :param vol: float or array, the Volatility(ies) ... if float, all strikes get same vol, else a vol smile
+    :param disc: float, the discount factor, e.g. 0.99
+    :param T: float, time to maturity of option
+    :return: float or np.ndarray, same shape as strikes
+    """
+    vol_st = vol * np.sqrt(T)
+    d_1 = (np.log(F / K) + 0.5 * vol_st ** 2) / vol_st
+    d_2 = d_1 - vol_st
+    phi = 1 if is_call else -1
+    return -T * disc * phi * (F * norm.cdf(phi * d_1) - K * norm.cdf(phi * d_2))
+
+
 def black76_vega(F: float,
                  K: Union[float, np.ndarray],
                  vol: Union[float, np.ndarray],