11.1 复杂季节性

到目前为止,我们已经考虑了相对简单的季节模式,如季度和月度数据。然而,高频时间序列数据往往表现出更为复杂的季节模式。例如,日数据可能同时存在周季节模式和年季节模式。小时数据通常包含三种类型的季节性:日季节模式、周季节模式和年季节模式。即使是周数据也经常会因存在周期为 \(365.25/7\approx 52.179\) 的年季节模式而导致预测困难。

这种多季节模式在记录的高频数据中越来越常见。例如,呼叫中心的呼叫次数、每日住院人数、自动取款机的现金提取金额、用电用水量以及网站访问量等数据中都可能存在多季节模式。

到目前为止,我们所考虑的大多数方法都无法处理这些复杂季节性问题。即使 R 中的ts方法也只能处理一种类型的季节性,且通常假设频率取整数值。

为了处理这样的时间序列数据,我们将会使用msts方法来处理复杂季节性时间序列。这一方法允许我们指定所有可能存在的季节模式,同时也能灵活地处理频率为非整数的情况。

尽管这种方法具有灵活性,但我们没有必要将所有可能相关的频率都包含在内——而是可以仅仅包含数据中可能出现的那些频率。例如,如果你只有180天的日数据,你可以忽略年度季节性。如果你的数据是对自然现象(例如,温度)的测量,你或许可以忽略周季节性。

11.1的上部分展示了在33周内,一家零售银行在每个工作日早上7点到晚上9点之间每5分钟的接收来电次数;图的下部分则详细展示了这一时间序列前三周的数据情况。不难看出,该时间序列数据存在一个较明显的频率为169的日季节模式(每天有169个长度为5分钟的区间),和一个较弱的频率为 \(169 \times 5=845\) 的周季节模式。(周一的来电次数往往高于本周其他时间。)如果该时间序列能够有更长跨度的数据,可能还会存在年度季节模式。

p1 <- autoplot(calls) +
  ylab("来电次数") + xlab("周") +
  scale_x_continuous(breaks=seq(1,33,by=2)) +
  theme(text = element_text(family = "STHeiti"))
p2 <- autoplot(window(calls, end=4)) +
  ylab("来电次数") + xlab("周") +
  scale_x_continuous(minor_breaks = seq(1,4,by=0.2))+
  theme(text = element_text(family = "STHeiti"))
gridExtra::grid.arrange(p1,p2)
北美一家大型商业银行工作日期间每天早7点至晚9点间每5分钟的来电次数。上图展示了2003年3月3日至2003年5月23日之间的来电数据,下图仅展示了前三周的来电数据。

图 11.1: 北美一家大型商业银行工作日期间每天早7点至晚9点间每5分钟的来电次数。上图展示了2003年3月3日至2003年5月23日之间的来电数据,下图仅展示了前三周的来电数据。

包含多个季节周期的 STL

函数mstl()stl()的变体,旨在处理多重季节性。它将返回多个季节组成成分、趋势和剩余组成成分。

calls %>% mstl() %>% 
  autoplot() + xlab("周") + 
  theme(text = element_text(family = "STHeiti"))
来电次数数据的多季节 STL。

图 11.2: 来电次数数据的多季节 STL。

图示存在两个季节模式,一个是日季节模式(第三个面板),一个是周季节模式(第四个面板),需要注意的是面板中垂直方向的尺度,这对于正确地解读该图是非常重要的。在这个例子中,趋势成分和周季节成分相比其他成分的变动范围相对较窄,这是因为数据中几乎没有趋势,且周季节性也较弱。

这一分解方法也可用于预测,使用季节朴素方法对包含的每个季节成分分别进行预测,在此基础上使用 ETS 方法(或其他用户指定的方法)进行季节调整并预测。函数stlf()可以自动执行此操作。

calls %>%  stlf() %>% 
  autoplot() + xlab("周") + 
   theme(text = element_text(family = "STHeiti"))
