8.5 非季节性ARIMA模型
当我们将差分和自回归模型以及移动平均模型结合起来的时候，我们可以得到一个非季节性 ARIMA 模型。ARIMA 是 AutoRegressive Integrated Moving Average 的简称。（在这里Integrated指的是差分的逆过程) ARIMA模型的表示如下： \[\begin{equation} y'_{t} = c + \phi_{1}y'_{t-1} + \cdots + \phi_{p}y'_{t-p} + \theta_{1}\varepsilon_{t-1} + \cdots + \theta_{q}\varepsilon_{t-q} + \varepsilon_{t}, \tag{8.1} \end{equation}\] 上式中 \(y'_{t}\) 是差分序列（它可能经过多次差分）。右侧的“预测变量”包括 \(y_t\) 的延迟值和延迟的误差。我们将这个模型称为ARIMA(\(p,d,q\)) 模型，参数含义如下：
| \(p=\) | 自回归模型阶数 |
| \(d=\) | 差分阶数 |
| \(q=\) | 移动平均模型阶数 |
自回归模型和于动平均模型中的平稳性和可逆性条件在ARIMA模型中依然适用。
我们之前讨论的很多模型其实都是ARIMA模型的特殊情况，如下表所示。
| 白噪声 | ARIMA(0,0,0) |
| 随机游走模型 | ARIMA(0,1,0) with no constant |
| 带漂移的随机游走模型 | ARIMA(0,1,0) with a constant |
| 自回归模型 | ARIMA(\(p\),0,0) |
| 移动平均模型 | ARIMA(0,0,\(q\)) |
#> # A tibble: 3 × 2
#> col1 col2
#> <chr> <chr>
#> 1 $p=$ 自回归模型阶数
#> 2 $d=$ 差分阶数
#> 3 $q=$ 移动平均模型阶数一旦我们开始组合不同模型来形成复杂的模型，延迟算子就会显得格外简便。比如,等式(8.1)可以被表示为： \[\begin{equation} \tag{8.2} \begin{array}{c c c c} (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\\ {\uparrow} & {\uparrow} & &{\uparrow}\\ \text{AR($p$)} & \text{$d$ differences} & & \text{MA($q$)}\\ \end{array} \end{equation}\]
R中的参数与上式稍有不同： \[\begin{equation} \tag{8.3} (1-\phi_1B - \cdots - \phi_p B^p)(y_t' - \mu) = (1 + \theta_1 B + \cdots + \theta_q B^q)\varepsilon_t, \end{equation}\] 此处 \(y_t' = (1-B)^d y_t\),\(\mu\) 是 \(y_t'\) 的均值。 为了表示成(8.2)中给出的形式，设\(c = \mu(1-\phi_1 - \cdots - \phi_p )\).
选择合适的 \(p\),\(d\) 和 \(q\) 是一个困难的过程。然而，R中的auto.arima()函数将会为你自动完成这一工作。在8.7节，我们将学习这些函数如何工作，还会介绍一些可供你自己选择取值的方法。
美国消费支出
图8.7中显示的是美国消费支出季度占比的变化情况。尽管这是一个季度时间序列，该数据并没有显示出季节性的特征，因此我们使用非季节性ARIMA模型来进行拟合。
autoplot(uschange[,"Consumption"]) +
xlab("年份") + ylab("季度占比变化")+
theme(text = element_text(family = "STHeiti"))+
theme(plot.title = element_text(hjust = 0.5))
图 8.7: 美国消费支出的季度占比变化数据.
下列R代码可以用来自动选择模型：
fit <- auto.arima(uschange[,"Consumption"], seasonal = FALSE)#> Series: uschange[, "Consumption"]
#> ARIMA(1,0,3) with non-zero mean
#>
#> Coefficients:
#> ar1 ma1 ma2 ma3 mean
#> 0.589 -0.353 0.085 0.174 0.745
#> s.e. 0.154 0.166 0.082 0.084 0.093
#>
#> sigma^2 = 0.35: log likelihood = -164.8
#> AIC=341.6 AICc=342.1 BIC=361这是一个ARIMA(1,0,3)模型： \[ y_t = c + 0.589y_{t-1} -0.353 \varepsilon_{t-1} + 0.0846 \varepsilon_{t-2} + 0.174 \varepsilon_{t-3} + \varepsilon_{t}, \] 上式中,\(c= 0.745 \times (1 - 0.589 + NA) = 0.307\),\(\varepsilon_t\) 是标准差为 \(0.592 = \sqrt{0.350}\) 的白噪声。 该模型的预测如图8.8所示。
fit %>% forecast(h=10) %>% autoplot(include=80)
图 8.8: 美国消费支出季度变化百分比预测
理解ARIMA模型
auto.arima()函数十分有用，但是完全依赖自动化的程序则是一种危险的行为。因此即使你已经使用了程序来自动选择模型，了解模型的特征和行为仍然是十分必要的。
常数 \(c\) 对这些模型的长期预测有重要的影响。
- 当 \(c=0\) 且 \(d=0\) 时，长期预测将会变为0。
- 当 \(c=0\) 且 \(d=1\) 时，长期预测将会变为一个非零常数。
- 当 \(c=0\) 且 \(d=2\) 时，长期预测将会变为直线。
- 当 \(c\ne0\) 且 \(d=0\) 时，长期预测将会变为数据的均值。
- 当 \(c\ne0\) 且 \(d=1\) 时，长期预测将会变为直线。
- 当 \(c\ne0\) 且 \(d=2\) 时，长期预测将会存在一个二次的趋势。
除此之外,\(d\) 的取值也会影响预测区间——\(d\)取值越大，预测区间(prediction interval)将会增加的更快。当 \(d=0\) 时，长期预测的标准差将会与历史数据的标准差相同，因此预测区间将会一直不变。
这种情况在图8.8中 \(d=0\) 且 \(c\ne0\) 中得以出现。在图中，预测区间在最后的几个预测范围（prediction horizon）中几乎一致，并且点预测（point forecast）等于数据的均值。
\(p\) 值的取值在数据存在周期时非常重要。为了进行周期性的预测, \(p\ge2\) 是必要条件之一，除此之外还应对系数加以额外限制。对于AR(2)模型而言，\(\phi_1^2+4\phi_2<0\)时模型将会出现周期性的特征。在这种情况下，周期的平均长度为14 \[ \frac{2\pi}{\text{arc cos}(-\phi_1(1-\phi_2)/(4\phi_2))}. \]
自相关图和偏自相关图
仅从时序图中我们通常无法判断 \(p\) 和 \(q\) 的合适取值。然而我们有时从自相关图和与其紧密相关的偏自相关图中，可以判断 \(p\) 和 \(q\) 的合适取值。
回顾之前我们所说过的，自相关图（ACF Plot）反映了自回归中 \(y_t\) 和 \(y_{t-k}\) 在不同 \(k\) 取值之下的关系。现在如果 \(y_t\) 和 \(y_{t-1}\) 已经存在相关性，则 \(y_{t-1}\) 和 \(y_{t-2}\) 必然存在相关性。然而 \(y_t\) 和 \(y_{t-2}\) 也同样必然相关，这是因为它们都与 \(y_{t-1}\) 相关而不是因为 \(y_{t-2}\) 包含新的信息可以用来预测\(y_t\)。
为了解决这个问题，我们可以使用偏自相关(partial autocorrelations)。它衡量的是在移除延迟 \(1,2,3,\dots,k - 1\) 对 \(y_t\)的影响的情况下, \(y_{t}\) 和 \(y_{t-k}\) 的关系。因此延迟一阶偏自相关系数和延迟一阶自相关系数是相同的，因为没有延迟需要移除。每个偏自相关系数都可以被估计为一个自回归模型中的末项系数。特别地，\(\alpha_k\) 作为第 \(k\) 个偏自相关系数，等于在一个AR(\(k\))模型中 \(\phi_k\) 的估计值。在实际应用中,除了拟合这些自回归，也有一些更高效的算法可以计算出 \(\alpha_k\) 的取值，但他们的结果其实是一样的。
图8.9和图8.10中显示的是图8.7中美国消费数据的自相关图和偏自相关图。偏自相关系数具有和自相关系数相同的临界值\(\pm 1.96/\sqrt{T}\)，这些在图8.9中可以看到。
ggAcf(uschange[,"Consumption"],main="")
#> Warning: Ignoring unknown parameters: main
图 8.9: 美国消费的季度变化百分比的自相关系数。
ggPacf(uschange[,"Consumption"],main="")
#> Warning: Ignoring unknown parameters: main
图 8.10: 美国消费的季度变化百分比的偏自相关系数。
如果数据来自于ARIMA(\(p\),\(d\),0)或者ARIMA(0,\(d\),\(q\))模型，则自相关图和偏自相关图在判定\(p\)或者\(q\)的取值时非常有用。15如果 \(p\) 和 \(q\) 都为正，则这些图在寻找最合适的\(p\)和\(q\)值时不再有用。
如果差分数据的自相关图和偏自相关图显示出如下特征，则他们可能来自于ARIMA(\(p\),\(d\),0)模型：
- 自相关系数呈现指数下降或者类似正弦型的波动；
- 偏自相关图中的延迟\(p\)中有明显突起，但延迟更大时不存在类似的突起。
如果差分数据的自相关图和偏自相关图数据显示出如下特征,则他们可能来自于ARIMA(0,\(d\),\(q\))模型：
- 偏自相关系数呈现指数下降或者类似正弦型的波动；
- 在自相关图中的延迟\(q\)中有明显突起,但延迟更大时不存在类似的突起。
在图8.9中，我们发现自相关图中出现了三个明显突起，之后在延迟为4时突起最为明显。在偏自相关图中，同样有三个明显的突起，之后没有任何明显突起。（除了延迟为22时出现的突起) 如果有明显突起但却不在前几个延迟期内，我们可以忽略它们。毕竟产生明显突起的概率接近二十分之一，而我们在每个图中都画出了22个突起。图中前几个突起符合ARIMA(3,0,0)的特征，偏自相关系数也逐渐下降。因此在这种情况下，拟合ARIMA(3,0,0)模型较合适。
(fit2 <- Arima(uschange[,"Consumption"],order=c(3,0,0)))
#> Series: uschange[, "Consumption"]
#> ARIMA(3,0,0) with non-zero mean
#>
#> Coefficients:
#> ar1 ar2 ar3 mean
#> 0.227 0.160 0.203 0.745
#> s.e. 0.071 0.072 0.071 0.103
#>
#> sigma^2 = 0.349: log likelihood = -165.2
#> AIC=340.3 AICc=340.7 BIC=356.5我们通过上述方法建立的模型其实拟合效果要好于使用auto.arima()建立的模型。（通过比较该模型的AIC值340.67和原模型的342.08） R中的auto.arima()函数没有找到我们建立的模型，这是因为它并没有遍历所有可能的模型。你可以通过stepwise=FALSE和approximation=FALSE等参数来增加它遍历的数量：
(fit3 <- auto.arima(uschange[,"Consumption"],seasonal=FALSE,
stepwise=FALSE,approximation=FALSE))
#> Series: uschange[, "Consumption"]
#> ARIMA(3,0,0) with non-zero mean
#>
#> Coefficients:
#> ar1 ar2 ar3 mean
#> 0.227 0.160 0.203 0.745
#> s.e. 0.071 0.072 0.071 0.103
#>
#> sigma^2 = 0.349: log likelihood = -165.2
#> AIC=340.3 AICc=340.7 BIC=356.5我们也可以设置参数seasonal=FALSE来防止程序在季节性ARIMA模型中进行选择；我们将在8.9节中探讨这些模型。
这一次，auto.arima()找到了我们通过观察自相关图和偏自相关图所得到的模型。这个ARIMA(3,0,0)的模型预测和图8.8中通过ARIMA(2,0,2)模型进行预测的结果类似，因此在这里我们不重复画图了。