import pandas as pd 
import numpy as np 
import matplotlib.pyplot as plt 
%matplotlib inline

raw = pd.read_csv('5Firms.csv.zip',parse_dates=[2,3])
raw.drop(columns=['secid','issuer','exercise_style','last_date'],inplace=True)
raw.head(5)

raw.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 4467804 entries, 14 to 5719675
Data columns (total 14 columns):
date             datetime64[ns]
exdate           datetime64[ns]
cp_flag          object
strike_price     float64
best_bid         float64
best_offer       float64
volume           int64
open_interest    int64
imvol            float64
delta            float64
optionid         int64
ticker           object
index_flag       int64
midquo           float64
dtypes: datetime64[ns](2), float64(6), int64(4), object(2)
memory usage: 511.3+ MB

# change the implied volatility name (it's too long)
raw.rename(columns={'impl_volatility': 'imvol'}, inplace=True)
# transform the strike price and calculate the mid quotes
raw['strike_price'] = raw['strike_price'] / 1000.0
raw['midquo'] = (raw['best_bid'] + raw['best_offer']) / 2
raw.head(5)

# drop contracts with NaN implied volatility or Nan delta
raw.dropna(subset=['imvol'], how='any', inplace=True)
raw.dropna(subset=['delta'], how='any', inplace=True)

# Exclude implied volatility above 150%
raw = raw[raw['imvol'] < 1.5]
# Exclude price below 0.05
raw = raw[raw['best_offer'] > 0.05]
raw.head(5)

# calculate maturity
raw['maturity'] = ((raw['exdate'] - raw['date']) / np.timedelta64(1, 'D')).astype(int)
# Exclude contracts with maturity longer than 240 days and shorter than 7 days
raw = raw[(raw['maturity'] <= 360) & (raw['maturity'] >= 7)]
raw.head(5)

# save the cleaned raw data, because it is time-consuming to processing, next time we just need to load it
hdf = pd.HDFStore('hdf.h5')
hdf['raw'] = raw
hdf.close()

# if next time we start from here
# hdf = pd.HDFStore('hdf.h5')
# raw = hdf['raw']
# hdf.close()

# select variable we want 
clean = raw[['date', 'ticker', 'imvol', 'delta', 'maturity', 'cp_flag']].reset_index(drop=True)
clean.head(5)

print('stocks ticker:\n',clean['ticker'].unique())
print('maturity range:\n',clean['maturity'].min(),clean['maturity'].max())
print('delta range:\n',clean['delta'].min(),clean['delta'].max())

stocks ticker:
 ['COF' 'DD' 'XOM' 'INTC' 'SPX' 'UNH']
maturity range:
 7 360
delta range:
 -0.999934 0.999999

stats = clean.groupby(['ticker','cp_flag'])[['imvol','delta']].describe()
stats

# constant
clean['const'] = 1
# scaled moneyness
clean['delta_n'] = clean['delta']
clean.loc[(clean['cp_flag']=='P'),'delta_n'] = 1 + clean.loc[(clean['cp_flag']=='P'),'delta']
clean['delta_n'] = clean['delta_n'] - 0.5
# scaled maturity
clean['matur_n'] = clean['maturity']/360 - 0.5
# quadratic term of moneyness
clean['delta2'] = clean['delta_n']**2
# interaction between maturity and moneyness
clean['inter'] = clean['delta_n']*clean['matur_n']
clean.head(5)

import statsmodels.api as sm

tickers = clean['ticker'].unique()
ticker = tickers[0]

tmp = clean[clean['ticker']==ticker]
dates = tmp['date'].unique()
tmp = tmp[tmp['date']==dates[0]]
tmp.head(5)

exogs = ['const','delta_n','matur_n','delta2','inter']
mB = np.linalg.lstsq(tmp[exogs],tmp['imvol'])[0]

mB

array([     0.272,      0.176,     -0.170,      0.920,      0.054])

def CrossReg(df,exogs):
    mod = sm.OLS(df['imvol'],df[exogs]).fit()
    ind = pd.MultiIndex.from_product([['params','tvalues'],exogs])
    ret = pd.Series(index=ind)
    ret.loc['params'] = mod.params.values
    ret.loc['tvalues'] = mod.tvalues.values
    ret.loc['r2'] = mod.rsquared
    return ret 

