# Option-implied PDFs

http://econweb.rutgers.edu/mizrach/pubs/[42]-2010_Handbook.pdf

https://www.mathworks.com/company/newsletters/articles/estimating-option-implied-probability-distributions-for-asset-pricing.html

## Black-Scholes

$$\begin{align}
C &= N(d_1) S - N(d_2) Ke^{-rt} \\
P &= - N(-d_1)S + N(-d_2)Ke^{-rt} 
\end{align}
$$

where 
$$\begin{align}
d_1 &= \frac{1}{\sigma \sqrt{t}}\left[\ln \frac S K + (r + \frac{\sigma^2}{2})t\right] \\
d_2 &=  d_1 - \sigma \sqrt t
\end{align}
$$

We will use a simple Black-Scholes model in order to generate payoff diagrams.


```python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy
from scipy.stats import norm
import yfinance as yf
```


```python
def option_value(option_type, S, K, sigma, t=0, r=0):
    """
    Calculate the value of an option using the Black-Scholes model
    
    :param option_type: "call"/"c" or "put"/"p"
    :type option_type: str
    :param S: price of the underlying
    :type S: float
    :param K: strike price of option
    :type K: float
    :param sigma: input implied volatility
    :type sigma: float
    :param t: time to expiration
    :type t: float, optional
    :param r: risk-free rate
    :type r: float, optional
    """
    with np.errstate(divide='ignore'):
        d1 = np.divide(1, sigma * np.sqrt(t)) * (np.log(S/K) + (r+sigma**2 / 2) * t)
        d2 = d1 - sigma * np.sqrt(t)
    if option_type.lower() in {"c", "call"}:
        return np.multiply(norm.cdf(d1),S) - np.multiply(norm.cdf(d2), K * np.exp(-r * t))
    elif option_type.lower() in {"p", "put"}:
        return -np.multiply(norm.cdf(-d1), S) + np.multiply(norm.cdf(-d2), K * np.exp(-r * t))
```


```python
# Construction of a butterfly spread
S = np.linspace(50, 150, 1000)
C1 = option_value("c", S, 90, sigma=0.20)
C2 = -option_value("c", S, 100, sigma=0.20)
C3 = option_value("c", S, 110, sigma=0.20)
butterfly = C1 + 2 * C2 + C3
```


```python
# (Gross) payoff diagram
fig, (ax, ax1) = plt.subplots(1,2, figsize=(12,4), sharey=True)
ax.plot(S, C1, S, C2, S, C3)
ax.set_xlabel("Stock price at expiration")
ax.set_ylabel("Gross payoff")
ax.legend(["long call, $90 strike", "2x short call, $100 strike", "long call, $110 strike"], loc="best")

ax1.plot(S, butterfly, c="m")
ax1.legend(["long butterfly"], loc="upper left")
# plt.show();
plt.savefig("long_butterfly.png", dpi=200)
```

## Implied distributions from butterfly spreads

Using data from my broker, Firstrade, we will attempt to use the prices of butterfly spreads to infer the risk-neutral probability of a certain move in the underlying.


```python
# Read in call and put data from excel
calls = pd.read_excel("SPY_191020exp_290920.xlsx", sheet_name="call")
calls
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>strike</th>
      <th>lastPrice</th>
      <th>change</th>
      <th>bid</th>
      <th>ask</th>
      <th>BxA</th>
      <th>volume</th>
      <th>openInterest</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>170</td>
      <td>0.00</td>
      <td>0.0</td>
      <td>163.70</td>
      <td>164.06</td>
      <td>146 x 146</td>
      <td>0</td>
      <td>0</td>
    </tr>
    <tr>
      <th>1</th>
      <td>175</td>
      <td>0.00</td>
      <td>0.0</td>
      <td>158.70</td>
      <td>159.06</td>
      <td>146 x 146</td>
      <td>0</td>
      <td>0</td>
    </tr>
    <tr>
      <th>2</th>
      <td>180</td>
      <td>0.00</td>
      <td>0.0</td>
      <td>153.69</td>
      <td>154.05</td>
      <td>146 x 102</td>
      <td>0</td>
      <td>0</td>
    </tr>
    <tr>
      <th>3</th>
      <td>185</td>
      <td>0.00</td>
      <td>0.0</td>
      <td>148.70</td>
      <td>149.04</td>
      <td>146 x 102</td>
      <td>0</td>
      <td>0</td>
    </tr>
    <tr>
      <th>4</th>
      <td>190</td>
      <td>0.00</td>
      <td>0.0</td>
      <td>143.70</td>
      <td>144.06</td>
      <td>146 x 146</td>
      <td>0</td>
      <td>0</td>
    </tr>
    <tr>
      <th>...</th>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
    </tr>
    <tr>
      <th>110</th>
      <td>405</td>
      <td>0.01</td>
      <td>0.0</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 11065</td>
      <td>0</td>
      <td>28</td>
    </tr>
    <tr>
      <th>111</th>
      <td>410</td>
      <td>0.02</td>
      <td>0.0</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 8620</td>
      <td>0</td>
      <td>200</td>
    </tr>
    <tr>
      <th>112</th>
      <td>415</td>
      <td>0.01</td>
      <td>0.0</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 8630</td>
      <td>0</td>
      <td>5</td>
    </tr>
    <tr>
      <th>113</th>
      <td>420</td>
      <td>0.02</td>
      <td>0.0</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 8641</td>
      <td>0</td>
      <td>55</td>
    </tr>
    <tr>
      <th>114</th>
      <td>425</td>
      <td>0.05</td>
      <td>0.0</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 9471</td>
      <td>0</td>
      <td>20</td>
    </tr>
  </tbody>
</table>
<p>115 rows × 8 columns</p>
</div>




```python
puts = pd.read_excel("SPY_191020exp_290920.xlsx", sheet_name="put")
puts
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>strike</th>
      <th>lastPrice</th>
      <th>change</th>
      <th>bid</th>
      <th>ask</th>
      <th>BxA</th>
      <th>volume</th>
      <th>openInterest</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>170</td>
      <td>0.01</td>
      <td>0.0</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 14192</td>
      <td>0</td>
      <td>12433</td>
    </tr>
    <tr>
      <th>1</th>
      <td>175</td>
      <td>0.02</td>
      <td>0.0</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 12853</td>
      <td>0</td>
      <td>1001</td>
    </tr>
    <tr>
      <th>2</th>
      <td>180</td>
      <td>0.01</td>
      <td>0.0</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 13285</td>
      <td>0</td>
      <td>29</td>
    </tr>
    <tr>
      <th>3</th>
      <td>185</td>
      <td>0.01</td>
      <td>0.0</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 12094</td>
      <td>0</td>
      <td>3</td>
    </tr>
    <tr>
      <th>4</th>
      <td>190</td>
      <td>0.01</td>
      <td>0.0</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 11473</td>
      <td>0</td>
      <td>311</td>
    </tr>
    <tr>
      <th>...</th>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
    </tr>
    <tr>
      <th>110</th>
      <td>405</td>
      <td>0.00</td>
      <td>0.0</td>
      <td>70.96</td>
      <td>71.29</td>
      <td>150 x 150</td>
      <td>0</td>
      <td>0</td>
    </tr>
    <tr>
      <th>111</th>
      <td>410</td>
      <td>0.00</td>
      <td>0.0</td>
      <td>75.97</td>
      <td>76.28</td>
      <td>150 x 150</td>
      <td>0</td>
      <td>0</td>
    </tr>
    <tr>
      <th>112</th>
      <td>415</td>
      <td>0.00</td>
      <td>0.0</td>
      <td>80.98</td>
      <td>81.27</td>
      <td>150 x 150</td>
      <td>0</td>
      <td>0</td>
    </tr>
    <tr>
      <th>113</th>
      <td>420</td>
      <td>0.00</td>
      <td>0.0</td>
      <td>85.97</td>
      <td>86.27</td>
      <td>150 x 150</td>
      <td>0</td>
      <td>0</td>
    </tr>
    <tr>
      <th>114</th>
      <td>425</td>
      <td>0.00</td>
      <td>0.0</td>
      <td>90.97</td>
      <td>91.28</td>
      <td>150 x 150</td>
      <td>0</td>
      <td>0</td>
    </tr>
  </tbody>
