import quandl
import cmocean
import warnings
import numpy as np
import pandas as pd
import seaborn as sns
import statsmodels.api as sm
from pytrends.request import TrendReq
from scipy import stats as ss
from datetime import datetime as dt
from arch.unitroot import PhillipsPerron, ADF
from matplotlib import pyplot as plt
from matplotlib import dates as mdates
from statsmodels.graphics.api import qqplot
pytrend = TrendReq('[email protected]', 'aaaaaaaaaa', hl='en-US', tz=360, custom_useragent=None)
warnings.filterwarnings('ignore')
quandl.ApiConfig.api_key = 'aaaaaaaaaaaaaaa'

pd.set_option("display.width",999)
plt.rcdefaults()
plt.style.use('seaborn-paper')
plt.rcParams['axes.facecolor']='w'
plt.rcParams['axes.grid']=False
plt.rcParams['font.family'] = 'serif'
plt.rcParams['figure.figsize'] = [12, 4]

exchangelist = [
    'BTCE',
    'KRAKEN',
    'BITFINEX',
    'BITSTAMP',
#     'COINBASE',
    'ROCK',
    'HITBTC',
    'ITBIT',
    'BITBAY',
    'BITKONAN',
    'CBX',
    'VCX'
]
n = len(exchangelist)
print('there are {} exchanges in total'.format(n))

exchanges = [quandl.get('BCHARTS/' + e + 'USD') for e in exchangelist]

for i in range(n):
    daterange = exchanges[i].index.tolist()
    print('{:>8}'.format(exchangelist[i]), daterange[0].date(), daterange[-1].date())

start_date = pd.to_datetime('2014-07-01')
mid_date = pd.to_datetime('2016-01-01')
end_date = pd.to_datetime('2016-07-01')

price = pd.concat([e.ix[start_date:end_date,'Close'] for e in exchanges], join='outer', axis=1)
volume = pd.concat([e.ix[start_date:end_date,'Volume (Currency)'] for e in exchanges], join='outer', axis=1)
price.columns = exchangelist
volume.columns = exchangelist
price = price[price > 0].fillna(method='ffill').fillna(method='bfill')
volume = price[volume > 0].fillna(method='ffill').fillna(method='bfill')
print(price.describe())

%config InlineBackend.figure_format = 'retina'
min_exchange = pd.Series([np.argmin(price.iloc[i], axis=0) for i in range(len(price))],
                         index=price.index).apply(lambda fn: exchangelist.index(fn))
max_exchange = pd.Series([np.argmax(price.iloc[i], axis=0) for i in range(len(price))],
                         index=price.index).apply(lambda fn: exchangelist.index(fn))
dmin = price.min(axis=1)
dmax = price.max(axis=1)

thermal = cmocean.cm.thermal

plt.close()
fig = plt.figure()

ax = fig.add_subplot(111)
for i in range(n):
    ax.plot((price.index[0], price.index[0]), (0, 0), c=thermal(i/n), label=exchangelist[i])
for i in range(len(price)):
    ax.plot((price.index[i], price.index[i]), (0, -dmin[i]), lw=1.5, c=thermal(min_exchange[i]/n))
    ax.plot((price.index[i], price.index[i]), (0, dmax[i]), lw=1.5, c=thermal(max_exchange[i]/n))
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=2))
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y\n %b'))
ax.set_xlim(start_date, end_date)
ax.set_ylim(-1500,1500)

