## Welcome to The QuantConnect Research PageÂ¶

#### Contribute to this template file https://github.com/QuantConnect/Lean/blob/master/Jupyter/BasicQuantBookTemplate.ipynbÂ¶

In [ ]:
import pandas as pd
from numpy import sqrt,mean,log,diff
import quandl
import scipy.stats as stats
import scipy
from scipy import interpolate
from scipy.interpolate import Rbf
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm, animation
import datetime
%pylab inline

In [2]:
def bsm_price(option_type, sigma, s, k, r, T, q):
# calculate the bsm price of European call and put options
sigma = float(sigma)
d1 = (np.log(s / k) + (r - q + sigma ** 2 * 0.5) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
if option_type == 'c':
price = np.exp(-r*T) * (s * np.exp((r - q)*T) * stats.norm.cdf(d1) - k *  stats.norm.cdf(d2))
return price
elif option_type == 'p':
price = np.exp(-r*T) * (k * stats.norm.cdf(-d2) - s * np.exp((r - q)*T) *  stats.norm.cdf(-d1))
return price
else:
print('No such option type %s') %option_type

In [3]:
def implied_vol(option_type, option_price, s, k, r, T, q):
# apply bisection method to get the implied volatility by solving the BSM function
precision = 0.00001
upper_vol = 500.0
max_vol = 500.0
min_vol = 0.0001
lower_vol = 0.0001
iteration = 0

while 1:
iteration +=1
mid_vol = (upper_vol + lower_vol)/2.0
price = bsm_price(option_type, mid_vol, s, k, r, T, q)
if option_type == 'c':

lower_price = bsm_price(option_type, lower_vol, s, k, r, T, q)
if (lower_price - option_price) * (price - option_price) > 0:
lower_vol = mid_vol
else:
upper_vol = mid_vol
if abs(price - option_price) < precision: break
if mid_vol > max_vol - 5 :
mid_vol = 0.000001
break
#             print("mid_vol=%f" %mid_vol)
#             print("upper_price=%f" %lower_price)

elif option_type == 'p':
upper_price = bsm_price(option_type, upper_vol, s, k, r, T, q)

if (upper_price - option_price) * (price - option_price) > 0:
upper_vol = mid_vol
else:
lower_vol = mid_vol
#             print("mid_vol=%f" %mid_vol)
#             print("upper_price=%f" %upper_price)
if abs(price - option_price) < precision: break
if iteration > 50: break

return mid_vol

In [4]:
# download option data for all expiry months from Yahoo Finance
# provide a formatted DataFrame with a hierarchical index
opt = Options('spy', 'yahoo')


NameErrorTraceback (most recent call last)
<ipython-input-4-1157c9d3cc8b> in <module>()
1 # download option data for all expiry months from Yahoo Finance
2 # provide a formatted DataFrame with a hierarchical index
----> 3 opt = Options('spy', 'yahoo')

NameError: name 'Options' is not defined
In [ ]:
opt.expiry_dates

In [ ]:
# plot the implied volatility of put/call options for specific expiration date
def IV_plot(opt,option_type,expiry_index):
expiry = opt.expiry_dates[expiry_index]
if option_type == 'c':
data = opt.get_call_data(expiry=expiry)
elif option_type == 'p':
data = opt.get_put_data(expiry=expiry)
r = 0.01 # risk free rate
d = 0.01 # continuous devidend yield
s = opt.underlying_price # data_call['Underlying_Price']  undelying price
expiry = data.index.get_level_values('Expiry')[0] # get the expiry
current_date = opt.quote_time#current_date = datetime.datetime.now() # get the current date
time_to_expire = float((expiry - current_date).days)/365 # compute time to expiration
strike = list(data.index.get_level_values('Strike')) # get the strike price
IV = []
for i in range(len(data)):
IV.append(implied_vol(option_type, premium.values[i], s, strike[i], r, time_to_expire, d))

plt.figure(figsize=(16, 7))
a = plt.scatter(strike,IV, c='r', label="IV by solving BSM")
b = plt.scatter(strike,data['IV'],c = 'b', label="IV from Yahoo Finance")
plt.grid()
plt.xlabel('strike')
if option_type == 'c':
plt.ylabel('Implied Volatility for call option')
plt.legend((a,b), ("IV(call) by solving BSM", "IV(call) from Yahoo Finance"))
elif option_type == 'p':
plt.ylabel('Implied Volatility for put options')
plt.legend((a,b), ("IV(put) by solving BSM", "IV(put) from Yahoo Finance"))

return strike,IV

In [ ]:
k_put, IV_put = IV_plot(opt,'p',7)
k_call, IV_call = IV_plot(opt,'c',7)

In [ ]:
plt.figure(figsize=(16, 7))
e = plt.scatter(k_call,IV_call, c ='red', label="IV(call options)")
f = plt.scatter(k_put,IV_put, c = 'black', label="IV(put options)")
plt.xlabel('strike')
plt.ylabel('Implied Volatility')
plt.legend((e,f), ("IV (call options)", "IV (put options)"))

In [ ]:
# generate the volatility scatter plot for all strikes and expiration (for call options)
r = 0.01 # risk free rate
d = 0.01 # continuous devidend yield
expiry_dates = [i for i in opt.expiry_dates if i > opt.quote_time.date()]
current_date = opt.quote_time.date()  ## get the current date
s = opt.underlying_price # undelying price
num_expiry = len(expiry_dates)
IV = [] # (num_expiry * 3)-dimension list with each row being [time_to_expire, strike, implied_vol]

for expiry_index in range(num_expiry-2):
# use call options as an example, for call options use  data = opt.get_call_data(expiry=expiry_dates[expiry_index])
data = opt.get_call_data(expiry=expiry_dates[expiry_index])
expiry = expiry_dates[expiry_index] # get the expiry
time_to_expire = float((expiry - current_date).days)/365 # compute time to expiration
strike = data.index.get_level_values('Strike') # get the strike price
num_strike = len(data)
for j in range(num_strike):
IV.append([time_to_expire, strike[j], implied_vol('c', premium.values[j], s, strike[j], r, time_to_expire, d)])

In [ ]:
# remove volatility with 0 values or greater than 3
IV = [i for i in IV if i[2]>0.01 and i[2]<3]

In [ ]:
x= [IV[i][0] for i in range(len(IV))]
y= [IV[i][1] for i in range(len(IV))]
z= [IV[i][2] for i in range(len(IV))]

In [ ]:
fig = plt.figure(figsize=(20,11))
ax.view_init(30,40)
ax.scatter(x,y,z)

In [ ]:
y_norm = (y-min(y))/(max(y)-min(y)) # normalize the strike price to [0,1]

In [ ]:
points = np.array([[i,j] for i,j in zip(x,y_norm)])
values = np.array(z)
grid_x, grid_y = np.mgrid[min(x):max(x):50j, min(y_norm):max(y_norm):50j]
znew = scipy.interpolate.griddata(points, values, (grid_x, grid_y), method='nearest')
xnew = grid_x
ynew = grid_y

In [ ]:
fig = plt.figure(figsize=(20,15))
ax.view_init(12,15)
ax.plot_wireframe(xnew, ynew, znew, rstride=3, cstride=3)
ax.plot_surface(xnew, ynew, znew, cmap=cm.jet, linewidth=0.001, rstride=1, cstride=1, alpha = 0.70)
ax.set_xlabel('Time to Expiration')
ax.set_ylabel('Strike')
ax.set_zlabel('Implied Volatility')
m = cm.ScalarMappable(cmap=cm.jet)
m.set_array(znew)
cbar = plt.colorbar(m)

In [ ]:
rbf_func = Rbf(xnew,ynew,znew, function='linear', smooth=2)
xnew = np.linspace(min(x),max(x),100)
ynew = np.linspace(min(y_norm),max(y_norm),100)
xnew, ynew = np.meshgrid(xnew, ynew)
znew = rbf_func(xnew, ynew)

In [ ]:
fig = plt.figure(figsize=(20,15))