</table>
<p>115 rows × 8 columns</p>
</div>




```python
# Find midprices from bid/asks
calls["midprice"] = (calls.bid + calls.ask)/2
calls = calls[calls.midprice > 0]
puts["midprice"] = (puts.bid + puts.ask)/2
puts = puts[puts.midprice > 0]
calls.tail(30)
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>strike</th>
      <th>lastPrice</th>
      <th>change</th>
      <th>bid</th>
      <th>ask</th>
      <th>BxA</th>
      <th>volume</th>
      <th>openInterest</th>
      <th>midprice</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>85</th>
      <td>353</td>
      <td>0.07</td>
      <td>0.00</td>
      <td>0.07</td>
      <td>0.08</td>
      <td>8069 x 12235</td>
      <td>0</td>
      <td>447</td>
      <td>0.075</td>
    </tr>
    <tr>
      <th>86</th>
      <td>354</td>
      <td>0.06</td>
      <td>0.00</td>
      <td>0.06</td>
      <td>0.07</td>
      <td>3377 x 15224</td>
      <td>0</td>
      <td>621</td>
      <td>0.065</td>
    </tr>
    <tr>
      <th>87</th>
      <td>355</td>
      <td>0.06</td>
      <td>0.01</td>
      <td>0.05</td>
      <td>0.06</td>
      <td>693 x 16589</td>
      <td>10</td>
      <td>1685</td>
      <td>0.055</td>
    </tr>
    <tr>
      <th>88</th>
      <td>356</td>
      <td>0.05</td>
      <td>-0.01</td>
      <td>0.04</td>
      <td>0.05</td>
      <td>6205 x 13597</td>
      <td>1</td>
      <td>133</td>
      <td>0.045</td>
    </tr>
    <tr>
      <th>89</th>
      <td>357</td>
      <td>0.05</td>
      <td>0.00</td>
      <td>0.03</td>
      <td>0.04</td>
      <td>12642 x 10876</td>
      <td>0</td>
      <td>217</td>
      <td>0.035</td>
    </tr>
    <tr>
      <th>90</th>
      <td>358</td>
      <td>0.05</td>
      <td>0.00</td>
      <td>0.02</td>
      <td>0.03</td>
      <td>24527 x 3035</td>
      <td>0</td>
      <td>931</td>
      <td>0.025</td>
    </tr>
    <tr>
      <th>91</th>
      <td>359</td>
      <td>0.03</td>
      <td>0.00</td>
      <td>0.02</td>
      <td>0.03</td>
      <td>17058 x 11058</td>
      <td>1</td>
      <td>166</td>
      <td>0.025</td>
    </tr>
    <tr>
      <th>92</th>
      <td>360</td>
      <td>0.02</td>
      <td>-0.01</td>
      <td>0.02</td>
      <td>0.03</td>
      <td>3404 x 14911</td>
      <td>13</td>
      <td>857</td>
      <td>0.025</td>
    </tr>
    <tr>
      <th>93</th>
      <td>361</td>
      <td>0.01</td>
      <td>-0.01</td>
      <td>0.01</td>
      <td>0.02</td>
      <td>24305 x 4582</td>
      <td>6</td>
      <td>867</td>
      <td>0.015</td>
    </tr>
    <tr>
      <th>94</th>
      <td>362</td>
      <td>0.02</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0.02</td>
      <td>18921 x 11277</td>
      <td>0</td>
      <td>156</td>
      <td>0.015</td>
    </tr>
    <tr>
      <th>95</th>
      <td>363</td>
      <td>0.02</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0.02</td>
      <td>15318 x 13311</td>
      <td>0</td>
      <td>37</td>
      <td>0.015</td>
    </tr>
    <tr>
      <th>96</th>
      <td>364</td>
      <td>0.01</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0.02</td>
      <td>6999 x 15467</td>
      <td>0</td>
      <td>123</td>
      <td>0.015</td>
    </tr>
    <tr>
      <th>97</th>
      <td>365</td>
      <td>0.02</td>
      <td>0.01</td>
      <td>0.01</td>
      <td>0.02</td>
      <td>1916 x 16969</td>
      <td>1</td>
      <td>647</td>
      <td>0.015</td>
    </tr>
    <tr>
      <th>98</th>
      <td>366</td>
      <td>0.01</td>
      <td>0.00</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 2241</td>
      <td>0</td>
      <td>40</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>99</th>
      <td>367</td>
      <td>0.01</td>
      <td>0.00</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 3176</td>
      <td>0</td>
      <td>71</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>100</th>
      <td>368</td>
      <td>0.01</td>
      <td>0.00</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 10160</td>
      <td>0</td>
      <td>6</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>101</th>
      <td>369</td>
      <td>0.01</td>
      <td>0.00</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 12132</td>
      <td>0</td>
      <td>21</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>102</th>
      <td>370</td>
      <td>0.01</td>
      <td>0.00</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 12649</td>
      <td>0</td>
      <td>571</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>103</th>
      <td>371</td>
      <td>0.01</td>
      <td>0.00</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 13388</td>
      <td>0</td>
      <td>11</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>104</th>
      <td>375</td>
      <td>0.01</td>
      <td>0.00</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 13630</td>
      <td>0</td>
      <td>424</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>105</th>
      <td>380</td>
      <td>0.01</td>
      <td>0.00</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 8050</td>
      <td>0</td>
      <td>1072</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>106</th>
      <td>385</td>
      <td>0.01</td>
      <td>0.00</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 9711</td>
      <td>0</td>
      <td>638</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>107</th>
      <td>390</td>
      <td>0.01</td>
      <td>0.00</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 8524</td>
      <td>0</td>
      <td>67</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>108</th>
      <td>395</td>
      <td>0.02</td>
      <td>0.00</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 9632</td>
      <td>0</td>
      <td>938</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>109</th>
      <td>400</td>
      <td>0.01</td>
      <td>0.00</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 9658</td>
      <td>0</td>
      <td>306</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>110</th>
      <td>405</td>
      <td>0.01</td>
      <td>0.00</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 11065</td>
      <td>0</td>
      <td>28</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>111</th>
      <td>410</td>
      <td>0.02</td>
      <td>0.00</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 8620</td>
      <td>0</td>
      <td>200</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>112</th>
      <td>415</td>
      <td>0.01</td>
      <td>0.00</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 8630</td>
      <td>0</td>
      <td>5</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>113</th>
      <td>420</td>
      <td>0.02</td>
      <td>0.00</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 8641</td>
      <td>0</td>
      <td>55</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>114</th>
      <td>425</td>
      <td>0.05</td>
      <td>0.00</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 9471</td>
      <td>0</td>
      <td>20</td>
      <td>0.005</td>
    </tr>
  </tbody>
</table>
</div>




```python
calls.iloc[70:83].reset_index(drop=True)
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>strike</th>
      <th>lastPrice</th>
      <th>change</th>
      <th>bid</th>
      <th>ask</th>
      <th>BxA</th>
      <th>volume</th>
      <th>openInterest</th>
      <th>midprice</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>338</td>
      <td>1.99</td>
      <td>-0.15</td>
      <td>2.01</td>
      <td>2.03</td>
      <td>138 x 754</td>
      <td>77</td>
      <td>1718</td>
      <td>2.020</td>
    </tr>
    <tr>
      <th>1</th>
      <td>339</td>
      <td>1.69</td>
      <td>-0.11</td>
      <td>1.66</td>
      <td>1.68</td>
      <td>185 x 550</td>
      <td>121</td>
      <td>970</td>
      <td>1.670</td>
    </tr>
    <tr>
      <th>2</th>
      <td>340</td>
      <td>1.38</td>
      <td>-0.14</td>
      <td>1.36</td>
      <td>1.37</td>
      <td>121 x 1364</td>
      <td>426</td>
      <td>14048</td>
      <td>1.365</td>
    </tr>
    <tr>
      <th>3</th>
      <td>341</td>
      <td>1.08</td>
      <td>-0.21</td>
      <td>1.09</td>
      <td>1.11</td>
      <td>250 x 1089</td>
      <td>37</td>
      <td>2455</td>
      <td>1.100</td>
    </tr>
    <tr>
      <th>4</th>
      <td>342</td>
      <td>0.87</td>
      <td>-0.14</td>
      <td>0.87</td>
      <td>0.88</td>
      <td>250 x 131</td>
      <td>191</td>
      <td>11182</td>
      <td>0.875</td>
    </tr>
    <tr>
      <th>5</th>
      <td>343</td>
      <td>0.74</td>
      <td>-0.07</td>
      <td>0.69</td>
      <td>0.70</td>
      <td>119 x 3308</td>
      <td>103</td>
      <td>10889</td>
      <td>0.695</td>
    </tr>
    <tr>
      <th>6</th>
      <td>344</td>
      <td>0.55</td>
      <td>-0.08</td>
      <td>0.54</td>
      <td>0.55</td>
      <td>171 x 3441</td>
      <td>56</td>
      <td>1694</td>
      <td>0.545</td>
    </tr>
    <tr>
      <th>7</th>
      <td>345</td>
      <td>0.44</td>
      <td>-0.06</td>
      <td>0.42</td>
      <td>0.43</td>
      <td>277 x 2330</td>
      <td>52</td>
      <td>2292</td>
      <td>0.425</td>
    </tr>
    <tr>
      <th>8</th>
      <td>346</td>
      <td>0.34</td>
      <td>-0.07</td>
      <td>0.32</td>
      <td>0.33</td>
      <td>5831 x 310</td>
      <td>69</td>
      <td>721</td>
      <td>0.325</td>
    </tr>
    <tr>
      <th>9</th>
      <td>347</td>
      <td>0.28</td>
      <td>-0.03</td>
      <td>0.25</td>
      <td>0.26</td>
      <td>6338 x 1368</td>
      <td>11</td>
      <td>922</td>
      <td>0.255</td>
    </tr>
    <tr>
      <th>10</th>
      <td>348</td>
      <td>0.21</td>
      <td>-0.01</td>
      <td>0.20</td>
      <td>0.21</td>
      <td>3149 x 5999</td>
      <td>21</td>
      <td>545</td>
      <td>0.205</td>
    </tr>
    <tr>
      <th>11</th>
      <td>349</td>
      <td>0.16</td>
      <td>-0.03</td>
      <td>0.16</td>
      <td>0.17</td>
      <td>315 x 8259</td>
      <td>24</td>
      <td>691</td>
      <td>0.165</td>
    </tr>
    <tr>
      <th>12</th>
      <td>350</td>
      <td>0.14</td>
      <td>0.00</td>
      <td>0.12</td>
      <td>0.13</td>
      <td>13692 x 334</td>
      <td>193</td>
      <td>1840</td>
      <td>0.125</td>
    </tr>
  </tbody>