来电次数数据的多季节 STL。

图 11.3: 来电次数数据的多季节 STL。

多季节模式动态谐波回归

就像在前面几章中所做的那样,我们可以使用傅里叶项处理包含多个季节模式的时间序列数据(5.4 节和 9.5 节)。由于包含多个季节模式,我们需要对每个季节周期添加傅里叶项。在这一例子中,季节周期分别为169和845,因此傅里叶项的形式为 \[ \sin\left(\frac{2\pi kt}{169}\right), \qquad \cos\left(\frac{2\pi kt}{169}\right), \qquad \sin\left(\frac{2\pi kt}{845}\right), \qquad \text{and} \cos\left(\frac{2\pi kt}{845}\right), \] 其中 \(k=1,2,\dots\) 。函数fourier()可以为你生成这些项。

我们将拟合一个误差具有 ARMA 结构的动态谐波回归模型。通过选择各季节周期对应的傅里叶项总数使得 AICc 最小化。我们将使用一个 log 变换(lambda=0)来确保预测值和预测区间均为正数。

fit <- auto.arima(calls, seasonal=FALSE, lambda=0,
         xreg=fourier(calls, K=c(10,10)))
fit %>%
  forecast(xreg=fourier(calls, K=c(10,10), h=2*169)) %>%
  autoplot(include=5*169) +
   ylab("来电次数") + xlab("周") +
   theme(text = element_text(family = "STHeiti"))
来电次数数据的动态谐波回归预测。

图 11.4: 来电次数数据的动态谐波回归预测。

这是一个非常大的模型,其中包含40个参数:4个 ARMA 系数,20个周期为169的傅里叶项系数和16个周期为845的傅里叶项系数。周期为845的傅里叶项中有一些与周期为169的傅里叶项重复,因此这部分重复项不再进入模型。

TBATS 模型

De Livera, Hyndman, & Snyder (2011) 提出的另一种方法是将傅里叶项与指数平滑状态空间模型、Box-Cox 转换组合在一起,以一种完全自动化的方式进行预测。与任何自动建模框架一样,尽管在某些条件下这可能是一种有用的方法,但仍存在一些该模型预测效果并不理想的情况。

TBATS 模型与动态谐波回归的不同之处在于, TBATS 模型允许季节性随着时间慢慢变化,而动态谐波回归的季节模式只能周期性地重复,而没有变化。TBATS 模型的一个缺点是估计速度很慢,特别是在长时间序列的情况下这一缺点尤为明显。因此,我们将考虑使用calls调用数据的子集,以节省时间。

calls %>%
  subset(start=length(calls)-2000) %>%
  tbats() -> fit2
fc2 <- forecast(fit2, h=2*169)
autoplot(fc2, include=5*169) +
  ylab("来电次数") + xlab("周") + 
  theme(text = element_text(family = "STHeiti"))
来电次数数据的 TBATS 模型预测。

图 11.5: 来电次数数据的 TBATS 模型预测。

这里的预测区间似乎太宽了 – 不幸的是,TBATS 模型经常会出现这种情况。

包含协变量的复杂季节性

TBATS 模型中不允许包含协变量,但协变量可以包含在动态谐波回归模型中。这种模型的一个常见应用是电力需求建模。

11.6展示了2014年澳大利亚维多利亚地区每半个小时的电力需求量,以及墨尔本(维多利亚最大的城市)同期的温度。

cbind("电力需求量" = elecdemand[, "Demand"],
      "温度" = elecdemand[, "Temperature"]) %>%
autoplot(elecdemand[,c("Demand","Temperature")], facet=TRUE) +
  scale_x_continuous(minor_breaks=NULL,
    breaks=2014+cumsum(c(0,31,28,31,30,31,30,31,31,30,31,30))/365,
    labels=month.abb) +
  xlab("时间") + ylab("") + 
  theme(text = element_text(family = "STHeiti"))
