前言

接上一篇博客，用python完成代码编写。
【学习笔记】时间序列模型(ARIMA)

一、示例

• 已知一个上市公司一段时期的开盘价，最高价，最低价，收盘价等信息，要求建立模型，预测股价。
• 这里只需要股票的收盘价(close)，我们可以把数据提取出来，并划分为训练集和测试集
• 本题我们把1-3月份的数据作为训练集，4-6月份的数据作为测试集

二、代码实现----python

全部数据的平稳性检验

# 差分法importpandas aspdimportmatplotlib.pyplot asplt# 导入数据ChinaBank =pd.read_csv('ChinaBank.csv',index_col ='Date',parse_dates=['Date'])#.diff(1)做一个时间间隔ChinaBank['diff_1']=ChinaBank['Close'].diff(1)#1阶差分#对一阶差分数据在划分时间间隔ChinaBank['diff_2']=ChinaBank['diff_1'].diff(1)#2阶差分fig =plt.figure(figsize=(12,10))#原数据ax1 =fig.add_subplot(311)ax1.plot(ChinaBank['Close'])#1阶差分ax2 =fig.add_subplot(312)ax2.plot(ChinaBank['diff_1'])#2阶差分ax3 =fig.add_subplot(313)ax3.plot(ChinaBank['diff_2'])plt.show()

运行结果：

结果分析：

可以看出，一阶差分和二阶差分后，平稳性变好。

ADF检验

importmatplotlib.pyplot aspltimportpandas aspdimportstatsmodels.api assmfromstatsmodels.tsa.seasonal importseasonal_decomposefromstatsmodels.tsa.stattools importadfuller asADF# 计算原始序列、一阶差分序列、二阶差分序列的单位根检验结果ChinaBank['diff_1']=ChinaBank['diff_1'].fillna(0)ChinaBank['diff_2']=ChinaBank['diff_2'].fillna(0)timeseries_adf =ADF(ChinaBank['Close'].tolist())timeseries_diff1_adf =ADF(ChinaBank['diff_1'].tolist())timeseries_diff2_adf =ADF(ChinaBank['diff_2'].tolist())# 打印单位根检验结果print('timeseries_adf : ',timeseries_adf)print('timeseries_diff1_adf : ',timeseries_diff1_adf)print('timeseries_diff2_adf : ',timeseries_diff2_adf)

运行结果：

timeseries_adf :(np.float64(0.5279198084831831),np.float64(0.9856974415734416),9,335,{ '1%':np.float64(-3.4500219858626227),'5%':np.float64(-2.870206553997666),'10%':np.float64(-2.571387268879483)},np.float64(-734.0738716811488))timeseries_diff1_adf :(np.float64(-6.17718554497899),np.float64(6.587109239761689e-08),8,336,{ '1%':np.float64(-3.449962981927952),'5%':np.float64(-2.870180642420163),'10%':np.float64(-2.5713734527352607)},np.float64(-735.8436797171294))timeseries_diff2_adf :(np.float64(-9.202545123160359),np.float64(1.9841232339614826e-15),13,331,{ '1%':np.float64(-3.4502615951739393),'5%':np.float64(-2.8703117734117742),'10%':np.float64(-2.5714433728242714)},np.float64(-717.2833732193096))

结果分析：

$ADF$检验的结果共有五个参数：

• 第一个值：表示 Test Statistic，即 T 检验，表示 T 统计量，假设检验值
• 第二个值：p-value，即 p 值，表示 T 统计量对应的概率值
• 第三/四个值：Lags Used，即表示延迟和测试的次数
• 第五个参数{ ‘1%’: xxx, ‘5%’: xxx, ‘10%’: xxx}：不同程度拒绝原假设的统计值

如何确定该序列是否平稳呢？

• 1%、5%、10%不同程度拒绝原假设的统计值和 ADF 假设检验值比较，ADF 假设检验值同时小于1%、5%、10%即说明非常好地拒绝该假设
• P-value 是否非常接近 0

图检验法

1. 原始数据

fig =plt.figure(figsize=(12,7))ax1 =fig.add_subplot(211)fig =sm.graphics.tsa.plot_acf(ChinaBank['Close'],lags=20,ax=ax1)ax1.xaxis.set_ticks_position('bottom')# 设置坐标轴上的数字显示的位置，top:显示在顶部  bottom:显示在底部ax2 =fig.add_subplot(212)fig =sm.graphics.tsa.plot_pacf(ChinaBank['Close'],lags=20,ax=ax2)ax2.xaxis.set_ticks_position('bottom')plt.show()

运行结果：

结果分析：

ACF中，大部分的值没有落在置信区间内，所以不具有平稳性。

2. 一次差分

fig =plt.figure(figsize=(12,7))ax3 =fig.add_subplot(211)fig =sm.graphics.tsa.plot_acf(ChinaBank['diff_1'],lags=20,ax=ax3)ax3.xaxis.set_ticks_position('bottom')ax4 =fig.add_subplot(212)fig =sm.graphics.tsa.plot_pacf(ChinaBank['diff_1'],lags=20,ax=ax4)ax4.xaxis.set_ticks_position('bottom')