</table>
</div>




```python
puts.head(10)
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>strike</th>
      <th>lastPrice</th>
      <th>change</th>
      <th>bid</th>
      <th>ask</th>
      <th>BxA</th>
      <th>volume</th>
      <th>openInterest</th>
      <th>midprice</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>170</td>
      <td>0.01</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.01</td>
      <td>0 x 14192</td>
      <td>0</td>
      <td>12433</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>1</th>
      <td>175</td>
      <td>0.02</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.01</td>
      <td>0 x 12853</td>
      <td>0</td>
      <td>1001</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>2</th>
      <td>180</td>
      <td>0.01</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.01</td>
      <td>0 x 13285</td>
      <td>0</td>
      <td>29</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>3</th>
      <td>185</td>
      <td>0.01</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.01</td>
      <td>0 x 12094</td>
      <td>0</td>
      <td>3</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>4</th>
      <td>190</td>
      <td>0.01</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.01</td>
      <td>0 x 11473</td>
      <td>0</td>
      <td>311</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>5</th>
      <td>195</td>
      <td>0.01</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.01</td>
      <td>0 x 11527</td>
      <td>0</td>
      <td>302</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>6</th>
      <td>200</td>
      <td>0.01</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.01</td>
      <td>0 x 10963</td>
      <td>0</td>
      <td>446</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>7</th>
      <td>205</td>
      <td>0.02</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.01</td>
      <td>0 x 7905</td>
      <td>0</td>
      <td>332</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>8</th>
      <td>210</td>
      <td>0.01</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.01</td>
      <td>0 x 8400</td>
      <td>0</td>
      <td>27</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>9</th>
      <td>215</td>
      <td>0.02</td>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.01</td>
      <td>0 x 10678</td>
      <td>0</td>
      <td>103</td>
      <td>0.005</td>
    </tr>
  </tbody>
</table>
</div>




```python
# Visualise put and call prices
fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(12,6))
ax0.scatter(calls.strike, calls.midprice);
ax1.scatter(puts.strike, puts.midprice);
plt.show()
```


```python
# Construct butterflies
data = []

