In [1]:
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sbn
import numpy as np
from warnings import filterwarnings
from helpers import HarModel, HarchModel
from itertools import product
from pandas.tools.plotting import autocorrelation_plot
from scipy.stats import expon, laplace, t
from arch.univariate import HARX, HARCH, GARCH, EGARCH, StudentsT, SkewStudent
from statsmodels.graphics.tsaplots import plot_pacf
In [2]:
filterwarnings('ignore')
%matplotlib notebook
sbn.set_style('whitegrid')
In [3]:
dat_file = './data/CommodityIndices.xlsx'
skiprows = (2, 3, 4, 5)
raw_data = pd.read_excel(dat_file, skiprows=skiprows, index_col=0, parse_cols='B:AE')
raw_data.columns = [column[4:-2].rstrip() for column in raw_data.columns]
raw_data.iloc[0] = 100.0
In [4]:
raw_data.columns.values.tolist()
Out[4]:
In [5]:
name = 'Gold' # Name of index as string
data = raw_data[name].resample('BM').mean()
fig, ax = plt.subplots(num=name+' Over the Aeons')
data.plot(ax=ax, label=name)
ax.legend(loc='best', frameon=True, fancybox=True, framealpha=0.8)
ax.set_xlabel('Time')
ax.set_ylabel('Relative index')
fig.tight_layout()
We don't want to do any smoothing or averaging to find a trend in the data because we know little of how to extrapolate such a trend to the future.
In [6]:
winsize = 12 # Window size (in business months) to average over
max_lag = 24 # Maximum lag to consider in (partial) autocorrelation
returns = 100 * data.pct_change(1).dropna()
rtnswin = returns.rolling(winsize, center=True) # or .ewm(com=winsize)
rtnmean = rtnswin.mean()
rtn_std = rtnswin.std()
fig = plt.figure(num='Returns on '+name+' and their (Partial) Autocorrelation', figsize=(9.7, 10))
return_ax = fig.add_subplot(311)
aucorr_ax = fig.add_subplot(312)
pacorr_ax = fig.add_subplot(313)
returns.plot(color='#4c72b0', style=['_'], marker='o', alpha=0.3, ax=return_ax)
rtnmean.plot(color='#4c72b0', ax=return_ax)
(rtnmean + rtn_std).plot(color='#c44e52', linewidth=1.0, ax=return_ax)
(rtnmean - rtn_std).plot(color='#c44e52', linewidth=1.0, ax=return_ax)
return_ax.set_xlabel('Time')
return_ax.set_ylabel('Monthly returns')
autocorrelation_plot(returns, linewidth=1.0, ax=aucorr_ax)
aucorr_ax.set_xlim(left=0, right=max_lag)
plot_pacf(returns, lags=max_lag, alpha=0.01, ax=pacorr_ax, linewidth=0.5)
pacorr_ax.set_ylim(aucorr_ax.get_ylim())
pacorr_ax.set_xlim(aucorr_ax.get_xlim())
pacorr_ax.set_title('')
pacorr_ax.set_xlabel('Lag')
pacorr_ax.set_ylabel('Partial autocorrelation')
fig.tight_layout()
We note that:
So let's look at the absolute returns and their autocorrelation next.
In [7]:
fig = plt.figure(num='Absolute Returns on '+name+' and their (Partial) Autocorrelation', figsize=(9.7, 10))
return_ax = fig.add_subplot(311)
aucorr_ax = fig.add_subplot(312)
pacorr_ax = fig.add_subplot(313)
returns.abs().plot(linewidth=0.5, ax=return_ax, label='Volatility')
return_ax.set_xlabel('Time')
return_ax.set_ylabel('Absolute monthly returns')
autocorrelation_plot(returns.abs(), linewidth=0.5, ax=aucorr_ax, label='Volatility');
aucorr_ax.set_xlim(left=0, right=max_lag)
plot_pacf(returns.abs(), lags=max_lag, alpha=0.01, ax=pacorr_ax, linewidth=0.5)
pacorr_ax.set_ylim(aucorr_ax.get_ylim())
pacorr_ax.set_xlim(aucorr_ax.get_xlim())
pacorr_ax.set_title('')
pacorr_ax.set_xlabel('Lag')
pacorr_ax.set_ylabel('Partial autocorrelation')
fig.tight_layout()
We observe volatility bunching and long-term autocorrelation of the absolute returns. The partial autocorrelation dies off at much shorter lags.
In [8]:
fig = plt.figure(num='Distribution of (Absolute) Returns on '+name, figsize=(9.7, 4))
returns_ax = fig.add_subplot(121)
absrtns_ax = fig.add_subplot(122)
sbn.distplot(returns,
bins=19,
kde=False,
fit=laplace,
fit_kws={'color': '#c44e52', 'label': 'Laplace'},
ax=returns_ax);
returns_ax.set_title('Monthly Returns on '+name)
returns_ax.set_xlabel('Value')
returns_ax.set_ylabel('Frequency')
sbn.distplot(returns.abs(),
bins=19,
kde=False,
fit=expon,
fit_kws={'color': '#c44e52'},
color='#55a868',
ax=absrtns_ax);
absrtns_ax.set_title('Absolute Monthly Returns on '+name)
absrtns_ax.set_xlabel('Absolute Value')
absrtns_ax.set_ylabel('Frequency')
x_lims = returns_ax.get_xlim()
x_vals = np.linspace(*x_lims, 200)
params = t.fit(returns[returns.abs() > 0.05], floc=0)
returns_ax.plot(x_vals, t.pdf(x_vals, *params), label='Students-t')
returns_ax.legend(loc='best', frameon=True, fancybox=True, framealpha=0.8)
fig.tight_layout()
They returns follow a Laplace- whereas the absolute returns follow an Exponential (or maybe also a Pareto) distribution. The tails are far too fat for a Gaussian process. Interestingly, a Students-t distribution appears to be quite a reasonable fit for the returns, if the pronounced peak at zero is neglected.
In [9]:
params = { 'y': returns,
'constant': False,
'hold_back': max_lag,
'distribution': StudentsT()}
har = HarModel(params)
In [10]:
while har.bic.gets_smaller():
har.refine()
har.extend()
result = har.model.fit()
result.summary()
Out[10]:
In [11]:
harch = HarchModel(har.model)
while harch.bic.gets_smaller():
harch.refine()
harch.extend()
result = harch.model.fit()
result.summary()
Out[11]:
In [12]:
fixed = { 'y': returns,
'lags': har.lags,
'hold_back': max_lag}
tests = { 'constant': [False, True],
'distribution': [StudentsT(), SkewStudent()],
'volatility': [HARCH(lags=harch.lags),
GARCH(p=1, o=0, q=1),
GARCH(p=1, o=1, q=1),
GARCH(p=1, o=0, q=1, power=1.0),
GARCH(p=1, o=1, q=1, power=1.0),
EGARCH(p=1, o=0, q=1),
EGARCH(p=1, o=1, q=1)]}
fit_kws = {'update_freq': 0,
'disp': 'off'}
params_grid = (dict(zip(tests.keys(), values)) for values in product(*tests.values()))
min_bic, min_params = min((HARX(**fixed, **params).fit(**fit_kws).bic, params) for params in params_grid)
print('Minimum BIC is: ', min_bic)
print('for model with: ', [str(param) for param in min_params.values()])
model = HARX(**fixed, **min_params)
result = model.fit(**fit_kws)
result.summary()
Out[12]:
In [13]:
fig = plt.figure(num='Residuals of HAR+Volatility Model and their (Partial) Autocorrelation', figsize=(9.7, 10))
residual_ax = fig.add_subplot(311)
autocorr_ax = fig.add_subplot(312)
partcorr_ax = fig.add_subplot(313)
returns.plot(ax=residual_ax, linewidth=1, label='Returns')
result.resid.plot(ax=residual_ax, linewidth=0.5, label='Residual')
(returns-result.resid).plot(ax=residual_ax, color='#c44e52', label='Model')
residual_ax.set_xlabel('Time')
residual_ax.set_ylabel('Residual of model')
residual_ax.legend(loc='best', frameon=True, fancybox=True, framealpha=0.8)
autocorrelation_plot(returns, ax=autocorr_ax, linewidth=1, label='Returns')
autocorrelation_plot(result.resid.dropna(), ax=autocorr_ax, linewidth=1, label='Residual')
autocorr_ax.set_xlim(left=0, right=max_lag)
autocorr_ax.legend(loc='best', frameon=True, fancybox=True, framealpha=0.8)
plot_pacf(returns, lags=max_lag, alpha=0.01, ax=partcorr_ax, linewidth=0.5, label='Returns')
plot_pacf(result.resid.dropna(), lags=max_lag, alpha=0.05, ax=partcorr_ax, linewidth=0.5, label='Residual')
partcorr_ax.set_ylim(autocorr_ax.get_ylim())
partcorr_ax.set_xlim(autocorr_ax.get_xlim())
partcorr_ax.set_title('')
partcorr_ax.set_xlabel('Lag')
partcorr_ax.set_ylabel('Partial autocorrelation')
fig.tight_layout()
In [14]:
resids = (returns.abs()-result.conditional_volatility).dropna()
resids.plot(linewidth=0.5, ax=return_ax, label='Residual')
result.conditional_volatility.plot(linewidth=1, ax=return_ax, label='Model')
return_ax.legend(loc='best', frameon=True, fancybox=True, framealpha=0.8)
autocorrelation_plot(resids, linewidth=1, ax=aucorr_ax, label='Residual')
aucorr_ax.legend(loc='best', frameon=True, fancybox=True, framealpha=0.8)
plot_pacf(resids, lags=max_lag, alpha=0.05, ax=pacorr_ax, linewidth=0.5)
In [15]:
n_sims = 1000 # Number of simulations to average over
future = 12 # Number of future business months to provide a forecast for
past = 6 # Number of business months in the past to start our the forecast at
start_d = returns.index[-past-1]
horizon = past + future
forecast = result.forecast(horizon=horizon,
start=start_d,
align='origin',
method='simulation',
simulations=n_sims)
fc_index = pd.DatetimeIndex(freq='BM', start=start_d + pd.tseries.offsets.BMonthEnd(1), periods=horizon)
fc_frame = pd.DataFrame(data=forecast.mean.iloc[-past-1:].values[0], index=fc_index, columns=['mean'])
fc_frame['+var'] = forecast.variance.iloc[-past-1:].values[0]
fc_frame['-var'] = -forecast.variance.iloc[-past-1:].values[0]
fig, ax = plt.subplots(num='Forecast on Returns of '+name)
returns.iloc[-horizon:].plot(ax=ax)
fc_frame['mean'].plot(ax=ax)
fc_frame['+var'].plot(ax=ax, color='#c44e52')
fc_frame['-var'].plot(ax=ax, color='#c44e52')
ax.fill_between(fc_index, fc_frame['+var'], fc_frame['-var'], color='#c44e52', alpha=0.2)
ax.set_ylabel('Monthly returns')
fig.tight_layout()
In [16]:
past_pad = max_lag # Number of business months in the past to include into the plot
fig, ax = plt.subplots(num='Forecast on Index of '+name, figsize=(9.7, 6))
data.iloc[-past_pad:].plot(ax=ax)
init_val = data[-past-1]
fc_frame = pd.DataFrame(data=np.zeros((horizon, n_sims)), index=fc_index)
fc_frame.loc[start_d] = init_val
fc_frame.sort_index(inplace=True)
for sim in range(n_sims):
fc_frame.loc[1:, sim] = init_val * (forecast.simulations.values[-past-1][sim] / 100 + 1).cumprod()
fc_frame[sim].plot(ax=ax, color='#55a868', alpha=0.01, linewidth=5)
confintvls = [fc_frame.iloc[0, :2].tolist()]
for day in range(1, horizon+1):
confintvls.append(np.percentile(fc_frame.iloc[day], (2.5, 97.5)).tolist())
ci_frame = pd.DataFrame(confintvls, index=fc_frame.index)
ci_frame.plot(ax=ax, color='#c44e52', legend=False)
ax.set_ylabel('Time')
ax.set_ylabel('Index')
fig.tight_layout()
In [20]:
max_months = 1000 # Number of business months to run the simulation for
init_val = 100. # Value that simulated index starts from
initial_vals = returns[-max(har.lags)]
initial_vola = result.conditional_volatility[-max(harch.lags)]
sim_returns = model.simulate(result.params,
max_months,
burn=1000,
initial_value=initial_vals,
initial_value_vol=initial_vola)
sim_trace = [init_val]
sim_trace.extend(init_val * (sim_returns.data / 100 + 1).cumprod())
fig = plt.figure(num='Simulated Traces for the '+name+' Index and its Volatility', figsize=(9.7, 6))
vals_ax = fig.add_subplot(211)
vola_ax = fig.add_subplot(212)
vals_ax.plot(sim_trace)
vals_ax.set_xlabel('Months')
vals_ax.set_ylabel('Index')
vals_ax.set_xlim(left=0, right=max_months)
sim_returns.volatility.plot(ax=vola_ax, color='#c44e52')
vola_ax.set_xlabel('Months')
vola_ax.set_ylabel('Volatility')
vola_ax.set_xlim(left=0, right=max_months)
fig.tight_layout()
In [ ]: