Python爬取中国银行外汇牌价(statsmodels预测分析)--(二)

• Author: sunhailin-Leo
• E-mail:379978424@qq.com
• 代码: 仓库地址 -> github.com/sunhailin-L…

简介

• 第一篇文章传送门: Python爬取中国银行外汇牌价(爬虫 + PyFlux简单预测分析)--(一)

• 如果数据还没获取到的话请移步看第一篇文章(本篇文章默认数据已存储在数据库中)

步入正题

• 上一篇文章末尾部分在讲使用Pyflux这个库对数据进行分析预测
• 总结两个方面:
  • 优点:
    • Pyflux模型文档"一针见血"(建立在对时序分析有一定基础的人, 能看懂部分核心公式)
  • 缺点:
    • 提供少量的数据分析API, 不像statsmodels提供了例如残差分析等方法进行模型验证调优的方法
• 本文将使用statsmodels对此前的数据进行分析。

数据前期准备

df['查询时间'] = df['查询时间'].apply(lambda x: x[:-9])
df['查询时间'] = pd.to_datetime(df['查询时间'], format="%Y-%m-%d")
df = df.groupby('查询时间')['现汇卖出价'].mean()
df = df.to_frame()
print(df)
复制代码

数据平稳性校验

• 直接上代码(差分画图看数据平稳性)
• 一共测试了1阶和2阶(不建议使用高阶数据, 容易数据造成破坏)

# 差分图
fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(111)
# 里面的1代表了差分阶数
diff1 = df.diff(1)
diff1.plot(ax=ax1)
plt.show()
# fig.savefig('./picture/1.jpg')
复制代码

• 一阶差分结果

• 二阶差分结果

• 结论：由观察可得出一阶差分波动相对较小，较为平稳 ( 之后建模的数据使用1阶差分的结果进行建模 )

自相关图和偏自相关图

• 介绍一些自相关性和非自相关性:

  • 自相关性(ACF): 自相关性是指随机误差项的各期望值之间存在着相关关系，称随机误差项之间存在自相关性（autocorrelation）或序列相关
  • 偏自相关性(PACF): 偏自相关是剔除干扰后时间序列观察与先前时间步长时间序列观察之间关系的总结（Partial autocorrelation）

• 接下来选择40个滞后做自相关和偏自相关的分析