for (_, left) ,(_,centre), (_, right) in zip(calls.iterrows(), calls.iloc[1:].iterrows(), calls.iloc[2:].iterrows()):
    # Filter out all zero volume
    if not any(vol > 0 for vol in {left.volume, centre.volume, right.volume}):
        continue
    # Filter out any zero open interest
    if not all(oi > 0 for oi in {left.openInterest, centre.openInterest, right.openInterest}):
        continue
    # Equidistant on either end
    if centre.strike - left.strike != right.strike - centre.strike:
        continue
    butterfly_price = left.midprice - 2* centre.midprice + right.midprice
    max_profit = centre.strike - left.strike
    data.append([centre.strike, butterfly_price, max_profit])

bflys = pd.DataFrame(data, columns=["strike", "price", "max_profit"])
bflys["prob"] = bflys.price / bflys.max_profit
bflys
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>strike</th>
      <th>price</th>
      <th>max_profit</th>
      <th>prob</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>322</td>
      <td>1.000000e-02</td>
      <td>1</td>
      <td>1.000000e-02</td>
    </tr>
    <tr>
      <th>1</th>
      <td>323</td>
      <td>1.500000e-02</td>
      <td>1</td>
      <td>1.500000e-02</td>
    </tr>
    <tr>
      <th>2</th>
      <td>324</td>
      <td>1.500000e-02</td>
      <td>1</td>
      <td>1.500000e-02</td>
    </tr>
    <tr>
      <th>3</th>
      <td>325</td>
      <td>3.000000e-02</td>
      <td>1</td>
      <td>3.000000e-02</td>
    </tr>
    <tr>
      <th>4</th>
      <td>326</td>
      <td>1.500000e-02</td>
      <td>1</td>
      <td>1.500000e-02</td>
    </tr>
    <tr>
      <th>5</th>
      <td>327</td>
      <td>2.500000e-02</td>
      <td>1</td>
      <td>2.500000e-02</td>
    </tr>
    <tr>
      <th>6</th>
      <td>328</td>
      <td>2.500000e-02</td>
      <td>1</td>
      <td>2.500000e-02</td>
    </tr>
    <tr>
      <th>7</th>
      <td>329</td>
      <td>3.000000e-02</td>
      <td>1</td>
      <td>3.000000e-02</td>
    </tr>
    <tr>
      <th>8</th>
      <td>330</td>
      <td>3.500000e-02</td>
      <td>1</td>
      <td>3.500000e-02</td>
    </tr>
    <tr>
      <th>9</th>
      <td>331</td>
      <td>3.000000e-02</td>
      <td>1</td>
      <td>3.000000e-02</td>
    </tr>
    <tr>
      <th>10</th>
      <td>332</td>
      <td>4.000000e-02</td>
      <td>1</td>
      <td>4.000000e-02</td>
    </tr>
    <tr>
      <th>11</th>
      <td>333</td>
      <td>4.500000e-02</td>
      <td>1</td>
      <td>4.500000e-02</td>
    </tr>
    <tr>
      <th>12</th>
      <td>334</td>
      <td>3.500000e-02</td>
      <td>1</td>
      <td>3.500000e-02</td>
    </tr>
    <tr>
      <th>13</th>
      <td>335</td>
      <td>4.500000e-02</td>
      <td>1</td>
      <td>4.500000e-02</td>
    </tr>
    <tr>
      <th>14</th>
      <td>336</td>
      <td>4.000000e-02</td>
      <td>1</td>
      <td>4.000000e-02</td>
    </tr>
    <tr>
      <th>15</th>
      <td>337</td>
      <td>5.000000e-02</td>
      <td>1</td>
      <td>5.000000e-02</td>
    </tr>
    <tr>
      <th>16</th>
      <td>338</td>
      <td>4.500000e-02</td>
      <td>1</td>
      <td>4.500000e-02</td>
    </tr>
    <tr>
      <th>17</th>
      <td>339</td>
      <td>4.500000e-02</td>
      <td>1</td>
      <td>4.500000e-02</td>
    </tr>
    <tr>
      <th>18</th>
      <td>340</td>
      <td>4.000000e-02</td>
      <td>1</td>
      <td>4.000000e-02</td>
    </tr>
    <tr>
      <th>19</th>
      <td>341</td>
      <td>4.000000e-02</td>
      <td>1</td>
      <td>4.000000e-02</td>
    </tr>
    <tr>
      <th>20</th>
      <td>342</td>
      <td>4.500000e-02</td>
      <td>1</td>
      <td>4.500000e-02</td>
    </tr>
    <tr>
      <th>21</th>
      <td>343</td>
      <td>3.000000e-02</td>
      <td>1</td>
      <td>3.000000e-02</td>
    </tr>
    <tr>
      <th>22</th>
      <td>344</td>
      <td>3.000000e-02</td>
      <td>1</td>
      <td>3.000000e-02</td>
    </tr>
    <tr>
      <th>23</th>
      <td>345</td>
      <td>2.000000e-02</td>
      <td>1</td>
      <td>2.000000e-02</td>
    </tr>
    <tr>
      <th>24</th>
      <td>346</td>
      <td>3.000000e-02</td>
      <td>1</td>
      <td>3.000000e-02</td>
    </tr>
    <tr>
      <th>25</th>
      <td>347</td>
      <td>2.000000e-02</td>
      <td>1</td>
      <td>2.000000e-02</td>
    </tr>
    <tr>
      <th>26</th>
      <td>348</td>
      <td>1.000000e-02</td>
      <td>1</td>
      <td>1.000000e-02</td>
    </tr>
    <tr>
      <th>27</th>
      <td>349</td>
      <td>0.000000e+00</td>
      <td>1</td>
      <td>0.000000e+00</td>
    </tr>
    <tr>
      <th>28</th>
      <td>350</td>
      <td>2.000000e-02</td>
      <td>1</td>
      <td>2.000000e-02</td>
    </tr>
    <tr>
      <th>29</th>
      <td>351</td>
      <td>-2.775558e-17</td>
      <td>1</td>
      <td>-2.775558e-17</td>
    </tr>
    <tr>
      <th>30</th>
      <td>352</td>
      <td>1.000000e-02</td>
      <td>1</td>
      <td>1.000000e-02</td>
    </tr>
    <tr>
      <th>31</th>
      <td>353</td>
      <td>-2.775558e-17</td>
      <td>1</td>
      <td>-2.775558e-17</td>
    </tr>
    <tr>
      <th>32</th>
      <td>354</td>
      <td>6.938894e-18</td>
      <td>1</td>
      <td>6.938894e-18</td>
    </tr>
    <tr>
      <th>33</th>
      <td>355</td>
      <td>0.000000e+00</td>
      <td>1</td>
      <td>0.000000e+00</td>
    </tr>
    <tr>
      <th>34</th>
      <td>356</td>
      <td>6.938894e-18</td>
      <td>1</td>
      <td>6.938894e-18</td>
    </tr>
    <tr>
      <th>35</th>
      <td>357</td>
      <td>-6.938894e-18</td>
      <td>1</td>
      <td>-6.938894e-18</td>
    </tr>
    <tr>
      <th>36</th>
      <td>358</td>
      <td>1.000000e-02</td>
      <td>1</td>
      <td>1.000000e-02</td>
    </tr>
    <tr>
      <th>37</th>
      <td>359</td>
      <td>0.000000e+00</td>
      <td>1</td>
      <td>0.000000e+00</td>
    </tr>
    <tr>
      <th>38</th>
      <td>360</td>
      <td>-1.000000e-02</td>
      <td>1</td>
      <td>-1.000000e-02</td>
    </tr>
    <tr>
      <th>39</th>
      <td>361</td>
      <td>1.000000e-02</td>
      <td>1</td>
      <td>1.000000e-02</td>
    </tr>
    <tr>
      <th>40</th>
      <td>362</td>
      <td>0.000000e+00</td>
      <td>1</td>
      <td>0.000000e+00</td>
    </tr>
    <tr>
      <th>41</th>
      <td>364</td>
      <td>0.000000e+00</td>
      <td>1</td>
      <td>0.000000e+00</td>
    </tr>
    <tr>
      <th>42</th>
      <td>365</td>
      <td>-1.000000e-02</td>
      <td>1</td>
      <td>-1.000000e-02</td>
    </tr>
    <tr>
      <th>43</th>
      <td>366</td>
      <td>1.000000e-02</td>
      <td>1</td>
      <td>1.000000e-02</td>
    </tr>
  </tbody>
