8.7 在R中建立ARIMA模型

auto.arima()函数是如何实现的?

R中的auto.arima()函数通过使用结合多种单位根检验的Hyndman-Khandakar算法(Hyndman & Khandakar, 2008),最小化AICc值及最大似然估计来建立合适的ARIMA模型。 它的算法由以下几个步骤组成:

ARIMA建模的Hyndman-Khandakar算法

  1. \(0 \le d\le 2\)时差分的取值可以使用重复KPSS检验来确定。
  2. \(p\)\(q\)的取值可以通过差分数据\(d\)次后最小化AICc值来得到.与测试所有\(p\)\(q\)的所有可能取值不同,本算法使用逐步搜索来遍历模型空间。
    1. 四个初试拟合模型:

      • ARIMA\((0,d,0)\),
      • ARIMA\((2,d,2)\),
      • ARIMA\((1,d,0)\),
      • ARIMA\((0,d,1)\).

      除非\(d=2\),不然模型将会包含一个常数.如果\(d \le 1\),则还需要添加一个附加模型来进行拟合:

      • 不含常数的ARIMA\((0,d,0)\)模型.
    2. 将在(a)中得到的AICc值最小的最优模型设定为“当前模型”.

    3. 通过当前模型生成衍生模型并用以比较:

      • 对当前模型的\(p\)和/或\(q\)进行\(\pm1\)的操作;
      • 清除或者添加当前模型的常量\(c\). 目前的最优模型将会在衍生模型和当前模型中选出(选择AICc最小的一个),将选出的模型作为 新的当前模型.
    4. 重复2(c)环节直到没有更小的AICc值出现.

auto.arima()的函数参数可以对算法进行修正和改进。上述步骤是在默认参数下的过程。

默认的步骤中采用了一些近似方法来加速搜索,可以设置参数approximation=FALSE不使用这些近似方法,可能会出现由于这些近似方法或逐步变量筛选(stepwise)而无法找到最小AICc的情况。如果设置参数 stepwise = FALSE,会有更多的模型被搜索。所有参数的详细描述请查看帮助文档。

选择自己的模型

如果你想自己选择模型,你可以使用R中的Arima()函数。R中的arima()函数也可以用来拟合ARIMA模型,但是它不适用于常数 \(c\) 存在的情况(除非 \(d=0\)),并且它和forecast包中的一些函数不兼容。最后,它的估计模型也不能被用于新数据集(新数据在检测预测准确率时非常有用)。因此我们强烈建议使用Arima()而不是arima()

建模过程

使用 ARIMA 模型对非季节性时间序列进行拟合时,下列过程是较为通用的方法:

  1. 画出数据时序图并检查有无异常观测。
  2. 如果必要的话,对数据进行变换(如 Box-Cox 变换)来稳定方差。
  3. 如果数据非平稳,对数据进行差分直到数据平稳。
  4. 检查自相关图和偏自相关图:数据是否符合ARIMA(\(p,d,0\))或者ARIMA(\(0,d,q\))模型的特征?
  5. 尝试拟合你选择的模型,并使用AICc来进行模型选择。
  6. 通过画出残差自相关图来检查你选择模型的残差,并且对残差进行一元混成检验(portmanteau test)。如果它们看起来不像白噪声,那么模型就需要修正。
  7. 如果残差看起来类似白噪声,进行模型预测。

Hyndman-Khandakar算法 只能完成第 3 到第 5 步骤。因此即使你使用了这个算法,你仍然需要自行完成其他步骤。

8.11概括了所需的过程。

使用ARIMA模型进行预测的通用过程.

图 8.11: 使用ARIMA模型进行预测的通用过程.

实例:季节调整后的电器订单

我们将要使用上文中所讨论的过程来处理季节调整后的电器订单数据,如图8.12所示:

elecequip %>% stl(s.window='periodic') %>% seasadj() -> eeadj
autoplot(eeadj)
欧元区季节调整后的电器订单指数.

图 8.12: 欧元区季节调整后的电器订单指数.

  1. 时序图中显示出一些突然的变化,尤其是2008/2009年的急剧下降。这些变化是全球经济形势的变化导致的,除此之外时序图中没有其他异常,因此不用对数据进行调整。
  2. 图中并没有方差变化的迹象,因此我们不需要对数据进行Box-Cox变换。
  3. 图中序列在长期上不断上升下降,显示出非常明显的非平稳性。因此,我们对数据进行一阶差分,图8.13显示的是差分后的数据。它们看起来较为平稳,因此我们不再进行进一步的差分。
eeadj %>% diff() %>% ggtsdisplay(main="")
差分后的欧元区季节调整后的电器订单指数的时序图,自相关图以及偏自相关图。

图 8.13: 差分后的欧元区季节调整后的电器订单指数的时序图,自相关图以及偏自相关图。

  1. 8.13中显示的偏自相关图类似于AR(3)模型,因此我们将最初模型设定为ARIMA(3,1,0),除此之外没有明显的备选模型。

  2. 我们使用 ARIMA(3,1,0)以及衍生模型ARIMA(4,1,0)、ARIMA(2,1,0)、ARIMA(3,1,1)等等进行拟合,在这些模型中,ARIMA(3,1,1)是AICc值最小的模型。

    fit <- Arima(eeadj,order=c(3,1,1))
    summary(fit)
    #> Series: eeadj 
    #> ARIMA(3,1,1) 
    #> 
    #> Coefficients:
    #>         ar1    ar2    ar3     ma1
    #>       0.004  0.092  0.370  -0.392
    #> s.e.  0.220  0.098  0.067   0.243
    #> 
    #> sigma^2 estimated as 9.58:  log likelihood=-492.7
    #> AIC=995.4   AICc=995.7   BIC=1012
    #> 
    #> Training set error measures:
    #>                   ME  RMSE   MAE      MPE  MAPE   MASE
    #> Training set 0.03288 3.055 2.357 -0.00647 2.482 0.2884
    #>                  ACF1
    #> Training set 0.008981
  3. ARIMA(3,1,1)模型的残差自相关图显示出所有的自相关系数都在置信域之内,这反映出残差类似于白噪声。一元混成检验得到了一个较大的p值,同样意味着残差类似白噪声。

    checkresiduals(fit)

    #> 
    #>  Ljung-Box test
    #> 
    #> data:  Residuals from ARIMA(3,1,1)
    #> Q* = 24, df = 20, p-value = 0.2
    #> 
    #> Model df: 4.   Total lags used: 24
  4. 8.14显示的是选择的模型的预测值。

    autoplot(forecast(fit))
    季节调整后的电器订单指数预测.

    图 8.14: 季节调整后的电器订单指数预测.

如果我们使用自动的算法的话,我们将会在默认设置下得到一个ARIMA(3,1,0)模型,如果我们approximation=FALSE,则会得到一个ARIMA(3,1,1)模型。

理解R中的常量

一个非季节性的ARIMA模型可以被表示为: \[\begin{equation} \tag{8.4} (1-\phi_1B - \cdots - \phi_p B^p)(1-B)^d y_t = c + (1 + \theta_1 B + \cdots + \theta_q B^q)\varepsilon_t, \end{equation}\] 或者也可以表示为: \[\begin{equation} \tag{8.5} (1-\phi_1B - \cdots - \phi_p B^p)(1-B)^d (y_t - \mu t^d/d!) = (1 + \theta_1 B + \cdots + \theta_q B^q)\varepsilon_t, \end{equation}\] 这里的 \(c = \mu(1-\phi_1 - \cdots - \phi_p )\)\(\mu\)\((1-B)^d y_t\) 的均值。R 使用的是等式(8.5)中的参数表示。

因此,在非平稳的 ARIMA 模型中加入常量等同于在预测函数中包含一个 \(d\) 阶的多项式趋势项(polynomial trend)。(如果我们忽略掉常数,预测函数将会包含一个 \(d-1\) 阶的多项式趋势项。)当 \(d=0\) 时,我们可以得到 \(\mu\) 等于 \(y_t\) 的均值的特例。

在默认设置下,Arima()函数会在 \(d>0\) 设置 \(c=\mu=0\) 并且在 \(d=0\) 时给出\(\mu\)的估计。它和样本均值很接近,但是通常并不等于样本均值,因为样本均值在 \(p+q>0\) 并不等于极大似然估计值。

include.mean参数只在 \(d=0\) 时发挥作用,并且默认为真,将它设置为假将会使得\(\mu=c=0\)

include.drift参数可以允许 \(d=1\)\(\mu\ne0\) 。当 \(d>1\) 时,常量是不能允许存在的,因为它引起的二阶或高阶趋势会使得预测受到过大影响。当 \(d=1\) 时,参数会在R的输出中被称为“drift”。

这里还有 一个参数include.constant,如果为真将会在 \(d=0\) 时设置include.mean=TRUE,在 \(d=1\) 时设置include.drift=TRUE。如果include.constant=FALSE,则include.meaninclude.drift都会被设置为 FALSE。一旦使用了include.constant,则include.mean=TRUEinclude.drift=TRUE都会被忽略。

auto.arima()函数 对常量的取舍进行了自动处理。在默认设置下,对于 \(d=0\)\(d=1\) 的情况,如果常量的引入可以降低AICc值,模型将包含这个常量。在\(d>1\)的情况下常量将被忽略。如果专门设置了allowdrift=FALSE,则只有在 \(d=0\) 时才会考虑常量。

绘制特征根

(这是更深入的章节,你可以自行决定是否跳过这一部分。)

我们可以改写等式(8.4)\[\phi(B) (1-B)^d y_t = c + \theta(B) \varepsilon_t\] 这里的 \(\phi(B)= (1-\phi_1B - \cdots - \phi_p B^p)\)\(B\) 中的一个 \(p\) 阶多项式,\(\theta(B) = (1 + \theta_1 B + \cdots + \theta_q B^q)\)\(B\) 的一个 \(q\) 阶多项式。

模型的平稳性条件为 \(\phi(B)\)\(p\) 复数根在单位圆之外,并且可逆性的条件为 \(\phi(B)\)\(q\) 复数根在单位圆之外,因此我们可以通过图中根与复数单位圆的关系判断模型的平稳性和可逆性。

画出逆根(inverse roots)是一种更简便的办法,因为它们肯定在单位圆之内,这在R里可以很容易完成。对于拟合季节调整后的电器指数的ARIMA(3,1,1)模型,我们可以得到图8.15

autoplot(fit)
拟合季节调整后的电器指数的ARIMA(3,1,1)模型的逆特征根

图 8.15: 拟合季节调整后的电器指数的ARIMA(3,1,1)模型的逆特征根

左图中的三个红点对应的是多项式 \(\phi(B)\) 的根,而右图红点对应的是多项式 \(\theta(B)\) 的根。正如我们所期待的那样,它们都在单位圆之内,因为 R 保证了拟合的模型一定是平稳且可逆的。任何接近单位圆边界的根都有可能是不稳定的,它所对应的模型将不适合做预测。

Arima()函数不会逆根在单位圆以外的模型。auto.arima()函数甚至更严格,连一个根接近单位圆的模型都不会选择。

参考文献

Hyndman, R. J., & Khandakar, Y. (2008). Automatic time series forecasting: The forecast package for R. Journal of Statistical Software, 27(1), 1–22. https://doi.org/10.18637/jss.v027.i03