plt.legend(loc='lower center', frameon=0, ncol=n//2)
plt.tight_layout()
plt.savefig('price1.pdf')

plt.show()

active_price = price * (price.diff() != 0)
active_price = active_price[active_price > 0]
err_count = np.isnan(active_price).sum()

%config InlineBackend.figure_format = 'retina'
plt.close()
fig = plt.figure()

ax1 = fig.add_subplot(111)
(100*max_exchange.value_counts().sort_index().reindex(range(len(exchangelist))).fillna(0)/len(price)).plot.bar(ax=ax1, color='c', align='center', alpha=.5, label='% being highest')
(100*min_exchange.value_counts().sort_index().reindex(range(len(exchangelist))).fillna(0)/len(price)).plot.bar(ax=ax1, color='r', align='center', alpha=.5, label='% being lowest')
ax1.set_ylim(-5,80)
ax1.set_ylabel('%')
plt.xticks(range(len(err_count)), err_count.index.tolist(), rotation=0)
ax1.legend(loc='upper left', frameon=0)

ax2 = ax1.twinx()
err_count.plot(ax=ax2, marker='o', color='r', alpha=.5, label='total number of error days')
ax2.set_xlim(-.5,len(err_count)-.5)
ax2.set_xticks(range(len(err_count)), err_count.index.tolist())
ax2.legend(loc='upper right', frameon=0)
ax2.set_ylim(-35,500)

# plt.grid(None)
plt.tight_layout()
plt.show()

irr_exchanges = [10,9,8,7,5,0]  # delete those that have over 200 error days, or dominate over 10% of all time range
n -= len(irr_exchanges)
active_price.drop(active_price.columns[[irr_exchanges]], axis=1, inplace=1)
for i in irr_exchanges:
    exchangelist.pop(i)

%config InlineBackend.figure_format = 'retina'
min_exchange = pd.Series([np.argmin(active_price.iloc[i], axis=0) for i in range(len(active_price))], index=active_price.index).apply(lambda fn: exchangelist.index(fn))
dmin = active_price.min(axis=1)
max_exchange = pd.Series([np.argmax(active_price.iloc[i], axis=0) for i in range(len(active_price))], index=active_price.index).apply(lambda fn: exchangelist.index(fn))
dmax = active_price.max(axis=1)

thermal = cmocean.cm.thermal

plt.close()
fig = plt.figure()

ax = fig.add_subplot(111)
for i in range(n):
    ax.plot((active_price.index[0], active_price.index[0]), (0, 0), c=thermal(i/n), label=exchangelist[i])
for i in range(len(active_price)):
    ax.plot((active_price.index[i], active_price.index[i]), (0, -dmin[i]), lw=1.5, c=thermal(min_exchange[i]/n))
    ax.plot((active_price.index[i], active_price.index[i]), (0, dmax[i]), lw=1.5, c=thermal(max_exchange[i]/n))
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=2))
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y\n %b'))
ax.set_xlim(start_date, end_date)
ax.set_ylim(-1500,1500)

plt.legend(loc='lower center', frameon=0, ncol=n)
plt.tight_layout()
plt.savefig('price2.pdf')

plt.show()

%config InlineBackend.figure_format = 'retina'

err_count = np.isnan(active_price).sum()

plt.close()
fig = plt.figure()

ax1 = fig.add_subplot(111)
(100*max_exchange.value_counts().sort_index().reindex(range(len(exchangelist))).fillna(0)/len(price)).plot.bar(ax=ax1, color='c', align='center', alpha=.5, label='% being highest')
(100*min_exchange.value_counts().sort_index().reindex(range(len(exchangelist))).fillna(0)/len(price)).plot.bar(ax=ax1, color='r', align='center', alpha=.5, label='% being lowest')
ax1.set_ylim(-5,80)
ax1.set_ylabel('%')
plt.xticks(range(len(err_count)), err_count.index.tolist(), rotation=0)
ax1.legend(loc='upper left', frameon=0)

ax2 = ax1.twinx()
err_count.plot(ax=ax2, marker='o', color='r', alpha=.5, label='total number of error days')
ax2.set_xlim(-.5,len(err_count)-.5)
ax2.set_xticks(range(len(err_count)), err_count.index.tolist())
ax2.legend(loc='upper right', frameon=0)
ax2.set_ylim(-7,100)

# plt.grid(None)
plt.tight_layout()
plt.show()

%config InlineBackend.figure_format = 'retina'
plt.close()

fig = plt.figure()

ax = fig.add_subplot(111)
for i in range(n):
    ax.plot(active_price.ix[:,i], c=thermal(i/n), label=exchangelist[i])
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=2))
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y\n %b'))
ax.axvspan(mid_date, end_date, color='c', alpha=0.05)
ax.set_xlim(start_date, end_date)
ax.set_ylim(0,1000)

plt.legend(loc='lower center', frameon=0, ncol=n)
plt.tight_layout()
# plt.title('Historical Prices')
plt.savefig('price3.pdf')
plt.show()

price_spread = active_price.max(axis=1) - active_price.min(axis=1)
price_spread.describe()

%config InlineBackend.figure_format = 'retina'
plt.close()

fig = plt.figure()

ax = fig.add_subplot(111)
ax.plot(price_spread, lw=1,c='k')
ax.axhline(price_spread.mean(), alpha=0.5, ls='--', c='c', label='$\mu$')
ax.axhline(price_spread.mean()+price_spread.std()*3, alpha=0.5, ls='--', c='r', label='$\pm3\sigma$')
ax.axhline(price_spread.mean()-price_spread.std()*3, alpha=0.5, ls='--', c='r')
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=2))
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y\n %b'))
ax.axvspan(mid_date, end_date, color='c', alpha=0.05)
ax.set_xlim(start_date, end_date)
ax.set_ylim(-60,100)

plt.tight_layout()
plt.legend(loc='lower center', framealpha=0, ncol=2)
# plt.title('Price Spreads')
plt.savefig('spread.pdf')
plt.show()

log_price_spread = np.log(price_spread)
log_price_spread.describe()

%config InlineBackend.figure_format = 'retina'
plt.close()
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(log_price_spread, c='k',lw=1,label='')
ax.axhline(log_price_spread.mean(), alpha=0.5, ls='--', c='c', label='$\mu$')
ax.axhline(log_price_spread.mean()+log_price_spread.std()*3, alpha=0.5, ls='--', c='r',label='$\pm3\sigma$')
ax.axhline(log_price_spread.mean()-log_price_spread.std()*3, alpha=0.5, ls='--', c='r')
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=2))
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y\n %b'))
ax.axvspan(mid_date, end_date, color='c', alpha=0.05)
ax.set_ylim(-3,6)
ax.set_xlim(start_date, end_date)
plt.tight_layout()
plt.legend(loc='lower center', framealpha=0, ncol=2)
# plt.title('Log Price Spreads')
plt.savefig('logspread.pdf')
plt.show()

print(ADF(log_price_spread, trend='c'), '\n')
print(PhillipsPerron(log_price_spread, trend='c'), '\n')

def z_norm(s):
    return (s - s.mean(axis=0)) / s.std(axis=0)

log_price_spread = z_norm(log_price_spread.rolling(5).mean().fillna(method='bfill'))
log_price_spread.describe()

print(ADF(log_price_spread, trend='nc'), '\n')
print(PhillipsPerron(log_price_spread, trend='nc'), '\n')

%config InlineBackend.figure_format = 'retina'
plt.close()
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(log_price_spread, c='k',lw=1,label='')
ax.axhline(log_price_spread.mean(), alpha=0.5, ls='--', c='c', label='$\mu$')
ax.axhline(log_price_spread.mean()+log_price_spread.std()*3, alpha=0.5, ls='--', c='r',label='$\pm3\sigma$')
ax.axhline(log_price_spread.mean()-log_price_spread.std()*3, alpha=0.5, ls='--', c='r')
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=2))
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y\n %b'))
ax.axvspan(mid_date, end_date, color='c', alpha=0.05)
ax.set_ylim(-5,5)
ax.set_yticks(range(-5,6))
ax.set_xlim(start_date, end_date)
plt.tight_layout()
plt.legend(loc='lower center', framealpha=0, ncol=2)
# plt.title('Log Price Spreads')
plt.savefig('logspreadnorm.pdf')
plt.show()

train = log_price_spread.loc[:mid_date]
test = log_price_spread.loc[mid_date:]
train.name = 'LogPs'

%config InlineBackend.figure_format = 'retina'
plt.close()
fig = plt.figure(figsize=(8,4))

