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.


```python
from IPython.display import HTML
```


```python
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>"""
```


```python
greeks = HTML(table); greeks
```




<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>



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

### Python Implementation


```python
import numpy as np
import scipy.stats as si
import sympy as sy
import sympy.statistics as systats
```

#### Delta $\Delta$


```python
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
```


```python
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
```


```python
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
```


```python
delta(100, 50, 1, 0.05, 0.25, option = 'put')
```




    -0.00097550992808483267



#### Theta $\Theta$


```python
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
```


```python
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
```


```python
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
```


```python
theta(110, 100, 2, 0.05, 0.25, option = 'put')
```




    -1.4189854326915405



#### Gamma


```python
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
```


```python
gamma(110, 100, 1, 0.05, 0.25)
```




    0.013760721161229799



#### Vega


```python
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
```


```python
vega(110, 105, 100, 1, 0.05, 0.25)
```




    39.734082353051043



#### Rho


```python
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
```


```python
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
```


```python
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


```python
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>
"""
```


```python
greeks_div = HTML(table_div); greeks_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>




### Python Implementation

#### Delta $\Delta$


```python
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
```


```python
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
```


```python
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$


```python
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
```


```python
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
```


```python
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$


```python
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$


```python
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$


```python
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
```


```python
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
```


```python
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
```

## References

[Sherbin, A. (2015). How to price and trade options: identify, analyze, and execute the best trade probabilities. 
    Hoboken, NJ: John Wiley & Sons, Inc.](https://amzn.to/37ajBnM)

[Ursone, P. (2015). How to calculate options prices and their Greeks: exploring the Black Scholes model from Delta 
    to Vega. Chichester: Wiley.](https://amzn.to/2UzXDrD)

c7dc847278371c896631a4f58ee38e020667030b
ipynb
Jupyter Notebook
content/posts/The Greeks.ipynb
aschleg/aaronschlegel.me
2f2e143218445da0b6298671c67f9c4afa055d59
MIT
2018-02-19T00:18:17.000Z
2020-01-17T15:11:31.000Z
content/posts/The Greeks.ipynb
aschleg/aaronschlegel.me
2f2e143218445da0b6298671c67f9c4afa055d59
MIT
2018-01-26T00:11:30.000Z
2018-01-26T00:11:30.000Z
content/posts/The Greeks.ipynb
aschleg/aaronschlegel.me
2f2e143218445da0b6298671c67f9c4afa055d59
MIT
Qwen/Qwen-72B
1. YES
2. YES