QuantConnect Logo

Welcome to The QuantConnect Research Page¶

Refer to this page for documentation https://www.quantconnect.com/docs#Introduction-to-Jupyter¶

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
from pandas_datareader.data import Options
import pandas_datareader.data as web
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
    premium = (data['Ask'] + data['Bid'])/2 # option premium
    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    
    premium = (data['Ask'] + data['Bid'])/2.0 # option premium
    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 = fig.add_subplot(111, projection='3d')
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 = fig.add_subplot(111, projection='3d')
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))
ax = fig.add_subplot(111, projection='3d')
ax.view_init(12,15)
ynew = ynew*(max(y)-min(y))+min(y)
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)