ax1 = fig.add_subplot(211)
sm.graphics.tsa.plot_acf(train, lags=40, ax=ax1, label='ACF')
ax1.set_ylim(-.55,1.25)
ax1.set_yticks(np.arange(-.5,1.5,.5))
ax1.set_xticks(range(0,40,5))
ax1.set_title('')
handles, labels= ax1.get_legend_handles_labels()
handles = handles[:-len(handles)//3][1::2]
labels = labels[:-len(handles)//3][1::2]
ax1.legend(handles=handles, labels=labels,loc='best', frameon=0, numpoints=2)

ax2 = fig.add_subplot(212)
sm.graphics.tsa.plot_pacf(train, lags=40, ax=ax2, label='PACF')
ax2.set_ylim(-.55,1.25)
ax2.set_yticks(np.arange(-.5,1.5,.5))
ax2.set_xticks(range(0,40,5))
ax2.set_title('')
handles, labels= ax2.get_legend_handles_labels()
handles = handles[:-len(handles)//3][1::2]
labels = labels[:-len(handles)//3][1::2]
ax2.legend(handles=handles, labels=labels,loc='best', frameon=0, numpoints=2)

plt.tight_layout()
plt.savefig('uniacfpacf.pdf')
plt.show()

min_aic = min_bic = 1e5
min_aic_order = min_bic_order = (0,0)
for p in range(1,10):
    for q in range(6):
        model = sm.tsa.ARMA(train, (p, q))
        try:
            result = model.fit(trend='nc')
            stat = abs(sm.stats.durbin_watson(result.resid)-2)
            stat_score = (stat <= .01) * 3 + (.01 < stat <= .05) * 2 + (.05 < stat <= .10) * 1
            if np.abs(np.concatenate([result.arroots, result.maroots])).min() > 1.05:
                aic = result.aic
                bic = result.bic
                if aic < min_aic: 
                    min_aic = aic
                    min_aic_order = (p,q)
                if bic < min_bic: 
                    min_bic = bic
                    min_bic_order = (p,q)
                print('({},{}) AIC = {:>7.2f}, BIC = {:>7.2f}, |DW-2| = {:.8f} {}'.format(p,q,aic,bic,stat,stat_score*'*'))
        except:
            continue
print()
print('Minimum AIC = {:.2f} at {}, minimum BIC = {:.2f} at {}'.format(min_aic, min_aic_order, min_bic, min_bic_order))

model = sm.tsa.ARMA(train, (2,0))
result = model.fit(trend='nc')
print(result.summary())

sm.stats.jarque_bera(result.resid)[1]

sm.stats.acorr_ljungbox(result.resid, lags=1)[1][0]

%config InlineBackend.figure_format = 'retina'
plt.close()
fig = plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
fig = qqplot(result.resid, line='45', ax=ax, fit=True)
lines = fig.findobj(lambda x: hasattr(x, 'get_color') and x.get_color() == 'r')
dots = fig.findobj(lambda x: hasattr(x, 'get_color') and x.get_color() == 'b')
[l.set_linewidth(2) for l in lines]
[d.set_alpha(0.3) for d in dots]
# [d.set_color('c') for d in dots]
[d.set_markersize(5) for d in dots]
ax.set_ylim(-3,3)
ax.set_yticks(range(-3,4))
ax.set_xlim(-3,3)
ax.set_xticks(range(-3,4))
plt.tight_layout()
# plt.title('p-value = {}'.format(ss.normaltest(result.resid)[1]))
plt.savefig('uniqqplot.pdf')
plt.show()

phi = result.arparams * (result.pvalues[:result.k_ar]<=0.05).values

def yhat(y):  # y is now nparray of shape (k_ar,2)
    ylag = y[::-1]
    return ylag.T @ phi

fit = result.fittedvalues

pred = test * np.nan
for i in range(len(phi)):
    pred[i] = yhat(train[-1-i-len(phi):-1-i])
for i in range(len(phi), len(test)):
    pred[i] = yhat(test[i-len(phi):i])

fcst = result.forecast(len(test))
forecast = pd.Series(fcst[0], index=test.index)
conf_interval = pd.DataFrame(fcst[2], index=test.index)

%config InlineBackend.figure_format = 'retina'
plt.close()
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(log_price_spread, c='k',lw=1, label='true data')
ax.plot(fit, 'ro-', markevery=10, alpha=0.5, label='insample fitted data')
# ax.plot(result.fittedvalues, 'b-')
ax.plot(pred, 'co-', markevery=10, alpha=0.5, label='outsample point forecast')
ax.plot(forecast, 'c--', label='outsample forecast')
ax.axvspan(mid_date, end_date, color='c', alpha=0.05)
ax.fill_between(conf_interval.index, conf_interval.ix[:,0], conf_interval.ix[:,1], color='k', alpha=.1, label='95% conf. interval')
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=2))
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y\n %b'))
ax.set_xlim(start_date, end_date)
ax.set_ylim(-5,5)
ax.set_yticks(range(-5,6))
plt.tight_layout()
plt.legend(loc='lower center', framealpha=0, ncol=5)
# plt.title('First Difference Log Price Spreads')
plt.savefig('unifit1.pdf')
plt.show()

temp_date = pd.to_datetime('2015-12-01')

%config InlineBackend.figure_format = 'retina'
plt.close()
fig = plt.figure()

ax = fig.add_subplot(111)
ax.plot(log_price_spread, c='k',lw=1, label='true data')
ax.plot(fit, 'ro-', markevery=5, alpha=0.5, label='insample fitted data')
ax.plot(pred, 'co-', markevery=5, alpha=0.5, label='outsample point forecast')
ax.plot(forecast, 'c--', alpha=0.5, label='outsample forecast')
ax.axvspan(mid_date, end_date, color='c', alpha=0.05)
ax.fill_between(conf_interval.index, conf_interval.ix[:,0], conf_interval.ix[:,1], color='k', alpha=.1, label='95% conf. interval')
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=1))
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y\n %b'))
ax.set_xlim(temp_date, end_date)
ax.set_ylim(-5,5)
ax.set_yticks(range(-5,6))
plt.tight_layout()
plt.legend(loc='lower center', framealpha=0, ncol=5)
# plt.title('First Difference Log Price Spreads')
plt.savefig('unifit2.pdf')
plt.show()