</table>
</div>




```python
# SPY was trading around 335 when this data was collected
plt.rcParams.update({'font.size': 16})
plt.figure(figsize=(9,6))
plt.scatter(bflys.strike, bflys.prob);
plt.xlabel("Strike")
plt.ylabel("Probability")
plt.show()
# plt.savefig("SPY_raw_bfly_prob.png", dpi=300)
```


```python
from scipy.ndimage import gaussian_filter1d

smoothed_prob = gaussian_filter1d(bflys.prob, 2)

plt.figure(figsize=(9,6))
plt.plot(bflys.strike, bflys.prob, "o", bflys.strike, smoothed_prob, "rx")
plt.legend(["raw prob", "smoothed prob"], loc="best")
plt.xlabel("Strike")
plt.ylabel("Probability")
plt.show()
# plt.savefig("SPY_smooth_bfly_prob.png", dpi=300)
```


```python
plt.figure(figsize=(9,6))
pdf = scipy.interpolate.interp1d(bflys.strike, smoothed_prob, kind="cubic",
                                 fill_value="extrapolate")
x_new = np.linspace(bflys.strike.min(), bflys.strike.max(), 100)
plt.plot(bflys.strike, smoothed_prob, "rx", x_new, pdf(x_new), "k-");
plt.legend(["smoothed prob", "fitted PDF"], loc="best")
plt.xlabel("K")
plt.ylabel("f(K)")
plt.tight_layout()
plt.show()
# plt.savefig("SPY_bfly_pdf.png", dpi=300)
```