2014年澳大利亚维多利亚地区每半个小时的电力需求量及同期的温度。

图 11.6: 2014年澳大利亚维多利亚地区每半个小时的电力需求量及同期的温度。

根据电力需求量与温度的散点图(图11.7)可知,两者之间存在非线性关系:低温状态下电力需求量增加(由于需要供暖),高温状态下电力需求量增加(由于需要制冷)。

elecdemand %>%
  as.data.frame() %>%
  ggplot(aes(x=Temperature, y=Demand)) + geom_point() +
    xlab("温度 (摄氏度)") +
    ylab("电力需求量 (GW)") + 
  theme(text = element_text(family = "STHeiti"))
维多利亚地区每半个小时的电力需求量及维多利亚最大的城市墨尔本的同期温度。

图 11.7: 维多利亚地区每半个小时的电力需求量及维多利亚最大的城市墨尔本的同期温度。

我们将使用相对温度的分段线性函数(以18摄氏度为分段节点)和谐波回归项来拟合一个回归模型,以此体现数据中出现的日季节模式。

cooling <- pmax(elecdemand[,"Temperature"], 18)
fit <- auto.arima(elecdemand[,"Demand"],
         xreg = cbind(fourier(elecdemand, c(10,10,0)),
               heating=elecdemand[,"Temperature"],
               cooling=cooling))

使用这种模型进行预测是困难的,因为我们需要预测变量的未来值。傅里叶项的未来值很容易计算,但未来时刻的温度显然是未知的。如果我们只对未来一周的预测感兴趣,我们可以利用气象模型得到温度预报。或者,我们可以使用情景预测(4.5 节)并插入可能的温度值。在下面的例子中,我们通过重复过去两天的温度值来生成未来可能的电力需求量数值。

temps <- subset(elecdemand[,"Temperature"],
                start=NROW(elecdemand)-2*48+1)
fc <- forecast(fit, xreg=cbind(fourier(temps, c(10,10,0)),
                        heating=temps,cooling=pmax(temps,18)))
#> Warning in forecast.forecast_ARIMA(fit, xreg =
#> cbind(fourier(temps, c(10, : xreg contains different column
#> names from the xreg used in training. Please check that the
#> regressors are in the same order.
autoplot(fc, include=14*48) +
  xlab("时间") + ylab("电力需求量") +
  theme(text = element_text(family = "STHeiti"))
每半个小时电力需求量数据的动态谐波回归预测。

图 11.8: 每半个小时电力需求量数据的动态谐波回归预测。

尽管短期预测看起来合理,但对一个复杂过程来说,这是一个非常粗糙的模型。残差表明数据中有大量信息没有被该模型捕捉到。

checkresiduals(fc)
动态谐波回归模型的残差诊断。

图 11.9: 动态谐波回归模型的残差诊断。

#> 
#>  Ljung-Box test
#> 
#> data:  Residuals from Regression with ARIMA(5,1,4) errors
#> Q* = 738412, df = 3455, p-value <2e-16
#> 
#> Model df: 49.   Total lags used: 3504

Hyndman & Fan (2010)Fan & Hyndman (2012) 中给出了可以提供更好预测效果的更复杂模型的相关描述。

参考文献

De Livera, A. M., Hyndman, R. J., & Snyder, R. D. (2011). Forecasting time series with complex seasonal patterns using exponential smoothing. J American Statistical Association, 106(496), 1513–1527. https://robjhyndman.com/publications/complex-seasonality/

Fan, S., & Hyndman, R. J. (2012). Short-term load forecasting based on a semi-parametric additive model. IEEE Transactions on Power Systems, 27(1), 134–141. https://robjhyndman.com/publications/stlf/

Hyndman, R. J., & Fan, S. (2010). Density forecasting for long-term peak electricity demand. IEEE Transactions on Power Systems, 25(2), 1142–1153. https://robjhyndman.com/publications/peak-electricity-demand/