Measuring Sensitivity to Derivatives Pricing Changes with the "Greeks" and Python

The Greeks are used as risk measures that represent how sensitive the price of derivatives are to change. This is useful as risks can be treated in isolation and thus allows for tuning in a portfolio to reach a desired level of risk. The values are called 'the Greeks' as they are denoted by Greek letters. Each will be presented in turn as an introduction:

  • Delta ($\Delta$) is the first derivative of the spot price $S, \frac{\partial V}{\partial S}$

    • Denotes the rate of change of a portfolio's value
  • Gamma ($\Gamma$) is the second derivative of the spot price $S, \frac{\partial^2 V}{\partial S^2}$

    • Rate of change of the Delta
  • Theta ($\Theta$) is the first derivative of time $t, \frac{\partial V}{\partial t}$

    • Rate of change of the portfolio's value with respect to time $t$ (Not the maturity time $T$)
  • Rho ($\rho$) is the first derivative with respect to the risk-free rate $r, \frac{\partial V}{\partial r}$

    • Denotes the rate of change of the portfolio's value with respect to the risk-free interest rate
  • Vega is the first deriative of volatility $\sigma, \frac{\partial V}{\partial \sigma}$

    • Rate of change of value with respect to $\sigma$

General Derivations to Aid in Greek Computation

Black-Scholes Formula for a Call Option

$$C(S,t) = SN(d_1) - Ke^{-r(T - t)} N(d_2)$$

Probability Density Function of Standard Normal Variable

The standard normal variable will come up often in the derivation of the Greeks, so it is important to note this prior to proceeding. It is denoted by $z$ and is stated as the following:

$$f(x) = \frac{1}{\sqrt{2\pi}} e^{-\frac{x^2}{2}}$$

Derive $S \cdot N'(d_1) = Ke^{-r(T-t)}N'(d_2)$

We know the probability density function as presented above and that:

$$d_1 = \frac{\ln \left(\frac{S}{K} \right) + \left(r + \frac{\sigma^2}{2} \right)(T)}{\sigma \sqrt{T}} + \frac{\sigma \sqrt{T}}{2}$$ $$d_2 = d_1 - \sigma \sqrt{T - t} = \frac{\ln \left(\frac{S}{K} \right) + \left(r - \frac{\sigma^2}{2}\right)(T)}{\sigma \sqrt{T}} - \frac{\sigma \sqrt{T}}{2}$$

Therefore,

$$N'(d_1) = N'(d2 + \sigma \sqrt{T})$$ $$N'(d_1) = \frac{1}{\sqrt{2\pi}} e^{\frac{-d_1^2}{2}}$$ $$N'(d_2) = \frac{1}{\sqrt{2\pi}} e^{\frac{-d_2^2}{2}}$$

To prove the derivation above, insert the equivalent of $N'(d_1)$ and $N'(d_2)$ as shown above.

$$Se^{-\frac{d_1^2}{2}} = Ke^{-r(T-t)} e^{-\frac{d_2^2}{2}}$$

Recall $exp(ln(x)) = x$.

$$S = Ke^{-r(T-t)} \frac{S}{K}e^{r(T-t)}$$

$$S = Ke^{-r(T-t)} exp \left(ln \left(\frac{S}{K} \right) + r(T-t) \right)$$

We can rearrange $d_1$ above to get:

$$ln \left(\frac{S}{K} \right) + r(T-t) = d_1 \sigma \sqrt{T-t} - \frac{\sigma^2 (T-t)}{2}$$

Which then leads us to:

$$S = Ke^{-r(T-t)} exp \left(d_1 \sigma \sqrt{T-t} - \frac{\sigma^2 (T-t)}{2} \right)$$

And knowing $d_2 = d_1 - \sigma \sqrt{T-t}$

$$Se^{-\frac{d_1^2}{2}} = Ke^{-r(T-t)} exp \left(-\frac{d_1^2}{2} + d_1 \sigma \sqrt{T-t} - \frac{\sigma^2(T-t)}{2} \right)$$

$$Se^{-\frac{d_1^2}{2}} = Ke^{-r(T-t)} exp \left(-\frac{\left(d_1 - \sigma \sqrt{T-t} \right)^2}{2} \right)$$

$$Se^{-\frac{d_1^2}{2}} = Ke^{-r(T-t)} e^{-\frac{d_2^2}{2}}$$

Which translates to:

$$SN'(d_1) = Ke^{-r(T-t)} N'(d_2)$$

Taking the derivative with respect to the spot price $S$ of $d_1$ and $d_2$, we see:

$$\frac{\partial d_1}{\partial S} = \frac{\partial d_2}{\partial S} = \frac{1}{\sigma S \sqrt{T-t}}$$

The last two formulas above will be used frequently in the derivation of the Greeks.

Greeks of Non-Dividend Paying Assets

The implementation of the Greeks will be introduced for both non-dividend and dividend paying assets.

In [3]:
from IPython.display import HTML
In [16]:
table = """
Greek Call Put
Delta $\Delta$ $N(d_1)$ $-N(-d_1)$
Theta $\Theta$ $\\frac{\sigma SN'(d_1)}{2\sqrt{T-t}} - rKe^{-r(T-t)} N(d_2)$ $\\frac{-\sigma SN'(d_1)}{2\sqrt{T-t}} + rKe^{-r(T-t)} N(d_2)$
Gamma $\Gamma$ $\\frac{N'(d_1)}{S\sigma \sqrt{T}}$ $\\frac{N'(d_1)}{S\sigma \sqrt{T}}$
Vega $v$ $S_0 N'(d_1)\sqrt{T}$ $S_0 N'(d_1)\sqrt{T}$
Rho $\\rho$ $TKe^{-r(T)} N(d_2)$ $-TKe^{-r(T)} N(-d_2)$
"""
In [17]:
greeks = HTML(table); greeks
Out[17]:
Greek Call Put
Delta $\Delta$ $N(d_1)$ $-N(-d_1)$
Theta $\Theta$ $\frac{\sigma SN'(d_1)}{2\sqrt{T-t}} - rKe^{-r(T-t)} N(d_2)$ $\frac{-\sigma SN'(d_1)}{2\sqrt{T-t}} + rKe^{-r(T-t)} N(d_2)$
Gamma $\Gamma$ $\frac{N'(d_1)}{S\sigma \sqrt{T}}$ $\frac{N'(d_1)}{S\sigma \sqrt{T}}$
Vega $v$ $S_0 N'(d_1)\sqrt{T}$ $S_0 N'(d_1)\sqrt{T}$
Rho $\rho$ $TKe^{-r(T)} N(d_2)$ $-TKe^{-r(T)} N(-d_2)$

Above is a table of all the derivations of the Greek calls and puts for non-dividend paying assets.

Python Implementation

In [13]:
import numpy as np
import scipy.stats as si
import sympy as sy
import sympy.statistics as systats

Delta $\Delta$

In [27]:
def delta_call(S, K, T, r, sigma):
    
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    delta_call = si.norm.cdf(d1, 0.0, 1.0)
    
    return delta_call
In [30]:
def delta_put(S, K, T, r, sigma):
    
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    delta_put = si.norm.cdf(-d1, 0.0, 1.0)
    
    return -delta_put
In [34]:
def delta(S, K, T, r, sigma, option = 'call'):
    
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    if option == 'call':
        result = si.norm.cdf(d1, 0.0, 1.0)
    if option == 'put':
        result = -si.norm.cdf(-d1, 0.0, 1.0)
        
    return result
In [35]:
delta(100, 50, 1, 0.05, 0.25, option = 'put')
Out[35]:
-0.00097550992808483267

Theta $\Theta$