```python
# Find area under curve
raw_total_prob = scipy.integrate.trapz(smoothed_prob, bflys.strike)
print(f"Raw total probability: {raw_total_prob}")

# Don't need to normalise because there is mass in the left tail that we are ignoring

# # Normalise
# normalised_prob = smoothed_prob / raw_total_prob
# total_prob = scipy.integrate.trapz(normalised_prob, bflys.strike)
# print(f"Normalised total probability: {total_prob}")
# # should be less than 1
```

    Raw total probability: 0.872620109383161


### Put butterflies

We can repeat the process with put butterflies instead of call butterflies.


```python
# Repeating the same with put butterflies
from scipy.ndimage import gaussian_filter1d

data = []

for (_, left) ,(_,centre), (_, right) in zip(puts.iterrows(), puts.iloc[1:].iterrows(), puts.iloc[2:].iterrows()):
    # Filter out all zero volume
    if not any(vol > 0 for vol in {left.volume, centre.volume, right.volume}):
        continue
    # Filter out any zero open interest
    if not all(oi > 0 for oi in {left.openInterest, centre.openInterest, right.openInterest}):
        continue
    # Equidistant on either end
    if centre.strike - left.strike != right.strike - centre.strike:
        continue
    butterfly_price = left.midprice - 2* centre.midprice + right.midprice
    max_profit = centre.strike - left.strike
    data.append([centre.strike, butterfly_price, max_profit])

put_bflys = pd.DataFrame(data, columns=["strike", "price", "max_profit"])
put_bflys["prob"] = put_bflys.price / put_bflys.max_profit
smoothed_prob_put = gaussian_filter1d(put_bflys.prob, 2)
plt.plot(put_bflys.strike, put_bflys.prob, "o", put_bflys.strike, smoothed_prob_put, "rx")
plt.legend(["raw prob", "smoothed prob"], loc="best")
plt.show()
```


```python
put_pdf = scipy.interpolate.interp1d(put_bflys.strike, smoothed_prob_put, kind="cubic", 
                                     fill_value=0.0)
x_new = np.linspace(bflys.strike.min(), put_bflys.strike.max(), 100)
plt.plot(x_new, pdf(x_new), "m-", x_new, put_pdf(x_new), "k-");
plt.legend(["call PDF", "put PDF"], loc="best")
plt.show()
```

### Refactor butterfly PDF code


```python
def construct_pdf(calls_df, make_plot=True, fill_value="extrapolate"):
    if "midprice" not in calls_df.columns:
        calls_df["midprice"] = (calls_df.bid + calls_df.ask) /2
    
    # Construct butterflies
    data = []

    for (_, left) ,(_,centre), (_, right) in zip(calls_df.iterrows(), calls_df.iloc[1:].iterrows(), calls_df.iloc[2:].iterrows()):
        # Filter out all zero volume
        if not any(vol > 0 for vol in {left.volume, centre.volume, right.volume}):
            continue
        # Filter out any zero open interest
        if not all(oi > 0 for oi in {left.openInterest, centre.openInterest, right.openInterest}):
            continue
        # Equidistant on either end
        if centre.strike - left.strike != right.strike - centre.strike:
            continue
        butterfly_price = left.midprice - 2* centre.midprice + right.midprice
        max_profit = centre.strike - left.strike
        data.append([centre.strike, butterfly_price, max_profit])

    bflys = pd.DataFrame(data, columns=["strike", "price", "max_profit"])
    bflys["prob"] = bflys.price / bflys.max_profit

    smoothed_prob = gaussian_filter1d(bflys.prob, 2)
    pdf = scipy.interpolate.interp1d(bflys.strike, smoothed_prob, kind="cubic",
                                     fill_value=fill_value)
    if not make_plot:
        return pdf
    
    plt.figure(figsize=(9,6))
    x_new = np.linspace(bflys.strike.min(), bflys.strike.max(), 100)
    plt.plot(bflys.strike, smoothed_prob, "rx", x_new, pdf(x_new), "k-");
    plt.legend(["smoothed prob", "fitted PDF"], loc="best")
    plt.xlabel("K")
    plt.ylabel("f(K)")


```

## Implied PDF from Breeden-Litzenberger

$$ \frac{\partial^2 C}{\partial K^2} = e^{r\tau} \frac{\partial^2 C}{\partial K^2} $$

So we need to interpolate. The choice is between price space and volatility space.


