数学建模：控制预测类——时间序列ARIMA模型

见贤思齐

目录1.时间序列ARIMA模型2.ARIMA模型大纲3.模型详解1）自回归模型AR(p)2）移动平均模型MA(q)3）自回归移动平均模型ARMA(p,q)4）差分自回归移动平均模型ARIMA(p,d,q)4.ARIMA模型建模步骤5.建模步骤名词解释1）平稳性2）ADF检验3）差分法4）拖尾和截尾5）自相关系数（ACF）6）偏自相关系数（PACF）7）AIC、BIC准则6.ARIMA模型案例演示1）导入数据2）提取要预测的列3）划分训练测试集4）差分法5）ADF检验6）对训练集通过ACF、PACF确定参数7）AIC、BIC准则寻找最优参数p、q8）模型检验（绘制残差ACF图）9）模型预测1.时间序列ARIMA模型时间序列也称动态序列，是指将某种现象的指标数值按照时间顺序排列而成的数值序列。时间序列可分成三大部分：描述过去、分析规律和预测未来。 2.ARIMA模型大纲 3.模型详解1）自回归模型AR(p)描述当前值和历史值之间的关系，用变量自身的历史数据对自身进行预测，其必须要满足平稳性要求，只适用于预测与自身前期相关的现象（时间序列的自相关性）。p阶自回归过程公式定义：，p表示用几期的历史值来预测。（其中：当前值；：常数项；p：阶数；：自相关系数；：误差项）2）移动平均模型MA(q)移动平均模型关注的是自回归模型中误差项的累计。q阶自回归过程公式定义：，即时间序列当前值与历史值没有关系，而只依赖于历史白噪声（随机误差）的线性组合。移动平均法能有效消除预测中的随机波动。3）自回归移动平均模型ARMA(p,q)公式定义：，该式表明：一个随机时间序列可以通过一个自回归移动平均模型来表示，即该序列可以由其自身的过去或滞后值以及随机扰动项来解释。如果该序列是平稳的，即它的行为（变化规律、波动）不会随着时间的推移而变化，那就可通过该序列过去的行为来预测未来。4）差分自回归移动平均模型ARIMA(p,d,q)p：自回归项；q：移动平均项数；d：时间序列成为平稳时所做的差分次数。原理：将非平稳时间序列转换为平稳时间序列，然后将因变量仅对它的滞后值（自回归模型）以及随机误差项的现值和滞后值（移动平均模型）进行回归所建立的模型。 4.ARIMA模型建模步骤1）对序列绘图，进行平稳性检验（ADF检验），观察序列是否平稳；对于非平稳时间序列要先进行d阶差分，转化为平稳时间序列。2）对平稳时间序列分别求其自相关系数（ACF）和偏自相关系数（PACF），通过对自相关图和偏自相关图的分析或通过AIC/BIC搜索，得到最佳阶数p、q。3）通过以上的d、p、q得到ARIMA模型。然后对得到的模型进行模型校验。 5.建模步骤名词解释1）平稳性平稳性要求经由样本时间序列所得到的拟合曲线在未来的一段时间内仍能够按照现有形态延续下去。平稳性要求序列的均值和方差不发生明显变化。严平稳：系列所有统计性质（期望、方差）都不随时间推移而变化。宽平稳：期望与相关系数（依赖性）不变，就是说t时刻的值X依赖于过去的信息。2）ADF检验检验时间序列是否满足平稳性要求。ADF大致思想就是基于随机游走（不平稳的一个特殊序列）的，对其进行回归，如果发现p=1，说明序列满足随机游走，就是非平稳的。ADF检验结果共有5个参数（第6个参数不做要求）：第1个值：TestStatistic，表示T统计量，即假设检验值；第2个值：p-value，即p值，表示T统计量对应的概率值；第3/4个值：表示延迟和测试的次数；第5个参数{‘1%’：xxx,‘5%’：xxx,‘10%’：xxx}：表示不同程度拒绝原假设的统计值。那么该如何确定该序列是否平稳呢？将1%、5%、10%不同程度拒绝原假设的统计值和ADF假设检验值进行比较，如果ADF假设检验值同时小于 1%、5%、10%，即说明非常好地拒绝原假设，该序列是平稳的。看p-value是否非常接近0，如果非常接近0，则拒绝原假设，确定该序列是平稳的。3）差分法时间序列在t和t-1时刻的差值。，（x=0,1,2...）例：一组数列[0,1,2,3,4,5,6,7]进行差分后就会得到新数列[1,1,1,1,1,1,1]。由此可以得到d（时间序列成为平稳时所做的差分次数）。4）拖尾和截尾拖尾指序列以指数率单调递减或震荡衰减，而截尾指序列从某个时点变得非常小。具体判断如下：拖尾：①有超过5%的样本(偏)自相关系数都落入2倍标准差范围之外；②或者由显著非0的(偏)自相关系数衰减为小值波动的过程比较缓慢或非常连续。截尾：①在最初的d阶明显大于2倍标准差范围；②之后几乎95%的(偏)自相关系数都落在2倍标准差范围以内；③且由非零自相关系数衰减为在零附近小值波动的过程非常突然。5）自相关系数（ACF）有序的随机变量序列与其自身相比较。自相关系数度量的是同一事件在两个不同时期之间的相关程度。公式：，取值范围[-1,1]。（即：对于时间序列，与的相关系数称为间隔k的自相关系数）。差分数据自相关图（ACF）从上图可以看出，趋势序列ACF有1阶截尾，因此可以选p=1。6）偏自相关系数（PACF）把其他要素的影响视为常数，即暂不考虑其他要素影响，单独研究两个要素之间的相互关系的密切程度，称为偏相关。 公式：差分数据偏自相关图（PACF）从上图可以看出，PACF有1阶截尾，因此可以选q=1。7）AIC、BIC准则AIC准则即最小化信息量准则，BIC准则即贝叶斯信息准则。这两个指标越小越好。然而由于通过自相关分析和偏自相关分析来判断ARIMA的参数存在人为主观性，为避免这一影响，那么就可基于AIC、BIC信息准则自动寻找最优参数。 6.ARIMA模型案例演示1）导入数据importpandasaspdimportmatplotlib.pyplotaspltChinaBank=pd.read_csv("D:\ChinaBank.csv",index_col='Date',parse_dates=['Date'])#parse_dates=['Date']:将Date列解析为日期时间类型（datetime）2）提取要预测的列这里提取Close列：ChinaBank.index=pd.to_datetime(ChinaBank.index)sub=ChinaBank.loc['2014-01':'2014-06','Close']#选取索引值从'2014-01'到'2014-06'的所有行，以及列名为'Close'的列sub.head()3）划分训练测试集train=sub.loc['2014-01':'2014-03']test=sub.loc['2014-04':'2014-06']#查看训练集的时间序列与数据（只包含训练集）plt.figure(figsize=(12,6))plt.plot(train)plt.xticks(rotation=45)#旋转45度plt.show()训练集的时间序列与数据如下：4）差分法#.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)#add_subplot方法用于在这个图形窗口中添加一个子图网格#311用于指定子图的布局,3表示行，1表示列，1表示位置编号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数值进行判断），进行差分之后，均值和方差就不发生明显变化了，满足平稳性要求。5）ADF检验importmatplotlib.pyplotaspltimportpandasaspdimportstatsmodels.apiassm#statsmodels提供了许多用于估计和检验统计模型的类和函数，以及进行统计测试的数据和结果的可视化。fromstatsmodels.tsa.seasonalimportseasonal_decompose#seasonal_decompose函数用于时间序列的季节性分解，将时间序列数据分解为趋势（trend）、季节性（seasonal）和残差（residual）成分。fromstatsmodels.tsa.stattoolsimportadfullerasADF#计算原始序列、一阶差分序列、二阶差分序列的单位根检验结果ChinaBank['diff_1']=ChinaBank['diff_1'].fillna(0)#.fillna(0)：填充缺失值ChinaBank['diff_2']=ChinaBank['diff_2'].fillna(0)timeseries_adf=ADF(ChinaBank['Close'].tolist())#.tolist()：将ChinaBank['Close']转换为一个Python列表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)运行以上代码，结果如下，可以发现进行1阶差分后数据很好地满足平稳性要求：timeseries_adf：(0.527919808483182,0.9856974415734416,9,335,{'1%':-3.4500219858626227,'5%':-2.870206553997666,'10%':-2.571387268879483},-734.0738716811488)timeseries_diff1_adf：(-6.177185544979001,6.58710923976123e-08,8,336,{'1%':-3.449962981927952,'5%':-2.870180642420163,'10%':-2.5713734527352607},-735.8436797171294)timeseries_diff2_adf：(-9.202545123160352,1.9841232339615405e-15,13,331,{'1%':-3.4502615951739393,'5%':-2.8703117734117742,'10%':-2.5714433728242714},-717.2833732193085)6）对训练集通过ACF、PACF确定参数importstatsmodels.apiassm#绘制fig=plt.figure(figsize=(12,7))ax1=fig.add_subplot(211)fig=sm.graphics.tsa.plot_acf(train,lags=20,ax=ax1)#lags=20:这个参数指定了要计算并绘制自相关性的滞后阶数。它会计算并显示从0阶（即当前值与自身的相关性，总是1）到20阶的自相关性#ax=ax1:这个参数允许你将绘制的图形添加到指定的Matplotlib轴（axis）对象上。ax1.xaxis.set_ticks_position('bottom')#设置坐标轴上的数字显示的位置#fig.tight_layout()ax2=fig.add_subplot(212)fig=sm.graphics.tsa.plot_pacf(train,lags=20,ax=ax2)ax2.xaxis.set_ticks_position('bottom')#设置坐标轴上的数字显示的位置#fig.tight_layout()plt.show()可以观察到序列是比较平稳的，基本都落在置信区间内，这里训练集d就取0，其中ACF图为2阶拖尾，PACF图为1阶截尾，因此可以确定p=2,q=1，为避免人为主观性影响，这里我们还是通过AIC、BIC准则来寻找最优参数p、q。7）AIC、BIC准则寻找最优参数p、q#寻找最优参数模型建立#遍历，寻找适宜的参数importitertools#itertools包含了一组创建迭代器的函数importnumpyasnpimportseabornassns#seaborn是一个基于matplotlib的Python数据可视化库#确定pq的取值范围p_min=0d_min=0q_min=0p_max=5d_max=0q_max=5#初始化dataframe用来存储结果，以BIC准则results_bic=pd.DataFrame(index=['AR{}'.format(i)foriinrange(p_min,p_max+1)],columns=['MA{}'.format(i)foriinrange(q_min,q_max+1)])forp,d,qinitertools.product(range(p_min,p_max+1),range(d_min,d_max+1),range(q_min,q_max+1)):ifp==0andd==0andq==0:#将results_bicDataFrame中对应的单元格设置为np.nanresults_bic.loc['AR{}'.format(p),'MA{}'.format(q)]=np.nancontinuetry:#创建一个ARIMA模型实例model=sm.tsa.ARIMA(train,order=(p,d,q),#enforce_stationarity=False,#enforce_invertibility=False)results=model.fit()#调用模型的fit方法来拟合模型到训练数据上#将拟合结果的BIC值存储在results_bicDataFrame的相应单元格中。单元格的位置由'AR'后跟p的值和'MA'后跟q的值来确定。results_bic.loc['AR{}'.format(p),'MA{}'.format(q)]=results.bicexcept:continue#得到结果后进行浮点型转换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(),#mask参数设置为results_bic中的空值（NaN），这将在热力图中用白色（或透明，取决于cmap）表示这些位置ax=ax,annot=True,#在每个单元格上显示数值fmt='.2f',#表示数值将格式化为保留两位小数的浮点数cmap="Purples")ax.set_title('BIC')plt.show()#返回最小值索引results_bic.stack().idxmin()可以得到结果p=1,q=0： 也可通过statsmodels库获取p、q最优值：#通过AIC和BIC来求p,q最优值#statsmodels库中的tsa.arma_order_select_ic函数基于AIC和BIC选择自回归移动平均（ARMA）模型的最佳p（自回归项的阶数）和q（移动平均项的阶数）值train_results=sm.tsa.arma_order_select_ic(train,ic=['aic','bic'],trend='n',max_ar=5,max_ma=5)#train:要拟合ARMA模型的时间序列数据#ic=['aic','bic']:指定使用的信息准则列表#trend='n':指定模型中是否包含趋势项。'n'表示不包含趋势项（即只考虑ARMA模型，不考虑趋势）。#max_ar=8,max_ma=8:最大p为5和最大q为5print('AIC',train_results.aic_min_order)print('BIC',train_results.bic_min_order)8）模型检验（绘制残差ACF图）#根据以上求得p=1d=0q=0#创建一个ARIMA模型实例model=sm.tsa.ARIMA(sub,order=(p,d,q))results=model.fit()resid=results.resid#获取残差#绘制残差ACF图#查看所选取的整个数据的时间序列与数据fig,ax=plt.subplots(figsize=(12,5))ax=sm.graphics.tsa.plot_acf(resid,lags=40,ax=ax)plt.show()可以发现大多数自相关系数都落在置信区间内，说明残差是独立的，时间序列ARIMA模型适用。9）模型预测predict_sunspots=results.predict(dynamic=False)#results.predict(...):模型拟合结果对象results的一个方法，用于生成未来时间点的预测值。#dynamic=False:预测是基于历史数据（即直到最后一个观测值为止的所有数据）进行的，而不是基于之前步骤的预测结果。print(predict_sunspots)#查看时间序列与数据（蓝色代表原始数据，黄色代表预测数据）plt.figure(figsize=(12,6))plt.plot(sub)plt.xticks(rotation=45)#旋转45度plt.plot(predict_sunspots)plt.show()