【Python小竞赛】ARIMA算法预测三日后招商银行收盘价

介绍

本文整理记录了参与的一次小型数据分析竞赛【数据游戏】，竞争目标是预测2019年5月15日A股闭市时招商银行600036的股价。

主要思路是利用ARIMA算法做时间序列预测。

使用的数据是公开的数据集 tushare。

拿到题目和数据之后，首先结合既往经历，觉得想要预测准股价，本身是一个不可能的事情，尤其是A股。

影响股价的因素非常复杂而且不透明，以及金融投资领域具有的反身性理论，使得这次预测更多偏向于实验性质，同时对竞争结果不要有过高的期望。

预测得准，是你的运气；预测的烂，也不会影响你从中学到什么。学习第一，比赛第二吧。

鉴于以上，本次预测只使用了close的时间序列。更多的数据其实并没有什么用。

首先了解下本文的 ARIMA 建模过程：

获取时间序列数据；

观察数据是否为平稳时间序列；对于非平稳时间序列要先进行d阶差分运算，转化为平稳时间序列；

对处理得到的平稳时间序列，求它的阶数p, q

根据ARIMA算法建模，ARIMA(data, order=(p, d, q))

模型检验和调优

预测

初学的小伙伴们可能对这个过程并不熟悉，没关系，先背下来。

对Python不熟悉的小伙伴们，我给大家总结了几句车轱辘话，大家先强行了解下：

获取时间序列数据；
data = pd.read_excel(‘600036.xlsx’, index=None)
train = data[‘close’]

观察数据是否为平稳时间序列；对于非平稳时间序列要先进行d阶差分运算，转化为平稳时间序列；
adf_data = sts.adfuller(train)
diff = train.diff(1)

对处理得到的平稳时间序列，求它的阶数p, q
sm.tsa.arma_order_select_ic(train, max_ar=8, max_ma=8, ic=[‘aic’, ‘bic’, ‘hqic’])

根据ARIMA算法建模
ARMAModel = sm.tsa.ARIMA(train, order=(4,1,2)).fit()

模型检验和调优
train_shift = train.shift(1)
pred_recover = predicts.add(train_shift)
np.sqrt( sum( (pred_recover -train) ** 2)/train.size )

预测
f = ARMAModel.forecast(3)

以上，是本文的核心代码，大家如果一时看不懂，可以跳过。
可以看看下面更详细的步骤。

导入数据并处理

# 导入必须的模块
import tushare as ts  #使用的公开的数据
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
import statsmodels.tsa.stattools as stsimport wconfig as wc  #自定义打印输出控制模块，与本次竞赛无关
wc.dispy() # 导入数据，从去年1月开始即可
data = ts.get_hist_data('600036', start='2018-01-08', end='2019-05-10').sort_index()
#data = ts.get_hist_data('600036', start='2018-04-01').sort_index().reset_index()# 数据安全和源数据备份
#data.to_excel('600036.xlsx')# 源数据可靠性检查。12日发现tushare的数据有错误，需要手工矫正（13日该数据恢复正常）。
print('tushare中5月10日close股价为: %.2f元，与实际不符！' % (data['close']['2019-05-10']))
data['close']['2019-05-10']=33.61
print('当日实际收盘价应为: %.2f元' % data['close']['2019-05-10'])# 只取 close 字段作为训练数据
train = data['close']
train.index = pd.to_datetime(train.index)  # 将字符串索引转换成时间索引
train.tail()
train.tail().index

tushare中5月10日close股价为: 33.48元，与实际不符！
当日实际收盘价应为: 33.61元
DatetimeIndex([‘2019-05-06’, ‘2019-05-07’, ‘2019-05-08’, ‘2019-05-09’, ‘2019-05-10’],
dtype=‘datetime64[ns]’, name=‘date’, freq=None)

检验时间序列的稳定性