```python
calls = pd.read_excel("SPY_191020exp_290920.xlsx", sheet_name="call")
calls["midprice"] = (calls.bid + calls.ask)/2
calls = calls[calls.midprice > 0]
calls
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>strike</th>
      <th>lastPrice</th>
      <th>change</th>
      <th>bid</th>
      <th>ask</th>
      <th>BxA</th>
      <th>volume</th>
      <th>openInterest</th>
      <th>midprice</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>170</td>
      <td>0.00</td>
      <td>0.0</td>
      <td>163.70</td>
      <td>164.06</td>
      <td>146 x 146</td>
      <td>0</td>
      <td>0</td>
      <td>163.880</td>
    </tr>
    <tr>
      <th>1</th>
      <td>175</td>
      <td>0.00</td>
      <td>0.0</td>
      <td>158.70</td>
      <td>159.06</td>
      <td>146 x 146</td>
      <td>0</td>
      <td>0</td>
      <td>158.880</td>
    </tr>
    <tr>
      <th>2</th>
      <td>180</td>
      <td>0.00</td>
      <td>0.0</td>
      <td>153.69</td>
      <td>154.05</td>
      <td>146 x 102</td>
      <td>0</td>
      <td>0</td>
      <td>153.870</td>
    </tr>
    <tr>
      <th>3</th>
      <td>185</td>
      <td>0.00</td>
      <td>0.0</td>
      <td>148.70</td>
      <td>149.04</td>
      <td>146 x 102</td>
      <td>0</td>
      <td>0</td>
      <td>148.870</td>
    </tr>
    <tr>
      <th>4</th>
      <td>190</td>
      <td>0.00</td>
      <td>0.0</td>
      <td>143.70</td>
      <td>144.06</td>
      <td>146 x 146</td>
      <td>0</td>
      <td>0</td>
      <td>143.880</td>
    </tr>
    <tr>
      <th>...</th>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
    </tr>
    <tr>
      <th>110</th>
      <td>405</td>
      <td>0.01</td>
      <td>0.0</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 11065</td>
      <td>0</td>
      <td>28</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>111</th>
      <td>410</td>
      <td>0.02</td>
      <td>0.0</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 8620</td>
      <td>0</td>
      <td>200</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>112</th>
      <td>415</td>
      <td>0.01</td>
      <td>0.0</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 8630</td>
      <td>0</td>
      <td>5</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>113</th>
      <td>420</td>
      <td>0.02</td>
      <td>0.0</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 8641</td>
      <td>0</td>
      <td>55</td>
      <td>0.005</td>
    </tr>
    <tr>
      <th>114</th>
      <td>425</td>
      <td>0.05</td>
      <td>0.0</td>
      <td>0.00</td>
      <td>0.01</td>
      <td>0 x 9471</td>
      <td>0</td>
      <td>20</td>
      <td>0.005</td>
    </tr>
  </tbody>
</table>
<p>115 rows × 9 columns</p>
</div>




```python
calls_sub = calls[(calls.strike > 250) & (calls.strike < 400)]
plt.figure(figsize=(12,6))
plt.plot(calls_sub.strike, calls_sub.midprice, ".");
plt.xlabel("strike")
plt.ylabel("price")
plt.savefig("call_prices.png", dpi=400)
plt.show()
```

### B-S IV

B-S doesn't have a closed-form inverse, which means to find the IV from the price we must resort to numerical methods. This in itself can be quite a deep topic, but we are going to be simplistic in our approach and just use Newton's method. To do, we need a function that computes the call price from inputs (already written above), and a function to evaluate the gradient of the call price function with respect to the volatility parameter - this is known as **vega**. 


```python
def call_value(S, K, sigma, t=0, r=0):
    # use np.multiply and divide to handle divide-by-zero
    with np.errstate(divide='ignore'):
        d1 = np.divide(1, sigma * np.sqrt(t)) * (np.log(S/K) + (r+sigma**2 / 2) * t)
        d2 = d1 - sigma * np.sqrt(t)
    return np.multiply(norm.cdf(d1),S) - np.multiply(norm.cdf(d2), K * np.exp(-r * t))

def call_vega(S, K, sigma, t=0, r=0):
    with np.errstate(divide='ignore'):
        d1 = np.divide(1, sigma * np.sqrt(t)) * (np.log(S/K) + (r+sigma**2 / 2) * t)
    return np.multiply(S, norm.pdf(d1)) * np.sqrt(t)

def bs_iv(price, S, K, t=0, r=0, precision=1e-4, initial_guess=0.2, max_iter=1000, verbose=False):
    iv = initial_guess
    for _ in range(max_iter):
        P = call_value(S, K, iv, t, r)
        diff = price - P
        if abs(diff) < precision:
            return iv
        grad = call_vega(S, K, iv, t, r)
        iv += diff/grad     
    if verbose:
        print(f"Did not converge after {max_iter} iterations")
    return iv 
```


```python
c_test = call_value(100, 110, 0.2, t=1)
print(c_test)
```

    4.292010941409885



```python
# Check that it works
bs_iv(c_test, 100,110, t=1)
```




    0.2




```python
S = 332
t = 3/52
calls["iv"] = calls.apply(lambda row: bs_iv(row.midprice, S, row.strike, t, max_iter=500), axis=1)
```

    <ipython-input-21-08b1dc83cea8>:21: RuntimeWarning: divide by zero encountered in double_scalars
      iv += diff/grad
    <ipython-input-21-08b1dc83cea8>:4: RuntimeWarning: invalid value encountered in double_scalars
      d1 = np.divide(1, sigma * np.sqrt(t)) * (np.log(S/K) + (r+sigma**2 / 2) * t)
    <ipython-input-21-08b1dc83cea8>:10: RuntimeWarning: invalid value encountered in double_scalars
      d1 = np.divide(1, sigma * np.sqrt(t)) * (np.log(S/K) + (r+sigma**2 / 2) * t)



```python
def plot_vol_smile(calls, savefig=False):
    plt.figure(figsize=(9,6))
    plt.plot(calls.strike, calls.iv, ".")
    plt.xlabel("Strike")
    plt.ylabel("IV")
    if savefig:
        plt.savefig("vol_smile.png",dpi=300)
    plt.show()
```


```python
calls_no_na = calls.dropna()