def mase(data, pred):
    error = data - pred
    ydiff = data.diff().dropna()
    return error.abs().mean() / ydiff.abs().mean()

mase_uni_is = mase(train, fit)
mase_uni_osp = mase(test, pred)
mase_uni_os = mase(test, forecast)
print('MASE insample is {:.3f}'.format(mase_uni_is))
print('MASE outsample point is {:.3f}'.format(mase_uni_osp))
print('MASE outsample is {:.3f}'.format(mase_uni_os))

data = [
    'BCHAIN/ATRCT',     # Bitcoin Median Transaction Confirmation Time
    'BCHAIN/CPTRA',     # Bitcoin Cost Per Transaction
    'BCHAIN/DIFF',      # Bitcoin Mining Difficulty
    'BCHAIN/TVTVR',     # Bitcoin Trade Volume vs Transaction Volume Ratio
]

ind = [quandl.get(d)[start_date:end_date] for d in data]

max_volume = ((volume.apply(lambda v:v.index, axis=1).T == max_exchange.apply(lambda i:exchangelist[i])).T * volume).sum(axis=1)
min_volume = ((volume.apply(lambda v:v.index, axis=1).T == min_exchange.apply(lambda i:exchangelist[i])).T * volume).sum(axis=1)
trade_volume = max_volume + min_volume
transaction_volume = (trade_volume * ind[-1].T).T
ind[-1] = transfer_volume = transaction_volume.Value - trade_volume

dataset = pd.concat([price_spread, trade_volume] + ind, axis=1, join='outer')
dataset.columns = ['Ps', 'Tr', 'Ct', 'Co', 'Di', 'Tf']
print(dataset.describe())

googletrends = []
kw_list = ['bitcoin robot', 'bitcoin arbitrage', 'bitcoin spread']
pytrend.build_payload(kw_list=kw_list)
googletrends = pytrend.interest_over_time().ix[start_date:end_date, :] + 1  # as googletrends can be 0 from time to time, which is troublesome for taking log

plt.close()
fig = plt.figure()
ax = fig.add_subplot(111)
for i in range(len(kw_list)):
    ax.plot(googletrends.ix[:,i], c=cmocean.cm.algae(i/len(kw_list)), label=kw_list[i])
ax.plot(googletrends.sum(axis=1), c='c', label='all')
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=2))
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y\n %b'))
ax.set_xlim(start_date, end_date)
plt.tight_layout()
plt.legend(loc='lower center', frameon=0, ncol=len(kw_list)+1)
plt.ylim(-10,100)
plt.yticks(range(0,101,10))
plt.savefig('multrend.pdf')
plt.show()

exchangetrends = []
kw_list = np.setdiff1d(exchangelist, ['ROCK'])
pytrend.build_payload(kw_list=kw_list)
exchangetrends = pytrend.interest_over_time().ix[start_date:end_date, :] + 1

plt.close()
fig = plt.figure()
ax = fig.add_subplot(111)
for i in range(len(kw_list)):
    ax.plot(exchangetrends.ix[:,i], color=cmocean.cm.algae(i/len(kw_list)), label=kw_list[i])