ARIMA算法要求时间序列稳定，所以在建模之前，要先检验时间序列的稳定性。
adfuller就是用来干这个的。
adfuller 全称 Augmented Dickey–Fuller test，即扩展迪基-福勒检验，用来测试平稳性。


#先做一个解释器, 让 adfuller 的输出结果更易读易理解：def tagADF(t):result = pd.DataFrame(index=["Test Statistic Value", "p-value", "Lags Used", "Number of Observations Used", "Critical Value(1%)", "Critical Value(5%)", "Critical Value(10%)"],columns=['value'])result['value']['Test Statistic Value']=t[0]result['value']['p-value']=t[1]result['value']['Lags Used']=t[2]result['value']['Number of Observations Used'] = t[3]result['value']['Critical Value(1%)']=t[4]['1%']result['value']['Critical Value(5%)']=t[4]['5%']result['value']['Critical Value(10%)']=t[4]['10%']print('t is:', t)return result

adfuller检验是检查时间序列平稳性的统计测试之一。这里的零假设是：序列 train 是非平稳的。
测试结果包括测试统计和差异置信水平的一些关键值。如果测试结果中的 P-value 小于临界值，我们可以拒绝原假设并说该序列是平稳的。

adf_data = sts.adfuller(train)
tagADF(adf_data)

检验结果显示，p-value=0.414, 远远大于5%的临界值，说明零假设是成立的，即序列 train 是非平稳的。

t is: (-1.7325346908056185, 0.4144323576685054, 0, 322, {‘1%’: -3.4508226600665037, ‘5%’: -2.870558121868621, ‘10%’: -2.571574731684734}, 523.9067372199033)

	value
Test Statistic Value	-1.733
p-value	0.414
Lags Used	0
Number of Observations Used	322
Critical Value(1%)	-3.451
Critical Value(5%)	-2.871
Critical Value(10%)	-2.572

时间序列平稳化

为了让时间序列平稳，需要对 train 序列做差分运算：

# df.diff 差分运算，默认是后一行减前一行
# http://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.diff.html?highlight=diff#pandas.DataFrame.diff
diff = train.diff(1).dropna()  # df.dropna 删除有缺失值的行
plt.figure(figsize=(11,6))  # 指定显示大小
plt.plot(diff, label='Diff')  # 绘制趋势图
plt.legend(loc=0)  # 显示图例，loc指定图例位置，0为最佳位置。

关于时间序列稳定性的判断标准，可参考这篇博客：https://blog.csdn.net/u012735708/article/details/82460962

检验差分数据的平稳性

检验差分后数据的平稳性，和第一次验证方法相同

adf_Data1 = sts.adfuller(diff)
tagADF(adf_Data1)  # p-value很小，零假设不成立，因此，diff数据序列符合平稳性要求。

确定ARIMA的阶数p,q

在这里，先简单解释下自己对 ARIMA 算法的理解。不正确的地方请大家指点。

ARIMA算法认为时间序列上 t 时刻的值由2部分构成，第一部分是由之前p项历史值决定的，比如15日的收盘价是12日，13日，14日的收盘价的线性回归，用AR§表示。

但是这个线性回归的输出值肯定和15日的实际收盘价有一个误差，假设是e。

所以第二部分就是如果表示这个 e 。

这个e 可以认为是之前q项误差的线性回归，用MA(q)表示。

两部分合起来就是ARMA。

然后为了让时间序列平稳，再加个d 阶差分操作功能，就变成了ARIMA算法。

所以在使用 ARIMA 算法之前，需要先确定（p, d, q）的值。

确定ARIMA的阶数p,q# ARMA(p,q)是AR(p)和MA(q)模型的组合，关于p和q的选择，一种方法是观察自相关图ACF和偏相关图PACF,
# 另一种方法是通过借助AIC、BIC等统计量自动确定。
ic = sm.tsa.arma_order_select_ic(train, max_ar=8, max_ma=8, ic=['aic', 'bic', 'hqic']
)
ic

{‘aic’: 0 1 2 3 4 5 6 7 8
0 1515.258 1156.947 946.795 843.736 760.826 690.958 682.101 640.248 633.681
1 560.380 561.763 563.331 563.563 564.110 561.817 562.557 559.347 558.893
2 561.718 561.732 565.059 565.206 563.484 563.309 558.487 557.715 559.437
3 563.278 565.709 563.834 562.475 562.920 556.913 556.418 559.442 558.503
4 564.659 564.667 562.689 565.922 567.899 563.612 558.405 nan 559.008
5 564.854 563.982 562.661 567.897 565.930 561.758 560.349 555.892 554.069
6 561.132 562.839 563.492 565.472 561.965 563.734 561.397 554.178 563.975
7 562.364 558.978 555.652 557.528 559.511 560.338 561.747 562.075 566.033
8 560.181 559.217 557.546 556.357 557.818 555.207 560.736 562.675 564.121,
‘bic’: 0 1 2 3 4 5 6 7 8
0 1522.813 1168.280 961.906 862.624 783.492 717.402 712.322 674.246 671.457
1 571.712 576.874 582.219 586.228 590.553 592.038 596.556 597.124 600.447
2 576.829 580.620 587.725 591.649 593.705 597.308 596.264 599.269 604.769
3 582.166 588.375 590.278 592.696 596.919 594.690 597.972 604.774 607.612
4 587.325 591.111 592.910 599.921 605.676 605.166 603.737 nan 611.895
5 591.298 594.203 596.660 605.673 607.484 607.090 609.459 608.780 610.734
6 591.353 596.838 601.269 607.027 607.297 612.843 614.284 610.843 624.418
7 596.363 596.755 597.206 602.860 608.620 613.225 618.412 622.518 630.253
8 597.958 600.771 602.878 605.466 610.705 611.872 621.179 626.895 632.119,
‘hqic’: 0 1 2 3 4 5 6 7 8
0 1518.273 1161.471 952.827 851.276 769.874 701.514 694.165 653.819 648.761
1 564.903 567.795 570.871 572.610 574.666 573.881 576.129 574.427 575.481
2 567.750 569.272 574.107 575.762 575.548 576.881 573.567 574.303 577.533
3 570.818 574.757 574.390 574.539 576.492 571.993 573.006 577.538 578.106
4 573.707 575.223 574.753 579.494 582.979 580.200 576.501 nan 580.120
5 575.410 576.045 576.233 582.977 582.518 579.854 579.953 577.004 576.689
6 573.196 576.411 578.572 582.060 580.061 583.338 582.509 576.798 588.103
7 575.936 574.058 572.240 575.624 579.115 581.450 584.367 586.203 591.669
8 575.261 575.804 575.642 575.961 578.930 577.827 584.864 588.311 591.265,
‘aic_min_order’: (5, 8),
‘bic_min_order’: (1, 0),
‘hqic_min_order’: (1, 0)}

建立模型并拟合数据

注意ARIMA的参数中，输入数据应该是原始数据 train，ARIMA 会根据 d 的值，对原始数据做 d 阶差分运算。
d的含义是，输入序列需要先经过一个d阶的差分，变成一个平稳序列后才能进行数据拟合。

ARMAModel = sm.tsa.ARIMA(train, order=(4,1,2)).fit()  # order=(p,d,q)
# fittedvalues和diff对比
plt.figure(figsize=(11, 6))
plt.plot(diff, 'r', label='Orig')
plt.plot(ARMAModel.fittedvalues, 'g',label='ARMA Model')
plt.legend()

模型评估和调优

# 样本内预测
predicts = ARMAModel.predict()# 因为预测数据是根据差分值算的，所以要对它一阶差分还原
train_shift = train.shift(1)  # shift是指series往后平移1个时刻
pred_recover = predicts.add(train_shift).dropna()  #这里add是指两列相加，按index对齐# 模型评价指标 1：计算 score
delta = ARMAModel.fittedvalues - diff
score = 1 - delta.var()/train.var()
print('score:\n', score)# 模型评价指标 2：使用均方根误差（RMSE）来评估模型样本内拟合的好坏。
#利用该准则进行判别时，需要剔除“非预测”数据的影响。
train_vs = train[pred_recover.index]  # 过滤没有预测的记录
plt.figure(figsize=(11, 6))
train_vs.plot(label='Original')
pred_recover.plot(label='Predict')
plt.legend(loc='best')
plt.title('RMSE: %.4f'% np.sqrt(sum((pred_recover-train_vs)**2)/train_vs.size))
plt.show()

score:
0.9482183377736877

# 局部数据观察
train_t = train_vs.tail(15)
pred_t = pred_recover.tail(15)
plt.figure(figsize=(11, 6))
train_t.plot(label='Original')
pred_t.plot(label='Predict')
plt.legend(loc='best')
plt.title('RMSE: %.4f'% np.sqrt(sum((pred_t-train_t)**2)/train_t.size))
plt.show()

预测目标

使用 forecast 对样本外的时间序列进行预测。

关于foreast和predict的区别：
predict 可以对样本内和样本外的进行预测，结果是一样的。
举例说明：forecast（10），表示对未来10个点进行预测，但是可以用model.fittedvalues查看样本内点的拟合值；
而predict(start,end)里面的参数0表示样本内的第一个数，以此类推。
如果想要预测样本外的数，需要将start设置为len(data)+1,即数据长度+1，才表示预测样本外的第一个数字。
而 forecast函数，是对样本外的数据进行预测。
但是这两个函数的预测结果是一样的。
另外，需要提到的是，ARIMA算法一般只能预测一点点，越长越不准确，即便是简单的正弦函数也不能准确预测。

# 预测15日close股价，即10日之后的第三个交易日的收盘价
# 但是通过上面的局部数据观察发现，预测的数据趋势会延迟1个交易日，所以就取f[0][1]
f = ARMAModel.forecast(3)  # 样本外预测
print('5月15日close时的股价为：%.2f 元' % f[0][1])#----结束-----

5月15日close时的股价为：33.82 元

祝我好运吧，明天15号了，就要开奖了，哈哈。

顺便安利下这次参与的小竞赛【数据游戏】，是由个人发起、异步社区（邮电出版社）赞助的活动。

这个比赛不需要：

深厚的 AI 背景
多年比赛经验
入场费用

这个比赛需要的只是：

基础的 Python 技能
入门级的机器学习知识
激情与时间

总体来说，这个比赛氛围适中，没有太大的心理负担。

所以感觉比较适合打算初入的小伙伴，参与之后，这样的经验对以后的学习进阶和职业发展都有好处。

对比赛有兴趣的可加QQ群：759677734

本文参考资料

tushare的官方文档
http://tushare.org/trading.html
https://tushare.pro/document/2?doc_id=27
利用ARIMA进行时间序列数据分析（Python 很详细）
https://blog.csdn.net/u012735708/article/details/82460962
python时间序列预测股票走势
https://www.jianshu.com/p/d78ce35e42c0
用Python做股市数据分析（一）
http://python.jobbole.com/87813/
Python基于皮尔逊系数实现股票预测
https://www.cnblogs.com/yszd/p/10058475.html
Python中利用LSTM模型进行时间序列预测分析
https://www.cnblogs.com/arkenstone/p/5794063.html
Python时间序列实战之ARIMA股票预测
https://blog.csdn.net/qq_40707407/article/details/81938941
Python预测股票走势
https://blog.csdn.net/dipizhong7224/article/details/79252029
Python 股票价格预测（以中兴通讯为例）
https://blog.csdn.net/weixin_43181891/article/details/82627766

----------------------------------【帝魔微著】小组------------------------------------------