fig, (ax1, ax2) = plt.subplots(2,1, figsize=(6,7), sharex=True)
ax1.plot(calls_no_na.strike, calls_no_na.midprice, "r.")
ax1.set_ylabel("Call price")
ax2.plot(calls_no_na.strike, calls_no_na.iv, ".")
ax2.set_ylabel("IV")
ax2.set_xlabel("Strike")
plt.tight_layout()
# plt.savefig("calls_to_iv.png", dpi=400)
plt.show()
```


```python
calls_clean = calls.dropna().copy()
calls_clean["iv"] = gaussian_filter1d(calls_clean.iv, 3)
```


```python
plot_vol_smile(calls_clean)
```


```python
calls_clean = calls_clean[(calls_clean.strike > 300) & (calls_clean.strike < 375)]
```


```python
plot_vol_smile(calls_clean, savefig=False)
```


```python
plt.figure(figsize=(9,6))
vol_surface = scipy.interpolate.interp1d(calls_clean.strike, calls_clean.iv, kind="cubic",
                                 fill_value="extrapolate")
x_new = np.arange(calls_clean.strike.min(), calls_clean.strike.max(), 0.1)
plt.plot(calls_clean.strike, calls_clean.iv, "bx", x_new, vol_surface(x_new), "k-");
plt.legend(["smoothed IV", "fitted smile"], loc="best")
plt.xlabel("Strike")
plt.ylabel("IV")
plt.tight_layout()
# plt.savefig("SPY_smile.png", dpi=300)
plt.show()
```


```python
C_interp = call_value(S, x_new, vol_surface(x_new), t)
```


```python
first_deriv = np.gradient(C_interp, x_new)
second_deriv = np.gradient(first_deriv, x_new)
```


```python
def pdf2(Krange, S, t=0, r=0):
    # x is a range of strikes
    Crange = call_value(S, Krange, vol_surface(Krange), t, r)
    first_deriv = np.gradient(Crange, x_new, edge_order=0)
    second_deriv = np.gradient(first_deriv, x_new, edge_order=0)
    return np.exp(r * t) * second_deriv
```


```python
fig, ax1 = plt.subplots(figsize=(9,6))

col="slategrey"
ax1.set_xlabel('Strike')
ax1.set_ylabel('Call price', color=col)
ax1.plot(x_new, C_interp, color=col)
ax1.tick_params(axis='y', labelcolor=col)
ax1.axvline(332, color="k", linestyle="--")

ax2 = ax1.twinx()
col = 'darkred'
ax2.set_ylabel('f(K)', color=col)
ax2.plot(x_new, pdf2(x_new, S, t), color=col)
ax2.tick_params(axis='y', labelcolor=col)

fig.tight_layout()
# plt.savefig("volspace_pdf.png", dpi=300)
plt.show()
```


```python
data = np.vstack([x_new, pdf(x_new), pdf2(x_new, S, t)]).T
data = data[(data[:, 0] < 360) & (data[:, 0] > 320)]
k, bfly_pdf, curve_pdf = np.hsplit(data,3)
```


```python
plt.figure(figsize=(9,6))
plt.plot(k, bfly_pdf, "steelblue", k, curve_pdf, "darkred")
plt.xlabel("Strike")
plt.ylabel("f(K)")
plt.legend(["Butterfly PDF", "Breeden-Litzenberger PDF"], loc="best")
# plt.savefig("butterfly_vs_breeden.png", dpi=300)
plt.show()
```

Note that the Breeden-Liztenberger PDF is much more well-behaved OTM (doesn't have the same kinks between 355-360).

Also, we have a lot more data (butterflies each require 3 prices), which is why we can extend the strikes further out as per the previous diagram – though note that even Breeden-Litzenberger starts to behave nonsensically.

BUt we have one other thing to compare it to: Breeden-Liztenberger applied directly to the interpolated calls (without interpolating in IV space).


```python
calls["price_clean"] = gaussian_filter1d(calls.midprice, 3)
price_interp = scipy.interpolate.interp1d(calls.strike, calls.price_clean, kind="cubic",
                                 fill_value="extrapolate")
C_interp_direct = price_interp(x_new)
```


```python
plt.plot(x_new, C_interp - C_interp_direct)
plt.xlabel("Strike")
plt.ylabel("Vol space - call space")
plt.show()
```


```python
def pdf3(calls, Krange, S, t=0, r=0):
    # x is a range of strikes
    calls["price_clean"] = gaussian_filter1d(calls.midprice, 3)
    price_interp = scipy.interpolate.interp1d(calls.strike, calls.price_clean, kind="cubic",
                                     fill_value="extrapolate")
    Crange = price_interp(Krange)
    first_deriv = np.gradient(Crange, x_new, edge_order=0)
    second_deriv = np.gradient(first_deriv, x_new, edge_order=0)
    return second_deriv
```


```python
plt.figure(figsize=(9,6))
plt.plot(x_new, pdf3(calls, x_new, S, t), "darkgreen", x_new, pdf2(x_new, S, t), "darkred")
plt.xlabel("Strike")
plt.ylabel("f(K)")
plt.legend(["Price-space PDF", "IV-space PDF"], loc="best")
# plt.savefig("price_vol_comparison.png", dpi=300)
plt.show()
```


```python
plt.plot(x_new, pdf2(x_new, S, t) - pdf3(calls, x_new, S, t), "r--");
plt.xlabel("Strike")
plt.ylabel("Vol-space PDF - price-space PDF")
plt.show()
```


```python

```

87ca584b08297bbde40a44680e4f27838138a7a6
ipynb
Jupyter Notebook
ProbabilisticValuation/OptionImpliedPDF.ipynb
vishalbelsare/pValuation
5ec4c5a717e007abc641a4b13066dbfcdb4ea38c
BSD-3-Clause
2020-04-27T06:24:34.000Z
2022-02-28T15:55:16.000Z
ProbabilisticValuation/OptionImpliedPDF.ipynb
vishalbelsare/pValuation
5ec4c5a717e007abc641a4b13066dbfcdb4ea38c
BSD-3-Clause
ProbabilisticValuation/OptionImpliedPDF.ipynb
vishalbelsare/pValuation
5ec4c5a717e007abc641a4b13066dbfcdb4ea38c
BSD-3-Clause
2020-05-03T03:13:43.000Z
2022-03-11T14:09:20.000Z
Qwen/Qwen-72B
1. YES
2. YES