8.9 季节性ARIMA模型

之前我们一直将注意力集中在非季节性的数据和非季节性的 ARIMA 模型上。然而,ARIMA 模型同样可以用于很多季节性数据的建模。

季节性的 ARIMA模型 在我们之前讨论的 ARIMA 模型多项式中引入了季节性的项。它可以表示如下:

ARIMA \(\underbrace{(p, d, q)}\) \(\underbrace{(P, D, Q)_{m}}\)
\({\uparrow}\) \({\uparrow}\)
模型中的 模型中的
非季节部分 季节部分

这里的 \(m =\) 每年的观测数量。我们用大写字符来标记模型中季节性的部分,用小写字符来标记非季节性的部分。

模型的季节性部分中含有和非季节性非常类似的项,但包含了季节性时段的回溯。比如,一个ARIMA(1,1,1)(1,1,1)\(_{4}\)模型(不含常数)对应的是季节数据(\(m=4\)),它可以被表示为: \[ (1 - \phi_{1}B)~(1 - \Phi_{1}B^{4}) (1 - B) (1 - B^{4})y_{t} = (1 + \theta_{1}B)~ (1 + \Theta_{1}B^{4})\varepsilon_{t}. \] 季节性项与非季节性项相乘以后得到上述式子。

自相关图/偏自相关图

AR 或者 MA 模型中的季节性特征其实可以从偏自相关图和自相关图中的季节性延迟中被观察到。比如,一个ARIMA(0,0,0)(0,0,1)\(_{12}\)模型会显示出:

  • 自相关图中延迟期数为 12 上有一个明显突起,但是没有其他明显的突起;
  • 偏自相关图中季节性延迟(即当延迟为12,24,36,…)上出现指数性的衰减.

同样的,一个ARIMA(0,0,0)(1,0,0)\(_{12}\)

  • 自相关图中季节性延迟上出现指数型衰减;
  • 偏自相关图中延迟为 12 上的一个明显突起。

在考虑季节性 ARIMA 模型的合适的季节阶数时,我们需要对季节性延迟做出约束。

除了选择季节性的 AR 和 MA 项以及其他非季节性的项之外,这个建模过程其实和对非季节性数据是几乎一致的。我们可以通过下面的例子更好得理解这个过程。

实例:欧洲季度零售贸易额

我们这里使用欧洲1996年到2011年的季度零售贸易额来描述季节性 ARIMA 模型的建模过程。贸易额数据如图8.16所示。

autoplot(euretail) + ylab("零售指数") + xlab("年份")
欧元区(17国)季度零售贸易指数,1996--2011,包括批发和零售以及机动车修理。(指数: 2005年指数 = 100)。

图 8.16: 欧元区(17国)季度零售贸易指数,1996–2011,包括批发和零售以及机动车修理。(指数: 2005年指数 = 100)。

这里的数据显然是非平稳的,并且拥有一些季节性的特征,因此我们先对进行季节性差分,季节性差分数据如图8.17所示。由于非平稳性仍然存在,我们继续进行一阶差分,结果如图8.18所示。

euretail %>% diff(lag=4) %>% ggtsdisplay()
季节性差分后的欧洲零售贸易指数.

图 8.17: 季节性差分后的欧洲零售贸易指数.

euretail %>% diff(lag=4) %>% diff() %>% ggtsdisplay()
两次差分后的欧洲零售贸易指数.

图 8.18: 两次差分后的欧洲零售贸易指数.

经过了差分变换,我们现在的目标是通过图8.18中的自相关图和偏自相关图来找到合适的 ARIMA 模型。ACF图中在延迟为1时的突起显示出数据包含非季节性 MA(1) 的成分,而延迟为4的突起则显示出季节性 MA(1) 的成分。因此,我们先从ARIMA(0,1,1)(0,1,1)\(_4\)模型开始,这表明数据经过了一次季节性差分和一次一阶差分,同时包含季节性和非季节性的 MA(1) 成分。图8.19中显示的是拟合模型的残差。(通过类似的逻辑,我们也可以根据偏自相关图中的信息建立一个ARIMA(1,1,0)(1,1,0)\(_4\)模型。)

euretail %>%
  Arima(order=c(0,1,1),seasonal=c(0,1,1)) %>%
  residuals() %>%
  ggtsdisplay()
通过欧洲零售贸易指数数据拟合出的ARIMA(0,1,1)(0,1,1)$_4$模型的残差。

图 8.19: 通过欧洲零售贸易指数数据拟合出的ARIMA(0,1,1)(0,1,1)\(_4\)模型的残差。

自相关图和偏自相关图都在延迟为2的地方出现了明显的突起,在延迟为3的地方出现了较为明显的突起,这反映出有一些额外的非季节性项需要被补充进入模型。ARIMA(0,1,2)(0,1,1)\(_4\) 模型的AICc值 为 74.36,而 ARIMA(0,1,3)(0,1,1)\(_4\) 模型的AICc值为 68.53。我们尝试了其他拥有 AR 项的模型,但是没有找到更小的AICc值。因此,我们最后选择了 ARIMA(0,1,3)(0,1,1)\(_4\) 模型。它的残差如图8.20所示。图中所有残差都在显著性临界值内,因此它们类似于白噪声。Ljung-Box 检验也显示残差之间不存在自相关性。