# statmodels 自相关图 偏相关图
diff1 = diff1.dropna()
fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(diff1, lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(diff1, lags=40, ax=ax2)
fig.show()
# fig.savefig('./picture/acf_pacf.jpg')
复制代码

• 接下来由看图观察定义ARMA模型的p(自相关系数), q(偏自相关系数)参数:

  • 自相关图显示滞后有三个阶超出了置信边界
  • 偏相关图显示在滞后2阶时的偏自相关系数超出了置信边界(其实波动还是没有完全消除, 只是2阶的时候比较突出)

• 可以尝试如下p, q参数的模型, 并输出AIC,BIC,HQIC的结果

  • AIC:赤池信息量 akaike information criterion --> -2ln(L) + 2k
  • BIC:贝叶斯信息量 bayesian information criterion --> -2ln(L) + ln(n) * k
  • HQIC:HQ准则 --> -2ln(L) + ln(ln(n)) * k
  • L是在该模型下的最大似然，n是数据数量，k是模型的变量个数。

• 尝试结果如下:

  • p=2, q=0 --> -194.77200890459767 -179.81283725588074 -188.79262385129292
  • p=2, q=1 --> -197.01722373566554 -178.31825917476937 -189.54299241903462
  • p=2, q=2 --> -195.11076756537022 -172.6720100922948 -186.1416899854131
  • p=3, q=2 --> -201.37730533090803 -175.19875494565338 -190.91338148762475

• 模型选择:

  • 在上述的结果看来选择p=2, q=2为最优办法, 虽然不一定是最适合, 但是在差分和根据自相关和偏自相关图的结果看来, 可以尝试使用(2, 2)参数进行建模

建模

• 建模的代码很短,就不做过多叙述,看下述代码

arma_mod22 = sm.tsa.ARMA(diff1, (3, 2)).fit()
# 输出AIC,BIC和HQ准则结果
print(arma_mod22.aic, arma_mod22.bic, arma_mod22.hqic)
复制代码

模型验证

1. 残差验证

# 残差(输出形式为DataFrame)
resid = arma_mod22.resid
# 残差的ACF和PACF图
fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(resid.values.squeeze(), lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(resid, lags=40, ax=ax2)
fig.show()
# fig.savefig('./picture/resid_acf_pacf.jpg')
复制代码

• 根据上图可以发现虽然少数部分的阶数下还有出现了超过置信区间的问题，但是总体看来序列残差基本为白噪声。

2. D-W验证

• D-W验证是检验自相关性的方法(只适用于检测1阶自相关性, 高阶需要另辟蹊径)

# 残差D-W检验
resid_dw_result = sm.stats.durbin_watson(arma_mod22.resid.values)
# 1.9933441709003574 接近于 2 所以残差序列不存在自相关性。
print(resid_dw_result)
复制代码

• 结果等于1.9933441709003574
• 结论：残差序列不存在自相关性
• 理由：DW结果接近2的时候序列几乎不存在自相关性(参考文章: www.doc88.com/p-543541848… )

3. 正态分布

# 正态校验 -> 基本符合正态分布
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111)
fig = sm.qqplot(resid, line='q', ax=ax, fit=True)
fig.show()
# fig.savefig('./picture/normal_distribution.jpg')
复制代码

• 除了右上角的部分异常点(这几个点的数据很直观的在接下来的Q检验中显现出来)
• 由图得出，绝大部分数据还是符合正态分布

4. Q检验

# 残差序列Ljung-Box检验，也叫Q检验
r, q, p = sm.tsa.acf(resid.values.squeeze(), qstat=True)
data = np.c_[range(1, 41), r[1:], q, p]
table = pd.DataFrame(data, columns=['lag', "AC", "Q", "Prob(>Q)"])
temp_df = table.set_index('lag')
print(temp_df)
# Prob(>Q)的最小值: 0.025734615668093132, 最大值: 0.9874705305611844, 均值: 0.2782013984159408
prob_q_min = temp_df['Prob(>Q)'].min()
prob_q_max = temp_df['Prob(>Q)'].max()
prob_q_mean = temp_df['Prob(>Q)'].mean()
print("Prob(>Q)的最小值: {}, 最大值: {}, 均值: {}".format(prob_q_min, prob_q_max, prob_q_mean))
复制代码

• 结果：Prob(>Q)的最小值: 0.025734615668093132, 最大值: 0.9874705305611844, 均值: 0.2782013984159408
• 由上图可以看到最后10个数据的值几乎在0.05左右，可以对应回正态分布的那些离散点的值
• 根据Q校验的方法, 在95%的置信区间内, 当Prob大于0.05当前残差序列不存在自相关性, 因此滞后40个数据的模型属于基本不存在自相关性。

5. 预测(预测从11月9日到11月14日的汇率变动)

• 由于此前将数据进行了1阶差分，因此预测结果的数据也是1阶差分的预测值

# 预测
predict_data = arma_mod22.predict('2018-11-09', '2018-11-14', dynamic=True)
# 画预测图
fig, ax = plt.subplots(figsize=(12, 8))
ax = diff1.ix['2018-01-01':].plot(ax=ax)
fig = arma_mod22.plot_predict('2018-11-09', '2018-11-14', dynamic=True, ax=ax, plot_insample=False)
fig.show()
# 结果预测
last_data = df['现汇卖出价'].values[-1]
# 还原结果
predict_data_list = predict_data.values.tolist()
restore_list = []
for d in predict_data_list:
    last_data = last_data + d
    restore_list.append(last_data)
predict_data = pd.DataFrame(restore_list, index=predict_data.index, columns=['现汇卖出价'])
print(predict_data)
复制代码

• 注:以上预测数据均为当天的汇率的均值

总结

• 以上的分析、建模、校验、预测的流程只是简单的一个流程.预测模型准确性还需要不断验证, 以上的关键代码仅提供了一套可以参考的流程.
• 时序数据预测与分析还有待学习，针对不同类型的数据和多变量的数据有着不同的模型, 以上步骤仅供参考.
• 目前使用过时序预测的库有：statsmodels和pyflux, 之后考虑使用机器学习或者深度学习的方法进行预测.