9.6 滞后预测变量

有些时候，预测变量对待预测变量的影响不是简单而迅速的。例如，投放广告之后一段时间才会影响销售接受，而当月的销售收入会取决于过去几个月的广告支出。类似的，当工厂的安全政策的变化会立即减少工厂事故的发生，但随着时间的推移，工人会降低对新出安全政策的重视程度。

在这些情况下，我们需要再模型汇总引入预测变量的滞后项。假如模型中只有一个预测变量，则假如该预测变量的滞后项之后，滞后模型可写为： \[ y_t = \beta_0 + \gamma_0x_t + \gamma_1 x_{t-1} + \dots + \gamma_k x_{t-k} + \eta_t, \] 其中，\(\eta_t\) 为 ARIMA 过程。并通过AICc中来选择 \(K\) 值和ARIMA过程中的 \(p\) 和 \(q\)。

示例：电视广告和投保数量

美国一家保险公司希望通过在国家电视台投放广告来增加保险的销售量（且通常为新出的保险种类）。图9.10中分别为从2002年一月到2005年四月中每月的保险销量和广告支出。

cbind("保险销量" = insurance[, "Quotes"],
      "广告支出" = insurance[, "TV.advert"]) %>%
  autoplot(facets=TRUE)+
  xlab("年份") + ylab("") +
  ggtitle("保险销量和广告支出")+
  theme(text = element_text(family = "STHeiti"))+
  theme(plot.title = element_text(hjust = 0.5))
每月的保险销量和广告支出。

图 9.10: 每月的保险销量和广告支出。

在滞后模型中，我们考虑当月的广告和之前三个月的广告对报销销量的影响。并采用不存在滞后项的模型与滞后模型进行对比。

# Lagged predictors. Test 0, 1, 2 or 3 lags.
Advert <- cbind(
  AdLag0 = insurance[,"TV.advert"],
  AdLag1 = stats::lag(insurance[,"TV.advert"],-1),
  AdLag2 = stats::lag(insurance[,"TV.advert"],-2),
  AdLag3 = stats::lag(insurance[,"TV.advert"],-3))[1:NROW(insurance),]

# Restrict data so models use same fitting period
fit1 <- auto.arima(insurance[4:40,1], xreg=Advert[4:40,1],   d=0)
fit2 <- auto.arima(insurance[4:40,1], xreg=Advert[4:40,1:2], d=0)
fit3 <- auto.arima(insurance[4:40,1], xreg=Advert[4:40,1:3], d=0)
fit4 <- auto.arima(insurance[4:40,1], xreg=Advert[4:40,1:4], d=0)

滞后我们通过各个模型的AICc值来选择预测变量的最优滞后阶数。

c(fit1$aicc,fit2$aicc,fit3$aicc,fit4$aicc)
#> [1] 68.50 60.02 62.83 65.46

由上述结果可以看出，当预测变量的滞后阶数为2时，AICc 值最小，因此最优模型中应当包含预测变量的二阶滞后。也就是说，应该采用当月和上月的广告支出对当月的保险销量进行建模。模型如下所示：

(fit <- auto.arima(insurance[,1], xreg=Advert[,1:2], d=0))
#> Series: insurance[, 1] 
#> Regression with ARIMA(3,0,0) errors 
#> 
#> Coefficients:
#>         ar1     ar2    ar3  intercept  AdLag0  AdLag1
#>       1.412  -0.932  0.359      2.039   1.256   0.162
#> s.e.  0.170   0.255  0.159      0.993   0.067   0.059
#> 
#> sigma^2 estimated as 0.217:  log likelihood=-23.89
#> AIC=61.78   AICc=65.4   BIC=73.43

模型的误差项为 AR（3）过程。模型可写为： \[ y_t = 2.04 + 1.26 x_t + 0.16 x_{t-1} + \eta_t, \] 其中，\(y_t\) 为第 \(t\) 月的保险销量，\(x_t\) 为第 \(t\) 月的广告支出， \[ \eta_t = 1.41 n_{t-1} -0.93 n_{t-2} + 0.36 n_{t-3} + \varepsilon_t, \] 其中，\(\varepsilon_t\) 为白噪声。

假若我们已经有了未来某月的广告支出，那么我们可以预测出该月的保险销量。假设未来每月广告支出都为8，那么我们可以得到以下预测：

fc8 <- forecast(fit, h=20,
  xreg=cbind(AdLag0=rep(8,20), AdLag1=c(Advert[40,1],rep(8,19))))
autoplot(fc8) + ylab("保险销量") +
  ggtitle("当广告支出为8个单位时保险销量的预测值")+
  theme(text = element_text(family = "STHeiti"))+
  theme(plot.title = element_text(hjust = 0.5))
当广告支出为8个单位时保险销量的预测值。

图 9.11: 当广告支出为8个单位时保险销量的预测值。