fit3 <- Arima(euretail,order=c(0,1,3),seasonal=c(0,1,1))
checkresiduals(fit3)
欧洲零售指数拟合出的 ARIMA(0,1,3)(0,1,1)$_4$ 模型的残差图。

图 8.20: 欧洲零售指数拟合出的 ARIMA(0,1,3)(0,1,1)\(_4\) 模型的残差图。

#> 
#>  Ljung-Box test
#> 
#> data:  Residuals from ARIMA(0,1,3)(0,1,1)[4]
#> Q* = 0.51, df = 4, p-value = 1
#> 
#> Model df: 4.   Total lags used: 8

于是,我们现在得到了通过各项检验的一个季节性的ARIMA模型,现在我们将其运用于预测。对未来三年的预测结果如图8.21所示。预测值延续了数据的当前趋势,这是二次差分所导致的。巨大且快速增加的预测区间显示出零售贸易指数在任何时段增长和减少的可能性都较大——虽然点预测是向下的,但是预测时段的数据呈现持续增长的情况依然在预测区间之内。

fit3 %>% forecast(h=12) %>% autoplot()
ARIMA(0,1,3)(0,1,1)\(_4\) 模型对于欧洲零售贸易指数数据的预测,这里显示了 80% 和 95% 的预测区间。

图 8.21: ARIMA(0,1,3)(0,1,1)\(_4\) 模型对于欧洲零售贸易指数数据的预测,这里显示了 80% 和 95% 的预测区间。

我们可以用auto.arima()来完成这里大部分的工作,它的结果如下所示:

auto.arima(euretail)
#> Series: euretail 
#> ARIMA(0,1,3)(0,1,1)[4] 
#> 
#> Coefficients:
#>         ma1    ma2    ma3    sma1
#>       0.263  0.369  0.420  -0.664
#> s.e.  0.124  0.126  0.129   0.155
#> 
#> sigma^2 estimated as 0.156:  log likelihood=-28.63
#> AIC=67.26   AICc=68.39   BIC=77.65

auto.arima()函数使用nsdiffs()来确定 \(D\)(季节性差分的次数),使用ndiffs()来确定 \(d\)(普通差分的次数)。和非季节性模型一样,其他模型参数的选择(\(p,d,q\)\(Q\))都是通过最小化 AICc 来确定的。

实例:澳大利亚皮质类固醇药物销量

我们的第二个例子更为复杂,在这里我们将预测澳大利亚每月的皮质类固醇药物的销量。该类药物在解剖,治疗化学品分类目录中也被记为 H02 药物。

lh02 <- log(h02)
cbind("H02 sales (million scripts)" = h02,
      "Log H02 sales"=lh02) %>%
  autoplot(facets=TRUE) + xlab("年份") + ylab("")
澳大利亚皮质类固醇药物销量(每月每百万剂)对数数据见下图。

图 8.22: 澳大利亚皮质类固醇药物销量(每月每百万剂)对数数据见下图。

8.22显示的是从 1991 年 7 月到 2008 年 6 月的数据。图中方差随着销量水平增加而略微增加,因此我们运用一些算法来平稳定方差。

数据拥有很强的季节性,并且显然非平稳,因此我们使用季节性差分,季节性差分后的数据如图8.23所示。现在我们不清楚是否应该再做一次差分,我们决定不继续差分,当然这不是唯一的选择。

最后几个观测值和之前的数据较为不同(变化更大),这有可能是因为在一些情况下(比如先前的销售额没有按时汇报)数据会被修改。

lh02 %>% diff(lag=12) %>%
  ggtsdisplay(xlab="年份",main="季节性差分后的H02药物销量")
季节性差分后的澳大利亚皮质类固醇药物销量(每月每百万剂)。

图 8.23: 季节性差分后的澳大利亚皮质类固醇药物销量(每月每百万剂)。

在季节性差分后的数据的图中,我们发现偏自相关图在延迟为 12 和 24 的地方存在突起,然而在自相关图中并不存在季节性延迟的突起。这有可能暗示出模型应含有季节性 AR(2) 的部分。在非季节性的延迟中,偏自相关图有三个明显突起,显示出可能存在的 AR(3) 部分。自相关图并没有任何简单模型的特征,因此我们不作考虑。

通过上面的分析,我们确定了数据对应的一个可能的模型 ARIMA(3,0,0)(2,1,0)\(_{12}\),我们用这个模型以及它的一些衍生模型进行拟合,并且计算出了各自的AICc值如下表所示:

Model AICc
ARIMA(3,0,1)(0,1,2)\(_{12}\) -485.5
ARIMA(3,0,1)(1,1,1)\(_{12}\) -484.2
ARIMA(3,0,1)(0,1,1)\(_{12}\) -483.7
ARIMA(3,0,1)(2,1,0)\(_{12}\) -476.3
ARIMA(3,0,0)(2,1,0)\(_{12}\) -475.1
ARIMA(3,0,2)(2,1,0)\(_{12}\) -474.9
ARIMA(3,0,1)(1,1,0)\(_{12}\) -463.4

在这些模型中,ARIMA(3,0,1)(0,1,2)\(_{12}\)模型拥有最小的AICc值,因而是最优模型。

(fit <- Arima(h02,order=c(3,0,1),seasonal=c(0,1,2),lambda=0))
#> Series: h02 
#> ARIMA(3,0,1)(0,1,2)[12] 
#> Box Cox transformation: lambda= 0 
#> 
#> Coefficients:
#>          ar1    ar2    ar3    ma1    sma1    sma2
#>       -0.160  0.548  0.568  0.383  -0.522  -0.177
#> s.e.   0.164  0.088  0.094  0.190   0.086   0.087
#> 
#> sigma^2 estimated as 0.00428:  log likelihood=250
#> AIC=-486.1   AICc=-485.5   BIC=-463.3
checkresiduals(fit,lag=36)
ARIMA(3,0,1)(0,1,2)$_{12}$拟合H02药物月度销量数据的残差。

图 8.24: ARIMA(3,0,1)(0,1,2)\(_{12}\)拟合H02药物月度销量数据的残差。

#> 
#>  Ljung-Box test
#> 
#> data:  Residuals from ARIMA(3,0,1)(0,1,2)[12]
#> Q* = 51, df = 30, p-value = 0.01
#> 
#> Model df: 6.   Total lags used: 36

模型的残差如图8.24所示,在自相关图中有一些明显的突起,并且模型并没有通过Ljung-Box检验。模型仍然可以用于预测,只不过预测区间可能因为残差的相关性而不太准确。

接下来我们将会使用自动的 ARIMA 算法,在默认设置下,auto.arima()函数返回了 ARIMA(2,1,3)(0,1,1)\(_{12}\) 模型。然而,这个模型同样没有通过 Ljung-Box 检验。有时候找到一个通过所有检验的模型是不可能的。

测试集评价:

我们将通过包含最近两年数据的测试集对之前拟合出的模型进行比较。我们使用 1991 年 7 月到 2006 年 6 月的数据进行拟合,然后对 2006 年 7 月到 2008 年 6 月的药物销量进行预测,结果如下表所示:

表 8.2: 各个ARIMA模型拟合H01药物月度销量数据的均方根误差(RMSE).
Model RMSE
ARIMA(3,0,1)(0,1,2)\(_{12}\) 0.0622
ARIMA(3,0,1)(1,1,1)\(_{12}\) 0.0630
ARIMA(2,1,4)(0,1,1)\(_{12}\) 0.0632
ARIMA(2,1,3)(0,1,1)\(_{12}\) 0.0634
ARIMA(3,0,3)(0,1,1)\(_{12}\) 0.0637
ARIMA(2,1,5)(0,1,1)\(_{12}\) 0.0640
ARIMA(3,0,1)(0,1,1)\(_{12}\) 0.0644
ARIMA(3,0,2)(0,1,1)\(_{12}\) 0.0644
ARIMA(3,0,2)(2,1,0)\(_{12}\) 0.0645
ARIMA(3,0,1)(2,1,0)\(_{12}\) 0.0646
ARIMA(4,0,2)(0,1,1)\(_{12}\) 0.0648
ARIMA(4,0,3)(0,1,1)\(_{12}\) 0.0648
ARIMA(3,0,0)(2,1,0)\(_{12}\) 0.0661
ARIMA(3,0,1)(1,1,0)\(_{12}\) 0.0679

我们手动选择和auto.arima()选择的模型都在基于均方根误差的比较中表现较好,位列前四。

当我们使用 AICc 值来比较模型时,所有的模型应该是同阶差分的,这一点很重要。然而,当我们使用测试集来比较时,如何进行预测就显得不那么重要了——比较都是有效的。因此,在上表中,我们可以对使用进行差分的模型进行比较,而对于之前使用 AICc 值的表格中,我们只对进行了季节性差分而没有进行一阶差分的模型进行比较。

这些模型中没有一个能够通过所有的残差检验,在实际应用中,我们仍然会使用我们找到的最优模型,即使它不能通过所有的检验。

ARIMA(3,0,1)(0,1,2)\(_{12}\) 模型的预测(拥有最小的均方根误差,并且在只进行季节性差分的模型中拥有最优的AICc值)如图8.25所示。

h02 %>%
  Arima(order=c(3,0,1),seasonal=c(0,1,2),lambda=0) %>%
  forecast() %>%
  autoplot() +
    ylab("H02销量(百万剂)") + xlab("年份")
ARIMA(3,0,1)(0,1,2)$_{12}$模型对H01月度销量数据的预测。

图 8.25: ARIMA(3,0,1)(0,1,2)\(_{12}\)模型对H01月度销量数据的预测。