## 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}$$

$$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 = """<table>
<tr>
<th>Greek</th>
<th>Call</th>
<th>Put</th>
</tr>
<tr>
<td>Delta $\Delta$</td>
<td>$N(d_1)$</td>
<td>$-N(-d_1)$</td>
</tr>
<tr>
<td>Theta $\Theta$</td>
<td>$\\frac{\sigma SN'(d_1)}{2\sqrt{T-t}} - rKe^{-r(T-t)} N(d_2)$</td>
<td>$\\frac{-\sigma SN'(d_1)}{2\sqrt{T-t}} + rKe^{-r(T-t)} N(d_2)$</td>
</tr>
<tr>
<td>Gamma $\Gamma$</td>
<td>$\\frac{N'(d_1)}{S\sigma \sqrt{T}}$</td>
<td>$\\frac{N'(d_1)}{S\sigma \sqrt{T}}$</td>
</tr>
<tr>
<td>Vega $v$</td>
<td>$S_0 N'(d_1)\sqrt{T}$</td>
<td>$S_0 N'(d_1)\sqrt{T}$</td>
</tr>
<tr>
<td>Rho $\\rho$</td>
<td>$TKe^{-r(T)} N(d_2)$</td>
<td>$-TKe^{-r(T)} N(-d_2)$</td>
</tr>
</table>"""

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 = """<table>
<tr>
<th>Greek</th>
<th>Call</th>
<th>Put</th>
</tr>
<tr>
<th>Delta $\Delta$</th>
<th>$e^{-qT}N(d_1)$</th>
<th>$-e^{-qT}N(-d_1)$</th>
</tr>
<tr>
<th>Theta $\Theta$</th>
<th>$-e^{-qT} \\frac{SN(d_1)\sigma}{2\sqrt{T}} - rKe^{-rT}N(d_2) + qSe^{-qT}N(d_1)$</th>
<th>$-e^{-qT} \\frac{SN(d_1)\sigma}{2\sqrt{T}} + rKe^{-rT}N(-d_2) - qSe^{-qT}N(-d_1)$</th>
</tr>
<tr>
<th>Gamma $\Gamma$</th>
<th>$e^{-qT}\\frac{N(d_1)}{S\sigma \sqrt{T}}$</th>
<th>$e^{-qT}\\frac{N(d_1)}{S\sigma \sqrt{T}}$</th>
</tr>
<tr>
<th>Vega $v$</th>
<th>$\\frac{1}{\sqrt{2\pi}} Se^{-q(T-t)} e^{-\\frac{d_1^2}{2}} \sqrt{T-t}$</th>
<th>$\\frac{1}{\sqrt{2\pi}} Se^{-q(T-t)} e^{-\\frac{d_1^2}{2}} \sqrt{T-t}$</th>
</tr>
<tr>
<th>Rho $\\rho$</th>
<th>$KTe^{-rT}N(d_2)$</th>
<th>$-KTe^{-rT}N(-d_2)$</th>
</tr>
</table>
"""

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