时间序列的R语言实现_Part3

ARIMA模型

前面讲到的指数平滑模型，对于预测误差要求是无相关性并且结果是满足以0为中心的正态分布的。对时间序列作分析是，有些情况下，把相关性考虑进去能做更好的预测。ARIMA模型就包含了对时间序列的无规律部分做预测的模型，在无规律部分，允许出现非零相关性。

 

对时间序列做差分

ARIMA模型是针对稳定时间序列的。若是你的分析对象是非稳定的时间序列，那么你首先就需要对时间序列做差分，直到你得到一个稳定的时间序列。需要做差分时，我们用到的模型是ARIMA(p,d,q)，其中d的值是一个大于等于0的整数值，它就是我们对原始时间序列所做差分的次数。

R中有做差分的方法diff()，例如对前面用到的数据，1866年到1911年女性裙子的边缘直径数据，不是一个稳定的时间序列，数据的水平值随时间变化很大。

 

对这个数据做差分，过程结果如下：

 

diff()差分方法中，用到的两个参数，第一个是做差分的对象，即原始时间序列，后面diferences是所做差分次数

从上面的结果图能看出，该时间序列一次差分的结果依然不是稳定的时间序列。所以对该时间序列做二次差分。

 

这一次得到的时间序列，随时间的变化存在不变方差，可以认为是稳定的。所以对1866年到1911年女性裙子的边缘直径数据这个时间序列，需要对其做两次差分以得到稳定的时间序列。它所适用的ARIMA模型就是ARIMA(p,2,q)，前面讲到过，2就是需要对时间序列所做的差分的次数，接下来我们需要做的是确定p和q的值。

 

另一个时间序列的例子是，前面用到过的英国连续42任国王的逝世年龄。下面是对该时间序列做差分的过程：

 

一次差分后的结果可以看出就基本稳定了，所以这个例子适用的ARIMA模型就是ARIMA(p,1,q)。

 

选择候选ARIMA模型

前面讲到了ARIMA(p,d,q)模型d的值就是要得到稳定时间序列所要对原始时间序列做差分的次数。选择合适的p和q值则需要对稳定的时间序列做相关和偏自相关的测试。

R中有acf()和pacf()方法能得到数据的自相关图和偏自相关图。如果要查看自相关和偏自相关结果的值，设置参数plot=FALSE即可。下面是实际实现过程：

 

英国连续42任国王的例子

自相关过程及结果：

自相关的结果图上，只有lag1的结果超过了意义界限，其余都在接线之内。

偏自相关过程及结果如下：

 

偏自相关的结果中，lag为1,2和3的时候，结果超过了意义界限，这三个点都是负值，并且逐渐递减(lag 1: -0.360, lag 2: -0.335, lag 3:-0.321)。偏自相关的结果是3阶截尾的。

这个例子时间序列的自相关图1阶截尾，偏自相关图3阶截尾，就是说，下面的ARMA模型对一次差分后的时间序列都是合适的：

1.ARMA(3,0)模型，相当于AR(3)模型，因为偏自相关图3阶截尾，自相关图可视为拖尾。

2.ARMA(0,1)模型，相当于MA(1)模型，因为自相关图1阶截尾，偏自相关图可视为拖尾。

4.ARMA(p,q)模型，因为自相关图和偏自相关图都可视为拖尾。

模型选择的原则是精简原则，是参数尽量少，所以ARMA(0,1)相比而言是最合适的一个模型。

一次差分后的英国国王逝世年龄时间序列适用ARMA(0,1)模型，也就是说原始时间序列，适用ARIMA(0,1,1)模型。

 

北半球火山尘的例子

这是一个新的例子，来看如何选择合适的ARIMA模型。数据是北半球从1500年到1969年的火山尘指数值。读取存储和查看这个时间序列数据的过程如下：

 

从上图可以看出，时间序列的随机波动幅度随时间保持在一个稳定的值，所以加性模型可以用来描述这个模型。且该时间序列的差分值大致稳定，所以对这个时间序列不需要做差分处理。也就是ARIMA模型的d值为0。

接下来看时间序列的滞后期1-20的自相关图和偏自相关图来确定ARIMA模型的选择：

 

自相关图能看出来，在滞后期为1,2和3时自相关值超出了意义界限，在3之后，自相关值接近0，自相关值在滞后期为1,2,3时值为正值，且逐渐递减。滞后期为19和20时，自相关值同样超出了意义界限，但超出很少，可以认为是偶然出现的结果。

 

偏自相关的结果是2阶截尾，同样，滞后期为19时，偏自相关值超出意义界限可视为偶然出现结果不考虑其影响。

综上，自相关图3阶截尾，偏自相关图2阶截尾，所以该时间序列适用的模型有下面几种：

1.ARMA(2,0)模型，因为偏自相关图2阶截尾，自相关图可视为拖尾。

2.ARMA(0,3)模型，因为自相关图3阶截尾，偏自相关图可视为拖尾。

3.ARMA(p,q)模型，自相关图和偏自相关图都可视为拖尾。

同样，根据精简原则，参数尽量少的话，ARMA(p,q)参数至少为p+q>=2，所以，ARMA(1,0)模型和ARMA(p,q)模型都可考虑。

对本半球火山尘时间序列数据，适合的ARIMA模型是ARIMA(2,0,0)或者是ARIMA(p,0,q)。

说明：ARIMA算法是根据自相关图和偏自相关图的拖尾和截尾的性质来判断如何选择具体的ARIMA模型。但是拖尾和截尾的判断方法我没有找到一个很明确的说法，照这里看来判断过程带有一定的主观性，这个例子中，多个模型都可以用。另外原文中tail off to zero我理解为拖尾，is zero after lag n我理解为n阶截尾，这个不知道确不确切，欢迎大家指正，谢谢！

 

用ARIMA模型做预测

前面用到了两个例子，并且已经确定了这两个例子比较适用的ARIMA模型，R中提供了arima()方法。

 

英国连续42任国王的例子

对这个时间序列，前面已经确认过ARIMA(0,1,1)是比较适合的模型。arima()方法中的order参数可定义p,d,q的值，arima()方法创建ARIMA对象过程如下：

 

上面显示的kingstimeseriesarima对象，下面的值是ARIMA模型中用到的一些参数的预估值。

forecast包中的forecastArima()方法可用来对arima对象做预测。

上面h=5，表示做了5期的预测，原始数据是连续42任英国国王的逝世年龄，所以预测结果是43-47这5期的值。包括80%置信区间的上限和下限，95%置信区间的上限和下限。

同样，可以在图中查看预测结果：

 

跟在用指数平滑法做预测一样，我们可以计算预测误差看其是否以0为中心呈正态分布，是否存在不变方差，预测误差是否存在相关性。

所以，我们先看预测误差的相关性图和Ljung-Box测试结果：

 

相关性图中可看到，1-20的滞后期中，没有自相关值超出意义界限，Ljung-Box测试p值为0.9，我们可以据此推断预测误差在1-20的滞后期中，不存在非零自相关性。

检验预测误差是否以0为中心呈正态分布，是否存在不变方差，可以在时序图和直方图中查看预测误差。

 

 

从预测误差的时序图可以看出，其方差随时间大致相同；直方图能看出预测误差也是存在不变方差，满足正态分布的。

所以，ARIMA(0,1,1)对连续42任英国国王的逝世年龄的预测是有效的。

 

北半球火山尘的例子

前面对这个例子可适用的ARIMA模型做了判断，ARIMA(2,0,0)是比较合适的一个。得到ARIMA对象的过程如下：

再对其做预测，过程及结果如下：

在图中查看预测结果：

 

结果中，有一个问题，预测存在负值，但是实际的火山尘指数不可能为负值。arima()和forecast.Arima()方法无法定义结果只为正值，这是这个方法不太让人满意的地方。

同样，我们需要检验预测误差，看其是否存在非0相关性，并且是否以0为中心呈正态分布。同样看预测误差的自相关图，并给预测误差做Ljung-Box测试：

 

结果可以确定预测误差存在非0自相关性可能很小。

接下来看预测误差是否存在不变方差满足正态分布。同样查看预测误差的时序图，并查看其与正态分布做对比的直方图。过程结果如下：

 

 

 

同样这两个图可以大致确认，预测误差是存在不变方差，并且满足以0为中心的正态分布。

所以ARIMA(2,0,0)模型对原时间序列所做的预测是有效的。

说明：前面判断知道这个时间序列可适用ARIMA(2,0,0)模型，也可选择ARIMA(p,0,q)模型，按精简参数原则可以尝试使用ARIMA(1,0,1)模型重复这个过程，得到的预测结果与ARIMA(2,0,0)非常接近。