ax.plot(exchangetrends.sum(axis=1), c='c', label='all')
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=2))
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y\n %b'))
ax.set_xlim(start_date, end_date)
plt.tight_layout()
plt.legend(loc='lower center', frameon=0, ncol=len(kw_list)+1)
plt.ylim(-10,100)
plt.yticks(range(0,101,10))
plt.savefig('multrend.pdf')
plt.show()

dataset = pd.concat([dataset, googletrends.sum(axis=1), exchangetrends.sum(axis=1)], axis=1, join='outer').fillna(method='ffill').fillna(method='bfill')
dataset.columns = list(dataset.columns[:-2]) + ['Gt1', 'Gt2']
dataset.Gt1.interpolate(inplace=1)
dataset.Gt2.interpolate(inplace=1)

log_dataset = z_norm(np.log(dataset).rolling(5).mean().fillna(method='bfill'))
log_dataset.columns = ['Log'+c for c in log_dataset.columns]
print(log_dataset.describe())

%config InlineBackend.figure_format = 'retina'
plt.close()

fig = plt.figure(figsize=(5,4))

ax = fig.add_subplot(111)
idx = log_dataset.corr().sort_values(by='LogPs', ascending=0).index
log_dataset = log_dataset.ix[:,idx]
corr = log_dataset.corr()
sns.heatmap(corr, square=1, ax=ax)

plt.yticks(rotation=0)
plt.xticks(rotation=90)

plt.tight_layout()
plt.savefig('mulcorr.pdf')
plt.show()

fig = plt.figure(figsize=(12,4*len(corr)))
for i in range(len(corr)-1):
    ax = fig.add_subplot(len(corr)-1,1,i+1)
    ax.plot(log_price_spread, c='r', label='LogPs', alpha=.5)
    ax.plot(log_dataset.ix[:,i+1], c='c', alpha=.5, label=log_dataset.columns[i+1])
    ax.legend(frameon=0,loc='lower center', ncol=2)
    ax.xaxis.set_major_locator(mdates.MonthLocator(interval=2))
    ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y\n %b'))
    ax.set_xlim(start_date, end_date)
    ax.set_ylim(-5,5)
    ax.set_yticks(range(-5,6))
plt.tight_layout()
plt.show()

endog = log_dataset.LogPs
exog = log_dataset.ix[:,np.setdiff1d(log_dataset.columns,['LogPs'])]
train_endog = endog.ix[start_date:mid_date]
train_exog = exog.ix[start_date:mid_date]
test_endog = endog.ix[mid_date:end_date]
test_exog = exog.ix[mid_date:end_date]

min_aic = min_bic = 1e5
min_aic_order = min_bic_order = (0,0,0)
for p in range(1,10):
    for q in range(3):
        for d in range(1):
            model = sm.tsa.SARIMAX(train_endog, train_exog, (p, d, q))
            try:
                result = model.fit(trend='nc')
                result.summary()
            except:
                continue
            stat = abs(sm.stats.durbin_watson(result.resid)-2)
            stat_score = (stat <= 0.01) * 3 + (0.01 < stat <= 0.05) * 2 + (0.05 < stat <= 0.10) * 1
            if np.abs(np.concatenate([result.arroots, result.maroots])).min() > 1.05:
                aic = result.aic
                bic = result.bic
                if aic < min_aic: 
                    min_aic = aic
                    min_aic_order = (p,d,q)
                if bic < min_bic: 
                    min_bic = bic
                    min_bic_order = (p,d,q)
                print('({},{},{}) AIC = {:>7.2f}, BIC = {:>7.2f}, |DW-2| = {:.8f} {}'.format(p,d,q,aic,bic,stat,stat_score*'*'))
print()
print('Minimum AIC = {:.2f} at {}, minimum BIC = {:.2f} at {}'.format(min_aic, min_aic_order, min_bic, min_bic_order))

model = sm.tsa.SARIMAX(train_endog, train_exog, (8,0,1), freq='D')
result = model.fit(trend='nc')
print(result.summary())

sm.stats.durbin_watson(result.resid)

sm.stats.jarque_bera(result.resid)[1]

sm.stats.acorr_ljungbox(result.resid, lags=1)[1][0]