运行结果：

结果分析：

由图形可以看出，大部分的值都落在了置信区间内。

划分训练集

# 提取Close列ChinaBank.index =pd.to_datetime(ChinaBank.index)sub =ChinaBank.loc['2014-01':'2014-06','Close']sub.head()# 划分训练测试集train =sub.loc['2014-01':'2014-03']test =sub.loc['2014-04':'2014-06']

平稳性检验

ADF检验

# ADF检验fromstatsmodels.tsa.stattools importadfuller asADFtimeseries_adf =ADF(train.tolist())# 打印单位根检验结果print('timeseries_adf : ',timeseries_adf)

运行结果：

timeseries_adf :(np.float64(-2.902677813259885),np.float64(0.04503748919120268),0,61,{ '1%':np.float64(-3.542412746661615),'5%':np.float64(-2.910236235808284),'10%':np.float64(-2.5927445767266866)},np.float64(-227.51423828600673))

图检验法

importstatsmodels.api assm# 绘制fig =plt.figure(figsize=(12,7))ax1 =fig.add_subplot(211)fig =sm.graphics.tsa.plot_acf(train,lags=20,ax=ax1)ax1.xaxis.set_ticks_position('bottom')# 设置坐标轴上的数字显示的位置，top:显示在顶部  bottom:显示在底部ax2 =fig.add_subplot(212)fig =sm.graphics.tsa.plot_pacf(train,lags=20,ax=ax2)ax2.xaxis.set_ticks_position('bottom')plt.show()

运行结果：

确定 p，q

1. 相关函数法

由训练集的 ACF 和 PACF 图可以看出，ACF拖尾、PACF截尾，所以是AR(1)。

2. AIC、BIC准则

以 BIC 准则为例，确定 p,q 的取值范围为 [0,5]，通过循环网格搜索所有组合的 BIC 的值

#遍历，寻找适宜的参数importitertoolsimportnumpy asnpimportseaborn assns#确定pq的取值范围p_min =0d_min =0q_min =0p_max =5d_max =0q_max =5#Initialize a DataFrame to store the results,，以BIC准则results_bic =pd.DataFrame(index=['AR{ }'.format(i)fori inrange(p_min,p_max+1)],columns=['MA{ }'.format(i)fori inrange(q_min,q_max+1)])forp,d,q initertools.product(range(p_min,p_max+1),range(d_min,d_max+1),range(q_min,q_max+1)):ifp==0andd==0andq==0:results_bic.loc['AR{ }'.format(p),'MA{ }'.format(q)]=np.nan        continuetry:model =sm.tsa.ARIMA(train,order=(p,d,q))results =model.fit()results_bic.loc['AR{ }'.format(p),'MA{ }'.format(q)]=results.bic    except:continueresults_bic

运行结果：

绘制热力图：

#得到结果后进行浮点型转换results_bic =results_bic[results_bic.columns].astype(float)#绘制热力图fig,ax =plt.subplots(figsize=(10,8))ax =sns.heatmap(results_bic,mask=results_bic.isnull(),ax=ax,annot=True,fmt='.2f',cmap="Purples")ax.set_title('BIC')plt.show()results_bic.stack().idxmin()

运行结果：

结果分析和模型检验

残差序列的随机性可以通过自相关函数法来检验,即做残差的自相关函数图

# 模型检验#根据以上求得p =1d =0q =0model =sm.tsa.ARIMA(train,order=(p,d,q))results =model.fit()resid =results.resid #获取残差#绘制#查看测试集的时间序列与数据(只包含测试集)fig,ax =plt.subplots(figsize=(12,5))ax =sm.graphics.tsa.plot_acf(resid,lags=40,ax=ax)plt.show()

运行结果：

结果分析：从 ACF 图中可以看出残差之间独立性比较高。

predict_sunspots =results.predict(dynamic=False)print(predict_sunspots)#查看测试集的时间序列与数据(只包含测试集)plt.figure(figsize=(12,6))plt.plot(train)plt.xticks(rotation=45)#旋转45度plt.plot(predict_sunspots)plt.show()

运行结果：

结果分析：

该模型拟合效果较好。

模型预测

# 进行预测，预测范围未来三个月start_index =len(train)end_index =len(train)+len(test)-1forecast =results.predict(start=start_index,end=end_index,dynamic=False)# 将预测结果与原始数据序列合并forecast.index =test.index# 绘制结果plt.figure(figsize=(12,6))plt.plot(train.index,train,label='Train')plt.plot(test.index,test,label='Test')plt.plot(forecast.index,forecast,label='Forecast',linestyle='--')# plt.plot(predict_sunspots, label='Fit')plt.legend()plt.xticks(rotation=45)plt.show()

运行结果：

向后预测了三个月的数据。

代码改进见博客：

时间序列模型(ARIMA) — — 预测未来(含 python 和 Matlab 的完整代码)