exogs = ['const','delta_n','matur_n','delta2','inter']
mB = clean.groupby(['ticker','date']).apply(lambda x:CrossReg(x,exogs))

mB.loc['SPX']['params'].plot()

<matplotlib.axes._subplots.AxesSubplot at 0x7f2d109103c8>

def CrossRegSpeed(df,exogs):
    mB = np.linalg.lstsq(df[exogs],df['imvol'])[0]
    return pd.Series(mB,index=exogs)
mB = clean.groupby(['ticker','date']).apply(lambda x:CrossRegSpeed(x,exogs))

import multiprocessing as mp
from functools import partial

def CrossRegParHelper(ticker):
    tmp = DATA[DATA['ticker']==ticker]
    mB = tmp.groupby(['date']).apply(lambda x:CrossRegSpeed(x,EXOGS))
    mB['ticker'] = ticker
    return mB

def CrossRegPar(data,exogs):
    global DATA, EXOGS
    DATA = data 
    EXOGS = exogs
    pool = mp.Pool(processes=12)
    results = pool.map(CrossRegParHelper, tickers)
    pool.close()
    pool.join()
    mB = pd.concat(results)
    return mB

mB = CrossRegPar(clean,exogs)

%timeit clean.groupby(['ticker','date']).apply(lambda x:CrossRegSpeed(x,exogs))

8 s ± 33.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit CrossRegPar(clean,exogs)

4.2 s ± 51.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

import numba as nb

clean.sort_values(by=['ticker','date'],inplace=True)
locs = clean.groupby(['ticker','date']).size().values.cumsum()
locs = np.insert(locs,0,0)
mY = clean[['imvol']].values
mX = clean[exogs].values

@nb.jit
def CrossRegJIT(mX,mY,locs):
    N = locs.shape[0]
    mB = np.zeros((mX.shape[1],N-1))
    for i in range(0,N-1):
        b = np.linalg.lstsq(mX[locs[i]:locs[i+1]],mY[locs[i]:locs[i+1]])[0]
        mB[:,[i]] = b
    return mB

mB = CrossRegJIT(mX,mY,locs)

%timeit CrossRegJIT(mX,mY,locs)

702 ms ± 6.65 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

mB = pd.DataFrame(mB.T,index=clean.groupby(['ticker','date']).size().index,columns=exogs)

mB.head(5)

mB.reset_index(inplace=True)
level = pd.pivot_table(mB,index='date',columns='ticker',values='const')

level.plot()

from sklearn import preprocessing
import scipy as sp
def pcaAnalysis(mX, iPC, useCov=1):
    if useCov:
        # mCov = np.cov(mX, rowvar = 0)
        de_meanX = mX - np.tile(np.mean(mX, axis=0), [mX.shape[0], 1])
        mCov = np.dot(de_meanX.transpose(), de_meanX) / mX.shape[0]

        eigval, eigvec = sp.linalg.eigh(mCov)
        eigidx = np.argsort(eigval)
        maxidx = eigidx[::-1]
        sorted_vec = eigvec[:, maxidx]

        mPC = np.dot(de_meanX, sorted_vec)
        return mPC[:, 0:iPC], (eigval[maxidx] / np.sum(eigval))[0:iPC]

    else:
        mCoef = np.corrcoef(mX, rowvar=0)
        scaled_X = preprocessing.scale(mX)
        eigval, eigvec = sp.linalg.eigh(mCoef)

        eigidx = np.argsort(eigval)
        maxidx = eigidx[::-1]
        sorted_vec = eigvec[:, maxidx]

        mPC = np.dot(scaled_X, sorted_vec)
        return preprocessing.scale(mPC[:, 0:iPC]), (eigval[maxidx] /
                                                    np.sum(eigval))[0:iPC]

mPC,pct = pcaAnalysis(level.values[:500],3)
pct

array([  7.22e-01,   9.53e-02,   6.82e-02])

plt.plot(mPC)

[<matplotlib.lines.Line2D at 0x7f054028cb70>,
 <matplotlib.lines.Line2D at 0x7f054028c518>,
 <matplotlib.lines.Line2D at 0x7f054028cf28>]