In [46]:
def theta_call(S, K, T, r, sigma):
    
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = (np.log(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    prob_density = 1 / np.sqrt(2 * np.pi) * np.exp(-d1 ** 2 * 0.5)
    
    theta = (-sigma * S * prob_density) / (2 * np.sqrt(T)) - r * K * np.exp(-r * T) * si.norm.cdf(d2, 0.0, 1.0)
    
    return theta
In [50]:
def theta_put(S, K, T, r, sigma):
    
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = (np.log(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    prob_density = 1 / np.sqrt(2 * np.pi) * np.exp(-d1 ** 2 * 0.5)
    
    theta = (-sigma * S * prob_density) / (2 * np.sqrt(T)) + r * K * np.exp(-r * T) * si.norm.cdf(-d2, 0.0, 1.0)
    
    return theta
In [52]:
def theta(S, K, T, r, sigma, option = 'call'):
    
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = (np.log(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    prob_density = 1 / np.sqrt(2 * np.pi) * np.exp(-d1 ** 2 * 0.5)
    
    if option == 'call':
        theta = (-sigma * S * prob_density) / (2 * np.sqrt(T)) - r * K * np.exp(-r * T) * si.norm.cdf(d2, 0.0, 1.0)
    if option == 'put':    
        theta = (-sigma * S * prob_density) / (2 * np.sqrt(T)) + r * K * np.exp(-r * T) * si.norm.cdf(-d2, 0.0, 1.0)
    
    return theta
In [54]:
theta(110, 100, 2, 0.05, 0.25, option = 'put')
Out[54]:
-1.4189854326915405

Gamma

In [57]:
def gamma(S, K, T, r, sigma):
    
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    prob_density = 1 / np.sqrt(2 * np.pi) * np.exp(-d1 ** 2 * 0.5)
    
    gamma = prob_density / (S * sigma * np.sqrt(T))
    
    return gamma
In [58]:
gamma(110, 100, 1, 0.05, 0.25)
Out[58]:
0.013760721161229799

Vega

In [60]:
def vega(S, S0, K, T, r, sigma):
    
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    prob_density = 1 / np.sqrt(2 * np.pi) * np.exp(-d1 ** 2 * 0.5)
    
    vega = S0 * prob_density * np.sqrt(T)
    
    return vega
In [61]:
vega(110, 105, 100, 1, 0.05, 0.25)
Out[61]:
39.734082353051043

Rho

In [62]:
def rho_call(S, K, T, r, sigma):
    
    d2 = (np.log(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    rho = T * K * np.exp(-r * T) * si.norm.cdf(d2, 0.0, 1.0)
    
    return rho
In [64]:
def rho_put(S, K, T, r, sigma):
    
    d2 = (np.log(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    rho = -T * K * np.exp(-r * T) * si.norm.cdf(-d2, 0.0, 1.0)
    
    return rho
In [65]:
def rho(S, K, T, r, sigma, option = 'call'):
    
    d2 = (np.log(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    if option == 'call':
        rho = T * K * np.exp(-r * T) * si.norm.cdf(d2, 0.0, 1.0)
    if option == 'put':
        rho = -T * K * np.exp(-r * T) * si.norm.cdf(-d2, 0.0, 1.0)
        
    return rho

Greeks of Dividend Paying Assets

In [5]:
table_div = """
Greek Call Put
Delta $\Delta$ $e^{-qT}N(d_1)$ $-e^{-qT}N(-d_1)$
Theta $\Theta$ $-e^{-qT} \\frac{SN(d_1)\sigma}{2\sqrt{T}} - rKe^{-rT}N(d_2) + qSe^{-qT}N(d_1)$ $-e^{-qT} \\frac{SN(d_1)\sigma}{2\sqrt{T}} + rKe^{-rT}N(-d_2) - qSe^{-qT}N(-d_1)$
Gamma $\Gamma$ $e^{-qT}\\frac{N(d_1)}{S\sigma \sqrt{T}}$ $e^{-qT}\\frac{N(d_1)}{S\sigma \sqrt{T}}$
Vega $v$ $\\frac{1}{\sqrt{2\pi}} Se^{-q(T-t)} e^{-\\frac{d_1^2}{2}} \sqrt{T-t}$ $\\frac{1}{\sqrt{2\pi}} Se^{-q(T-t)} e^{-\\frac{d_1^2}{2}} \sqrt{T-t}$
Rho $\\rho$ $KTe^{-rT}N(d_2)$ $-KTe^{-rT}N(-d_2)$
"""
In [6]:
greeks_div = HTML(table_div); greeks_div
Out[6]:
Greek Call Put
Delta $\Delta$ $e^{-qT}N(d_1)$ $-e^{-qT}N(-d_1)$
Theta $\Theta$ $-e^{-qT} \frac{SN(d_1)\sigma}{2\sqrt{T}} - rKe^{-rT}N(d_2) + qSe^{-qT}N(d_1)$ $-e^{-qT} \frac{SN(d_1)\sigma}{2\sqrt{T}} + rKe^{-rT}N(-d_2) - qSe^{-qT}N(-d_1)$
Gamma $\Gamma$ $e^{-qT}\frac{N(d_1)}{S\sigma \sqrt{T}}$ $e^{-qT}\frac{N(d_1)}{S\sigma \sqrt{T}}$
Vega $v$ $\frac{1}{\sqrt{2\pi}} Se^{-q(T-t)} e^{-\frac{d_1^2}{2}} \sqrt{T-t}$ $\frac{1}{\sqrt{2\pi}} Se^{-q(T-t)} e^{-\frac{d_1^2}{2}} \sqrt{T-t}$
Rho $\rho$ $KTe^{-rT}N(d_2)$ $-KTe^{-rT}N(-d_2)$

Python Implementation

Delta $\Delta$

In [25]:
def delta_call_div(S, K, T, r, q, sigma):
    
    d1 = (np.log(S / K) + (r - q + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    delta = np.exp(-q * T) * si.norm.cdf(d1, 0.0, 1.0)
    
    return delta
In [35]:
def delta_put_div(S, K, T, r, q, sigma):
    
    d1 = (np.log(S / K) + (r - q + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    delta = -np.exp(-q * T) * si.norm.cdf(-d1, 0.0, 1.0)
    
    return delta
In [46]:
def delta_div(S, K, T, r, q, sigma, option = 'call'):
    
    d1 = (np.log(S / K) + (r - q + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    if option == 'call':
        delta = delta = np.exp(-q * T) * si.norm.cdf(d1, 0.0, 1.0)
    if option == 'put':
        delta = -np.exp(-q * T) * si.norm.cdf(-d1, 0.0, 1.0)
        
    return delta

Theta $\Theta$

In [52]:
def theta_call_div(S, K, T, r, q, sigma):
    
    d1 = (np.log(S / K) + (r - q + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = (np.log(S / K) + (r - q - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    theta = -np.exp(-q * T) * (S * si.norm.cdf(d1, 0.0, 1.0) * sigma) / (2 * np.sqrt(T)) - r * K * np.exp(-r * T) * si.norm.cdf(d2, 0.0, 1.0) + q * S * np.exp(-q * T) * si.norm.cdf(d1, 0.0, 1.0)
    
    return theta
In [56]:
def theta_put_div(S, K, T, r, q, sigma):

    d1 = (np.log(S / K) + (r - q + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = (np.log(S / K) + (r - q - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    theta = -np.exp(-q * T) * (S * si.norm.cdf(d1, 0.0, 1.0) * sigma) / (2 * np.sqrt(T)) + r * K * np.exp(-r * T) * si.norm.cdf(-d2, 0.0, 1.0) - q * S * np.exp(-q * T) * si.norm.cdf(-d1, 0.0, 1.0)
    
    return theta
In [57]:
def theta_div(S, K, T, r, q, sigma, option = 'call'):
    
    d1 = (np.log(S / K) + (r - q + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = (np.log(S / K) + (r - q - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    if option == 'call':
        theta = -np.exp(-q * T) * (S * si.norm.cdf(d1, 0.0, 1.0) * sigma) / (2 * np.sqrt(T)) - r * K * np.exp(-r * T) * si.norm.cdf(d2, 0.0, 1.0) + q * S * np.exp(-q * T) * si.norm.cdf(d1, 0.0, 1.0)
    if option == 'put':
        theta = -np.exp(-q * T) * (S * si.norm.cdf(d1, 0.0, 1.0) * sigma) / (2 * np.sqrt(T)) + r * K * np.exp(-r * T) * si.norm.cdf(-d2, 0.0, 1.0) - q * S * np.exp(-q * T) * si.norm.cdf(-d1, 0.0, 1.0)
        
    return theta

Gamma $\Gamma$

In [59]:
def gamma_div(S, K, T, r, q, sigma):
    
    d1 = (np.log(S / K) + (r - q + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    gamma = np.exp(-q * T) * si.norm.cdf(d1, 0.0, 1.0) / S * sigma * np.sqrt(T)
    
    return gamma

Vega $v$

In [11]:
def vega_div(S, K, T, r, q, sigma):

    d1 = (np.log(S / K) + (r - q + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    vega = 1 / np.sqrt(2 * np.pi) * S * np.exp(-q * T) * np.exp(-d1 ** 2 * 0.5) * np.sqrt(T)
    
    return vega

Rho $\rho$

In [16]:
def rho_call_div(S, K, T, r, q, sigma):
    
    d2 = (np.log(S / K) + (r - q - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    rho = K * T * np.exp(-r * T) * si.norm.cdf(d2, 0.0, 1.0)
    
    return rho
In [17]:
def rho_put_div(S, K, T, r, q, sigma):
    
    d2 = (np.log(S / K) + (r - q - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    rho = -K * T * np.exp(-r * T) * si.norm.cdf(-d2, 0.0, 1.0)
    
    return rho
In [18]:
def rho_div(S, K, T, r, q, sigma, option = 'call'):
    
    d2 = (np.log(S / K) + (r - q - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    if option == 'call':
        rho = K * T * np.exp(-r * T) * si.norm.cdf(d2, 0.0, 1.0)
    if option == 'put':
        rho = -K * T * np.exp(-r * T) * si.norm.cdf(-d2, 0.0, 1.0)
        
    return rho