%config InlineBackend.figure_format = 'retina'
plt.close()
fig = plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
fig = qqplot(result.resid, line='45', ax=ax, fit=True)
lines = fig.findobj(lambda x: hasattr(x, 'get_color') and x.get_color() == 'r')
dots = fig.findobj(lambda x: hasattr(x, 'get_color') and x.get_color() == 'b')
[l.set_linewidth(2) for l in lines]
[d.set_alpha(0.3) for d in dots]
[d.set_markersize(5) for d in dots]
ax.set_ylim(-3,3)
ax.set_yticks(range(-3,4))
ax.set_xlim(-3,3)
ax.set_xticks(range(-3,4))
plt.tight_layout()
# plt.title('p-value = {}'.format(ss.normaltest(result.resid)[1]))
plt.savefig('mulqqplot.pdf')
plt.show()

forecast = result.forecast(steps=len(test_endog), exog=test_exog)
fit = result.fittedvalues

model = sm.tsa.statespace.SARIMAX(endog, exog, model.order)
result = model.filter(result.params)
predict = result.get_prediction()
pred = predict.predicted_mean[mid_date:]
conf_interval = [forecast+predict.se_mean[len(train_endog):]*1.96,forecast-predict.se_mean[len(train_endog):]*1.96]

%config InlineBackend.figure_format = 'retina'
plt.close()
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(log_price_spread, c='k',lw=1, label='true data')
ax.plot(fit, 'ro-', markevery=10, alpha=0.5, label='insample fitted data')
ax.plot(pred, 'co-', markevery=10, alpha=0.5, label='outsample point forecast')
ax.plot(forecast, 'c--', label='outsample forecast')
ax.axvspan(mid_date, end_date, color='c', alpha=0.05)
ax.fill_between(test_endog.index, conf_interval[0], conf_interval[1], color='k', alpha=.1, label='95% conf. interval')
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=2))
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y\n %b'))
ax.set_xlim(start_date, end_date)
ax.set_ylim(-5,5)
ax.set_yticks(range(-5,6))
plt.tight_layout()
plt.legend(loc='lower center', framealpha=0, ncol=5)
# plt.title('First Difference Log Price Spreads')
plt.savefig('mulfit1.pdf')
plt.show()

%config InlineBackend.figure_format = 'retina'
plt.close()
fig = plt.figure()

ax = fig.add_subplot(111)
ax.plot(log_price_spread, c='k',lw=1, label='true data')
ax.plot(fit, 'ro-', markevery=5, alpha=0.5, label='insample fitted data')
ax.plot(pred, 'co-', markevery=5, alpha=0.5, label='outsample point forecast')
ax.plot(forecast, 'c--', alpha=0.5, label='outsample forecast')
ax.axvspan(mid_date, end_date, color='c', alpha=0.05)
ax.fill_between(test_endog.index, conf_interval[0], conf_interval[1], color='k', alpha=.1, label='95% conf. interval')
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=1))
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y\n %b'))
ax.set_xlim(temp_date, end_date)
ax.set_ylim(-5,5)
ax.set_yticks(range(-5,6))
plt.tight_layout()
plt.legend(loc='lower center', framealpha=0, ncol=5)
# plt.title('First Difference Log Price Spreads')
plt.savefig('mulfit2.pdf')

plt.show()

mase_mul_is = mase(train_endog, fit)
mase_mul_osp = mase(test_endog, pred)
mase_mul_os = mase(test_endog, forecast)
print('MASE insample is {:.3f}'.format(mase_mul_is))
print('MASE outsample is {:.3f}'.format(mase_mul_osp))
print('MASE outsample is {:.3f}'.format(mase_mul_os))

uni = np.array([mase_uni_is, mase_uni_osp, mase_uni_os])
mul = np.array([mase_mul_is, mase_mul_osp, mase_mul_os])
print(uni)
print(mul)

plt.figure(figsize=(8,4))
plt.plot((0,0),(0,5),'#eeeeee',ls='--')
plt.plot((1,1),(0,5),'#eeeeee',ls='--')
plt.plot((2,2),(0,5),'#eeeeee',ls='--')
plt.plot(uni, 'ro-', label='univariate model')
plt.plot(mul, 'ko-', label='multivariate model')
plt.legend(loc='lower center', frameon=0, ncol=2)
plt.xticks(range(3),['insample', 'outsample point', 'outsample'])
# plt.title('MASE Scores')
plt.yticks(np.arange(0,6))
plt.ylim(0,5)
plt.tight_layout()
plt.savefig('mase.pdf')
plt.show()