前段时间整理了一个预测的基本思考框架和常见的方法，其中提到了ARIMA模型，在《大数据预测》那本书里，ARIMA是单独开辟一章讲的，比较复杂和难理解的一个模型，自己最近找了点资料粗浅学习了一下理论，并尝试用Python实现一下。
一、时间序列数据及其预处理
1.时序数据
时序数据顾名思义就是随着时间而变动的数据，是指某个个体在不同时间点上收集到的数据。已经被收集(或者叫观察)到的数据其实是时间序列变量的一个观察值，但由于时间的不可逆性，每一个时间点的变量有且仅能有一个观察值，我们用这些观察值拟合预测模型，用来预测未来时刻时间序列变量的值(此时，未发生的时间序列变量是一个随机变量)。
2.时序数据预处理
时序数据的预处理主要是平稳性和随机性检验，根据检验结果的不同，用不同的模型对数据进行建模。 纯随机序列，又称为白噪声序列，序列的变量之间无任何关系，因此没有分析价值；平稳非白噪声序列，其均值和方差是常数，目前有比较成熟的模型对其进行模拟，主要是AR(Autoregressive model)自回归模型，MA(Movingaverage model)移动平均模型、以及ARMA-一个AR和MA结合的模型。
(1)平稳性检验
平稳性检验的基本原理是检验不同时间点时序变量之间的相关关系，由于是一个时间序列数据的不同时刻之间的相互比较，因此取名为自协方差函数和自相关系数，但基本原理和协方差函数、相关系数类似，当一个时间序列有为常数的均值和方差，且其自协方差函数只跟时间间隔有关，即γ(t,s) =γ(k,k+s-t) ，其中γ是自协方差函数。 平稳性检验的方法分为两类，一类是偏主观和模糊的图示检验，一类是比较精确的数学检验。 第一类主要是用时序图和自相关图进行观察，如果一个时间序列的时序图围绕某一个常数上下波动，且波动范围有界，可以认为是平稳的，除此之外，如果自相关图中自相关系数比较快速衰减到0，并且在零附近震荡，可以认为序列是平稳序列，这是因为平稳序列通常具有短期相关性，随着延迟期数的增加，变量之间的相关关系会指数降低。
第二类数据的检验方法是单位根检验。单位根又称为单位圆，指模型的特征方程的根和1的关系，如果模型的特征根小于1，则序列非平稳，如果特征根大于1，则序列平稳。
(2)白噪声检验
白噪声序列的各个序列值之间没有任何关系，即γ(k)=0，在实际中，只要序列变量之间的协方差系数在0值附近波动，就可以认为是白噪声序列了。
二、常用模型
平稳时间序列的建模通常是AR、MA和ARMA模型，它们都是一种多元线性回归模型。
1.AR模型
这个模型的定义是t时刻的序列取值和前p期的序列值有关，且是前p期的线性组合。这个比较容易理解，在通常的印象中，我们大都觉得这次发生的事和最近几次发生的事有关联。 并非所有的AR模型都是平稳的，AR模型满足平稳性的条件是系数|φ|<1； 同时自相关系数ACF呈指数的速度衰减，在0附近波动，表现出拖尾性；其次，由于AR模型里在计算t和t-k期之间的相关性时会受到中间k期的影响，剔除k期影响后的自相关系数称为偏自相关系数PACF，平稳AR模型的PACF具有截尾性(尾巴突然断掉)，这是因为yt只跟y(t-1)和y(t-p)之间的序列值有关，当时间间隔k>p的时候，yt和y(t-k)之间就没有直接关系了。
2.MA模型
这个模型刚开始不太好理解，某个时刻的序列值跟过去历史q期的白噪声有关？are you kidding me？后来看网上有人说MA模型关注的是AR模型中的随机扰动项的累加，通过扰动项的移动平均能够有效消除预测中的随机波动。 MA模型的自相关系数在q阶后全部等于0，具有截尾性，偏自相关系数在q阶后拖尾，这是因为任何一个可逆的MA模型都是可以转换为无穷阶的AR模型(各种嵌套回去)，而无穷阶的AR模型的PACF是拖尾的。
3.ARMA模型
ARMA模型是AR和MA模型的线性组合，特别的，当q=0，它是一个AR(p)模型，如果p=0，它是一个MA(q)模型。 平稳条件下的三类模型的自相关系数ACF和偏自相关系数PACF特性如下：
利用这些特性可以通过图示的方法判断平稳性。
三、建模步骤
时间序列经过数据预处理被判定为平稳非白噪声序列后，就可以利用ARMA进行建模。 1.计算ACF、PACF； 2.ARMA模型定阶，即确定最优p、q值； 3.估计模型参数、进行参数检验； 4.模型检验； 5.模型预测及评价 在这里提及一下ARIMA模型，ARIMA是指时间序列经过I阶差分之后可以满足平稳非白噪声序列的条件，并利用ARMA为其进行建模。因此，ARIMA和ARMA最大的差异就是前者对原数据序列进行了差分，差分次数一般为1次，2次都不太推荐，因为2次差分之后很可能就成为白噪声序列，别问我怎么知道的(╯︵╰)。
四、Python试一试
会列出一些重点步骤和代码。 1.把数据分成了data_time的模型内数据和用于判断模型效果的data_time_test模型外数据
data_time_origin = pd.read_excel(r'C:\*\arima_data_gmv.xlsx',index_col = '日期')data_time = data_time_origin[:'2017-05-20'] #将5-20号之前的数据作为模型数据data_time_test = data_time_origin['2017-05-21':]
2.看一下时序图和自相关函数图判断是否平稳序列
data_time.plot(figsize=(20,8),fontsize=15)#导入相应的统计模型包from statsmodels.graphics.tsaplots import plot_acf #acf计算from statsmodels.tsa.stattools import adfuller as ADF #单位根检测plot_acf(data_time).show()
时序图上能看到还是有比较明显的趋势的，初步判断不是平稳序列。
ACF上看自相关系数长期大于0，且后期又长期小于0，呈现一种对称模式，可以判断不是平稳序列。单位根检验也证实不是平稳序列。
adf_test = ADF(data_time['gmv'])print(adf_test)'''统计值=临界值不能拒绝零假设，存在单位根，序列非平稳。可以看到p-value=0.6377显著大于0.05，另外统计量计算值-1.28显著大于1%、5%、10%的临界值，因此序列存在单位根，序列是非平稳的'''result:(-1.281137399884022, 0.6377571909319921, 5, 75, {'1%': -3.520713130074074, '5%': -2.9009249540740742, '10%': -2.5877813777777776}, 573.5874067738569)
3.进行一阶差分运算看其是否能达到平稳
data_time2 = data_time.diff(periods=1)data_time2 = data_time2.dropna()data_time2.columns = ['gmv差分'] data_time2.plot(figsize=(20,8),fontsize = 15)
可以看到基本消除了趋势，数据有稳定的均值。
#自相关图plot_acf(data_time2).show()#偏自相关图from statsmodels.graphics.tsaplots import plot_pacfplot_pacf(data_time2).show() #单位根检验adf_test_d = ADF(data_time2['gmv差分'])print(adf_test_d)result:(-8.458458355369716, 1.591087857743347e-13, 4, 75, {'1%': -3.520713130074074, '5%': -2.9009249540740742, '10%': -2.5877813777777776}, 565.8248717064124)
自相关图能看出有一些短期的相关性，单位根检验的结果是统计量-8.45小于临界值，p-value也远远小于0.05，差分之后是平稳序列。
#白噪声检验from statsmodels.stats.diagnostic import acorr_ljungboxprint(acorr_ljungbox(data_time2, lags=[3]))result:(array([13.84640456]), array([0.00312186]))
白噪声检验的p-value 0.003小于0.05，因此拒绝原假设，差分后的时序不是白噪声序列。因此上面的检验过程说明差分后的时间序列是平稳非白噪声序列，可以进行平稳时序建模。
4.模型定阶
上面的PACF图，能看到诡异的拖尾性，前半段截尾，后半段拖尾...看着像是1阶AR模型，具体是不是、行不行我也不知道。 暴力用BIC信息值算一下，找最优的pq值。
from statsmodels.tsa.arima_model import ARIMAdata_time['gmv'] = data_time['gmv'].astype(float)#定阶pmax = int(len(data_time)/10) #一般阶数不超过length/10qmax = int(len(data_time)/10) #一般阶数不超过length/10bic_matrix = [] #bic矩阵for p in range(pmax + 1): temp = [] for q in range(qmax + 1): try: temp.append(ARIMA(data_time,(p,1,q)).fit().bic) except: temp.append(None) bic_matrix.append(temp)bic_matrix = pd.DataFrame(bic_matrix)bic_matrix = bic_matrix.stack()p,q = bic_matrix.idxmin()print('BIC值最小的p和q分别为：p=%s, q=%s' %(p,q))result:BIC值最小的p和q分别为：p=3, q=3
所以定阶结果竟然是3和3，那就试着用这个模拟一下看看。
5.构建最终模型
arima_model = ARIMA(data_time,(p,1,q)).fit() arima_model.summary()
强烈怀疑这个模型的有效性... wierd，谁能给指条明路啊。用了80天的数据做训练，然后在p=3，q=3的时候，模型结果只给了系数，没给有效性的信息，摊手。 6.既然模型参数没法检验，只能用预测结果看看效果
forecast_n = pd.DataFrame(arima_model.forecast(11)[0],index = pd.date_range('2017-05-21',periods = 11))forecast_n.columns = ['gmv_fore']forecast_nerror_test = pd.merge(data_time_test,forecast_n,left_index = True, right_index = True)error_test['误差'] = error_test['gmv'] - error_test['gmv_fore']plt.figure(figsize=(15,6))plt.plot(error_test.index,error_test['gmv'],'-*',error_test.index,error_test['gmv_fore'],'-go')plt.legend(loc = 'best', fontsize=20)plt.show()print('预测值与真实值之间的绝对误差为：', np.abs(sum(error_test['gmv']) - sum(error_test['gmv_fore'])))result:预测值与真实值之间的绝对误差为：17.937944428429546
其中蓝线是原值，绿线是预测值，看结果还行，把上升和下降的趋势都已经预测出来了，用原值的和同预测值的和做差取绝对值，绝对值17.9，在10天的预测区间内，200多的基础数量级上，这个误差我个人觉得还挺nice啊，比我之前做的那些劳什子玩意强多了，而且还是到天啊，到天！以前都是月度的。
7.做一点模型评价
计量经济学里进行模型评价的指标有三种：一是均方误差MSE、一个是绝对平均误差MAE，另外一个是符号正确百分率。符号正确百分率基本用于经济领域，或者说对增长和下降的趋势预测，只是方向对了就算预测不错了，所以此处用前两种。
'''1.均方误差2.绝对平均误差'''#均方根误差&均方误差rmse = np.sqrt(sum(error_test['误差']**2)/len(error_test['误差']))mse = sum(error_test['误差']**2)/len(error_test['误差'])print('均方根误差为：%.2f' % rmse,'\n','均方误差为：%.2f' % mse)#绝对平均误差mae = sum(abs(error_test['误差']))/len(error_test['误差'])print('绝对平均误差为：%.2f' % mae)result:均方根误差为：18.69 均方误差为：349.30绝对平均误差为：13.10
当然这几个指标是需要在不同模型之间比较才能有优劣关系的，也就是说选用不同的p，q，或者甚至是非ARIMA模型的其它模型预测出来的结果之间的相互比较，但鉴于我今天已经不想搞了，所以就不比较了。但无论怎样，我觉得这个预测的结果还真是挺不错的，对于天的预测竟然趋势全对，而且在200的基数量级上最多是5-30日那天预测误差49，预测多了49万，其它时候模型的效果还是蛮好的，我很满意，嗯，就是这样！
参考资料： 1.中国大学MOOC-中南财经政法大学《时间序列分析》. 2.中国大学MOOC-对外经济贸易大学《计量经济学导论》. 3.张良均 等. Python数据分析与挖掘实战[M].北京:机械工业出版社,2015:119-134. 4.网上众多牛人文章.