作者：天琼，「数据游戏」优胜队伍成员
介绍
本文整理记录了参与的一次小型数据分析竞赛「数据游戏」，竞赛目标是预测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 npimport pandas as pdimport matplotlib.pyplot as pltimport statsmodels.api as smimport 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.61print('当日实际收盘价应为: %.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元
检验时间序列的稳定性
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)
时间序列平稳化
为了让时间序列平稳，需要对 train 序列做差分运算：
# df.diff 差分运算，默认是后一行减前一行# http://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.diff.html?highlight=diff#pandas.DataFrame.diffdiff = train.diff(1).dropna() # df.dropna 删除有缺失值的行plt.figure(figsize=(11,6)) # 指定显示大小plt.plot(diff, label='Diff') # 绘制趋势图plt.legend(loc=0) # 显示图例，loc指定图例位置，0为最佳位置。复制代码
关于时间序列稳定性的判断标准，可参考这篇博客：
检验差分数据的平稳性
检验差分后数据的平稳性，和第一次验证方法相同
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复制代码
建立模型并拟合数据
注意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：计算 scoredelta = ARMAModel.fittedvalues - diffscore = 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()复制代码 # 局部数据观察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 元。
顺便安利下这次参与的小，是由个人发起、异步社区（邮电出版社）赞助的活动。
总体来说，这个比赛氛围适中，没有太大的心理负担。
所以感觉比较适合打算初入的小伙伴，参与之后，这样的经验对以后的学习进阶和职业发展都有好处。
Ad Time
学习更多数据科学知识请关注微信公众号：read_csv
参与数据科学活动请加 QQ 群：759677734