时序分析 29 - 时序预测 - 格兰杰因果关系（下） python实践2
Posted Magic Ktwc37
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了时序分析 29 - 时序预测 - 格兰杰因果关系（下） python实践2相关的知识，希望对你有一定的参考价值。
时序分析 29 Granger-Cause 实践2
金融时序数据分析
本文主要搭建经济模型，从而分析和预测金融时序数据，重点关注黄金和原油价格。这两种金融资产对全球经济影响较大，我们将使用格兰杰因果检验来测试多个变量之间的关系，然后会使用VAR(Vector Auto Regression)通过历史数据来尝试预测黄金和原油的未来价格。我们使用的数据包括黄金、白银、原油、股票指数、利率和美元。
流程总结
下图大致概括了上一篇博文和本篇博文两次时间的基本流程。
探索式分析
dataset = dataset.loc['20000101':'20180501']
dataset.head()
dataset.fillna(method='pad',inplace=True)
dataset.dropna(inplace=True)
画出时序图
# Plot
fig, axes = plt.subplots(nrows=3, ncols=2, dpi=120, figsize=(10,6))
for i, ax in enumerate(axes.flatten()):
data = dataset[dataset.columns[i]]
ax.plot(data, color=’red’, linewidth=1)
ax.set_title(dataset.columns[i])
ax.xaxis.set_ticks_position(‘none’)
ax.yaxis.set_ticks_position(‘none’)
ax.spines[“top”].set_alpha(0)
ax.tick_params(labelsize=6)
plt.tight_layout();
观察上面的时序图可以看出，每个时序都含有单位根，也就是说是不平稳的。
为了进一步了解数据，我们对其进行正态分布检验
stat,p = stats.normaltest(dataset.Gold)
print("Statistics = %.3f, p=%.3f" % (stat,p))
alpha = 0.05
if p > alpha:
print('Data looks Gaussian (fail to reject null hypothesis)')
else:
print('Data looks non-Gaussian (reject null hypothesis)')
output: Statistics = 658.293, p=0.000 Data looks non-Gaussian (reject null hypothesis)
拒绝黄金价格数据为正态分布的假设。
峰度小于0说明数据具有长尾特性
plt.figure(figsize=(14,6))
plt.subplot(1,2,1)
dataset['Gold'].hist(bins=50)
plt.title('Gold')
plt.subplot(1,2,2)
stats.probplot(dataset['Gold'], plot=plt);
dataset.Gold.describe().T
从图上看也可以推断出黄金数据的分布离正态分布相距甚远。
看一下数据中各变量之间的相关性
corr = dataset.corr()
sns.heatmap(corr,xticklabels=corr.columns.values,
yticklabels=corr.columns.values, annot=True,
annot_kws='size':12)
heat_map = plt.gcf()
heat_map.set_size_inches(10,6)
plt.xticks(fontsize=10)
plt.yticks(fontsize=10)
plt.show()
划分测试数据和训练数据
n_obs=15
X_train, X_test = dataset[0:-n_obs], dataset[-n_obs:]
print(X_train.shape, X_test.shape)
(5114, 6) (15, 6)
平稳化转换
做一阶差分把数据转换为平稳序列。
transform_data = X_train.diff().dropna()
平稳性检验
对差分后的数据应用ADF检验其平稳性
def augmented_dickey_fuller_statistics(time_series):
result = adfuller(time_series.values)
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('Critical Values:')
for key, value in result[4].items():
print('\\t%s: %.3f' % (key, value))print('Augmented Dickey-Fuller Test: Gold Price Time Series')
augmented_dickey_fuller_statistics(X_train_transformed['Gold'])
print('Augmented Dickey-Fuller Test: Silver Price Time Series')
augmented_dickey_fuller_statistics(X_train_transformed['Silver'])
print('Augmented Dickey-Fuller Test: Oil Price Time Series')
augmented_dickey_fuller_statistics(X_train_transformed['Oil'])
print('Augmented Dickey-Fuller Test: Interest_rate Time Series')
augmented_dickey_fuller_statistics(X_train_transformed['Interest_rate'])
print('Augmented Dickey-Fuller Test: Stock_Index Time Series')
augmented_dickey_fuller_statistics(X_train_transformed['Stock_Index'])
print('Augmented Dickey-Fuller Test: USD_Index Time Series')
augmented_dickey_fuller_statistics(X_train_transformed['USD_Index'])
fig, axes = plt.subplots(nrows=3, ncols=2, dpi=120, figsize=(10,6))
for i, ax in enumerate(axes.flatten()):
d = X_train_transformed[X_train_transformed.columns[i]]
ax.plot(d, color='red', linewidth=1)
# Decorations
ax.set_title(dataset.columns[i])
ax.xaxis.set_ticks_position('none')
ax.yaxis.set_ticks_position('none')
ax.spines['top'].set_alpha(0)
ax.tick_params(labelsize=6)
plt.tight_layout();
格兰杰因果测试
检查各变量之间是否存在格兰杰因果关系，下面所使用的函数来自上一篇博文。
grangers_causality_matrix(dataset, variables = dataset.columns)
如上一篇博客所言，单元中的值表示两个变脸的格兰杰因果关系显著性检验的P值，如果该值<0.05，则表示我们在统计意义上接受列变量对行变量具备格兰杰因果关系。例如上表中第一行第二列就表示白银价格对黄金价格具备格兰杰因果关系，而第二行第五列代表股票指数对白银价格不具备格兰杰因果关系。
我们可以看到上表中的变量之间有很多存在格兰杰因果关系，这就使VAR模型的正确性成为可能。
VAR模型
mod = smt.VAR(X_train_transformed)
res = mod.fit(maxlags=15, ic='aic')
print(res.summary())
残差可能有相关性的，我们画图看一下:
y_fitted = res.fittedvalues
plt.figure(figsize = (15,5))
plt.plot(residuals, label='resid')
plt.plot(y_fitted, label='VAR prediction')
plt.xlabel('Date')
plt.xticks(rotation=45)
plt.ylabel('Residuals')
plt.grid(True)
情况不严重，可以暂时忽略。
Durbin-Watson统计量
采用Durbin-Watson统计量检查序列自相关性，此统计量是一个0到4的数值。值为2说明序列中没有自相关现象存在，02意味着很可能存在正自相关而2到4之间说明存在负自相关。一般情况下，我们认为该值在1.52.5之间都属正常，而如果超出这个范围就说明模型欠拟合，应该尝试调参、检查数据或者更换模型。
from statsmodel.stats.stattools import durbin_watson
output = durbin_watson(res.resid)
for col,val in zip(transform_data.columns,output):
print(col,':',round(val,2))
检验结果表明残差之间不存在自相关性。
应用VAR模型进行预测
# Get the lag order
lag_order = res.k_ar
print(lag_order)
# Input data for forecasting
input_data = X_train_transformed.values[-lag_order:]
print(input_data)
# forecasting
pred = res.forecast(y=input_data, steps=n_obs)
pred = (pd.DataFrame(pred, index=X_test.index, columns=X_test.columns + '_pred'))
print(pred)
反转差分转换
记得我们为了保证时序的平稳性做过了一阶差分，现在预测结果的时候我们需要将其转换回来。
# inverting transformation
def invert_transformation(X_train, pred):
forecast = pred.copy()
columns = X_train.columns
for col in columns:
forecast[str(col)+'_pred'] = X_train[col].iloc[-1] + forecast[str(col)+'_pred'].cumsum()
return forecast
output = invert_transformation(X_train, pred)
#combining predicted and real data set
combine = pd.concat([output['Gold_pred'], X_test['Gold']], axis=1)
combine['accuracy'] = round(combine.apply(lambda row: row.Gold_pred /row.Gold *100, axis = 1),2)
combine['accuracy'] = pd.Series(["0:.2f%".format(val) for val in combine['accuracy']],index = combine.index)
combine = combine.round(decimals=2)
combine = combine.reset_index()
combine = combine.sort_values(by='Date', ascending=False)
评估
#Forecast bias
forecast_errors = [combine['Gold'][i]- combine['Gold_pred'][i] for i in range(len(combine['Gold']))]
bias = sum(forecast_errors) * 1.0/len(combine['Gold'])
print('Bias: %f' % bias)
print('Mean absolute error:', mean_absolute_error(combine['Gold'].values, combine['Gold_pred'].values))
print('Mean squared error:', mean_squared_error(combine['Gold'].values, combine['Gold_pred'].values))
print('Root mean squared error:', sqrt(mean_squared_error(combine['Gold'].values, combine['Gold_pred'].values)))
