原文:https://otexts.com/fppcn/arima-r.html
auto.arima()
函数是如何实现的?The auto.arima()
function in R uses a variation of the Hyndman-Khandakar algorithm (Hyndman and Khandakar 2008), which combines unit root tests, minimisation of the AICc and MLE to obtain an ARIMA model. The arguments to auto.arima()
provide for many variations on the algorithm. What is described here is the default behaviour.
Hyndman-Khandakar algorithm for automatic ARIMA modelling |
---|
|
|
|
|
|
|
默认的步骤中采用了一些估计方法来加速搜索,可以设置参数approximation=FALSE
不使用这些估计方法,可能会出现由于因为这些估计方法或步进(stepwise)而无法找到最小的 AICc的情况。如果设置参数 stepwise = FALSE
,会有更多的模型被搜索。所有参数的详细描述请查看帮助文档。
R中的auto.arima()
函数通过使用结合多种单位根检验的Hyndman-Khandakar算法(Hyndman and Khandakar 2008),最小化AICc值以及最小二乘法来建立合适的ARIMA模型. 它的算法由以下几个步骤组成:
ARIMA建模的Hyndman-Khandakar算法
四个初试拟合模型:
除非\(d=2\),不然模型将会包含一个常数.如果\(d \le 1\),则还需要添加一个附加模型来进行拟合:
重复2(c)环节直到没有更小的AICc值出现.
auto.arima()
的函数参数可以对算法进行修正和改进.上述过程是在默认参数下的过程.
默认的过程使用了一些近似来提高搜索效率.approximation=FALSE
语句可以去除近似的部分.是近似或者逐步搜索都有可能忽略掉拥有最小AICc的模型,因此stepwise=FALSE
语句可以扩大搜索的范围.帮助文档中有更多参数设置的描述.
如果你想自己选择模型,你可以使用R中的Arima()
函数。R中的arima()
函数也可以用来拟合ARIMA模型,但是它不适用于常数 \(c\) 存在的情况(除非 \(d=0\)),并且它和forecast
包中的一些函数不兼容。最后,它的估计模型也不能被用于新数据集(新数据在检测预测准确率时非常有用)。因此我们强烈建议使用Arima()
而不是arima()
。
使用 ARIMA 模型对对非季节性时间序列进行拟合时,下列过程是较为通用的方法:
Hyndman-Khandakar算法 只能完成第 3 到第 5 步骤。因此即使你使用了这个算法,你仍然需要自行完成其他步骤。
图8.11概括了所需的过程。
图 8.11: 使用ARIMA模型进行预测的通用过程.
我们将要使用上文中所讨论的过程来处理季节调整后的电器订单数据,如图8.12所示:
elecequip %>% stl(s.window='periodic') %>% seasadj() -> eeadj
autoplot(eeadj)
图 8.12: 欧元区季节调整后的电器订单指数.
eeadj %>% diff() %>% ggtsdisplay(main="")
图 8.13: 差分后的欧元区季节调整后的电器订单指数的时序图,自相关图以及偏自相关图。
我们使用 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 ACF1
#> Training set 0.03288 3.055 2.357 -0.00647 2.482 0.2884 0.008981
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
图8.14显示的是选择的模型的预测值。
autoplot(forecast(fit))
图 8.14: 季节调整后的电器订单指数预测.
如果我们使用自动的算法的话,我们将会在默认设置下得到一个ARIMA(3,1,0)模型,如果我们approximation=FALSE
,则会得到一个ARIMA(3,1,1)模型。
一个非季节性的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.mean
和include.drift
都会被设置为 FALSE。一旦使用了include.constant
,则include.mean=TRUE
和include.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)
图 8.15: 拟合季节调整后的电器指数的ARIMA(3,1,1)模型的逆特征根
左图中的三个红点对应的是多项式 \(\phi(B)\) 的根,而右图红点对应的是多项式 \(\theta(B)\) 的根。正如我们所期待的那样,它们都在单位圆之内,因为 R 保证了拟合的模型一定是平稳且可逆的。任何接近单位圆边界的根都有可能是不稳定的,它所对应的模型将不适合做预测。
Arima()
函数不会逆根在单位圆以外的模型。auto.arima()
函数甚至更严格,连一个根接近单位圆的模型都不会选择。