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模型的特殊情况,如下表所示。

表 8.1: 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 x 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 estimated as 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 estimated as 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=FALSEapproximation=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 estimated as 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)模型进行预测的结果类似,因此在这里我们不重复画图了。


  1. arc cos是cosine函数的反函数。计算器上可以找到这个函数,它可能被表示为 acos 或者 cos\(^{-1}\)↩︎

  2. ggtsdisplay函数可以用一行命令画出时序图、自相关图和偏自相关图。↩︎