Forecasting practice

layout: true

---
class: title-slide
background-image: url("figs/titlebg.png")
background-position: 100% 50%
background-size: 100% 100%

.content-box-green-trans[
.pull-left-1[
![](figs/buaalogo.png)
]
.pull-right-2[
# Forecasting practice

### Yanfei Kang

]]

---
# 一个例子：耗电量预测

---

.pull-left[
# 尝试1

``` r
fit_lm1 <- tslm(Demand~Temperature, data = elecdaily)
autoplot(elecdaily[,'Demand'], series="真实值") +
  autolayer(fitted(fit_lm1), series="线性模型1") + 
  ylab("耗电量 (千兆瓦)") +
  xlab("时间") +
  theme(plot.title = element_text(hjust = 0.5))
```
]

.pull-right[
<img src="figure/eledaily-lm1-out-1.svg" width="504" style="display: block; margin: auto;" />
]

---

.pull-left[
# 尝试2

``` r
fit_lm2 <- tslm(Demand~Temperature + WorkDay, data = elecdaily)
autoplot(elecdaily[,'Demand'], series="真实值") +
  autolayer(fitted(fit_lm2), series="线性模型2") + 
  ylab("耗电量 (千兆瓦)") +
  xlab("时间") +
  theme(plot.title = element_text(hjust = 0.5))
```
]

.pull-right[
<img src="figure/eledaily-lm2-out-1.svg" width="504" style="display: block; margin: auto;" />
]

---
# 耗电量和气温

---
.pull-left[
# 尝试3

``` r
fit_lm3 <- tslm(Demand~Temperature + I(Temperature^2) + WorkDay, data = elecdaily)
autoplot(elecdaily[,'Demand'], series="真实值") +
  autolayer(fitted(fit_lm3), series="线性模型3") + 
  ylab("耗电量 (千兆瓦)") +
  xlab("时间") +
  theme(plot.title = element_text(hjust = 0.5)) 
```
]

.pull-right[
<img src="figure/eledaily-lm3-out-1.svg" width="504" style="display: block; margin: auto;" />
]

---
# 模型残差

---
class: inverse, center, middle
# 时间序列的自相关性

---
# 自相关性

---
# ARIMA模型的核心思路

ARIMA模型旨在描绘数据的自回归性（autocorrelations）。

---
# ARIMA

- 当我们将差分和自回归模型以及移动平均模型结合起来的时候，我们可以得到一个非季节性 ARIMA 模型。

- ARIMA 是 AutoRegressive Integrated Moving Averaging 的简称。

- ARIMA模型的表示如下:
`$$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},$$`
     
上式中 `$y^{'}_{t}$`是差分序列（它可能经过多次差分）。右侧的“预测变量”包括 `$y_{t}$`的延迟值和延迟的误差。我们将这个模型称为 **ARIMA(p,d,q) 模型**。

---
# SARIMA

季节性的ARIMA模型在我们之前讨论的ARIMA模型多项式中引入了季节性的项。

---
# 例： 澳洲啤酒产量

``` r
autoplot(ausbeer) + ylab("产量（兆升）") + xlab("年份")
```

---
# 拟合SARIMA

``` r
fit.beer <- auto.arima(ausbeer)
fit.beer
#> Series: ausbeer 
#> ARIMA(1,1,2)(0,1,1)[4] 
#> 
#> Coefficients:
#>          ar1      ma1     ma2     sma1
#>       0.0495  -1.0091  0.3746  -0.7434
#> s.e.  0.1959   0.1826  0.1530   0.0502
#> 
#> sigma^2 = 241.3:  log likelihood = -886.41
#> AIC=1782.82   AICc=1783.11   BIC=1799.63
```

---
# 预测

``` r
fit.beer %>% forecast() %>% autoplot()
```

---
# 检验

``` r
fit.beer %>% residuals() %>% ggtsdisplay()
```

---
# Your turn

继续考虑耗电量数据（数据集 `elecdaily`）。
使用 `auto.arima()`函数建模并对未来14天进行预测，通过模型评价（`accuracy()`）你得到什么样的结论？

---
# ARIMA

- 在 ARIMA 模型中，我们仅利用变量的实际观测值进行建模，却没有考虑到随时间的推移而产生的相关信息的影响。

- 在回归模型中，模型考虑了大量预测变量的相关信息，但并未采用ARIMA模型处理时间序列中的动态因素。

---
# 回归模型回顾

我们曾建立如下形式的回归模型：
`$$y_t = \beta_0 + \beta_1 x_{1,t} + \dots + \beta_k x_{k,t} + \varepsilon_t,$$`
其中， `$y_t$` 是 `$k$` 个预测变量（ `$x_{1,t},\dots,x_{k,t}$` ）的线性组合， `$\varepsilon_t$` 是不相关的误差项（例如，误差项为白噪声）。

---
class: inverse, center, middle
# 动态回归模型

---
# 动态回归模型

我们现在允许误差项中存在自相关。我们用 `$\eta_t$` 代替 `$\varepsilon_t$`。 `$\varepsilon_t$` 表示误差项，且利用 ARIMA 模型拟合残差 `$\eta_t$` 。例如，若残差 `$\eta_t$` 满足 ARIMA(1,1,1)，则回归模型可写为：
`\begin{align*}
  y_t &= \beta_0 + \beta_1 x_{1,t} + \dots + \beta_k x_{k,t} + \eta_t,\\
      & (1-\phi_1B)(1-B)\eta_t = (1+\theta_1B) \varepsilon_t,
\end{align*}`
其中 `$\varepsilon_t$` 为白噪声。模型中存在两个残差项：一是来自回归模型的误差，可用 `$\eta_t$` 表示；二是来自 ARIMA 模型的误差，可用 `$\varepsilon_t$` 表示。且假设来自 ARIMA 模型的误差为白噪声。

---
# 例：耗电量的预测

.pull-left[
数据集`elecdaily`中包含日度耗电量、工作日的指示变量（1表示为工作日，0表示非工作日）和当日最高温度。由于数据存在周度季节性，因此指定频率为7。

``` r
cbind("耗电量" = elecdaily[,"Demand"],
      "温度" = elecdaily[,"Temperature"]) %>%
  autoplot(facets=TRUE)+
  ylab("") + guides(colour="none")+
  theme(plot.title = element_text(hjust = 0.5))
```
]

.pull-right[
<img src="figure/electime-out-1.svg" width="504" style="display: block; margin: auto;" />
]

---
# 建模

``` r
xreg <- cbind(MaxTemp = elecdaily[, "Temperature"], 
              MaxTempSq = elecdaily[, "Temperature"]^2, 
              Workday = elecdaily[, "WorkDay"])
fit <- auto.arima(elecdaily[, "Demand"], xreg = xreg)
```

---
# 动态回归模型的误差
.scroll-output[

```
#> 
#> 	Ljung-Box test
#> 
#> data:  Residuals from Regression with ARIMA(2,1,2)(2,0,0)[7] errors
#> Q* = 28.229, df = 8, p-value = 0.0004326
#> 
#> Model df: 6.   Total lags used: 14
```

<img src="figure/unnamed-chunk-6-1.svg" width="504" style="display: block; margin: auto;" />
]

---
# Your turn

继续考虑耗电量数据（数据集 `elecdaily`）。
使用动态回归模型建模并对未来14天进行预测，并进行模型评价。

---
class: inverse, center, middle
#  动态谐波回归

---
# 傅里叶序列

- 季节性可以通过傅里叶项来捕捉：
`$$s_k(t) = \sin(\frac{2\pi kt}{m}),~c_k(t) = \cos(\frac{2\pi kt}{m}),$$`
`$$y_t = a + bt + \sum_{k=1}^K[\alpha_ks_k(t) + \beta_kc_k(t)] + \epsilon_t.$$`

- 任意周期函数都可以通过正弦余弦项来估计。

- 我们可以估计 `$K$`。

- `fourier()` 可以实现。

---
#  动态谐波回归

- 当数据存在较长季节性时，包含傅里叶项的动态回归模型通常效果更好。例如，日度数据存在长度为365的季节性，周度数据存在长度为52的季节性。

- ARIMA 模型可用于周期为12的月度数据或周期为4的季度数据。

- `auto.arima()`函数可以允许季节周期达到 `$m=350$`，但是实际上通常在季节性周期超过200时会导致内存不足。

- 对于这样的时间序列，我们更倾向于谐波回归方法，其中季节性模式使用傅里叶项建模，短期时间序列由 ARMA 误差项处理。

---
# 优点

- 它允许任何长度的季节性；

- 对于具有一个以上的季节性周期数据，可以包括不同频率的傅里叶项；

- 对于较小的 `$K$` 值（傅里叶项正弦余弦项的对数），季节性模式是平滑的（可以通过增加 `$K$` 来处理更多的波动的季节性）；

- 通过简单的 ARMA 误差可以轻松处理短期动态。

---
# 示例：澳大利亚月度人口外出就餐支出

- 我们组合了用于捕获季节性的傅里叶项和用于捕获数据中其它动态的 ARIMA 误差项。

- 数据集`auscafe`中，包含从2004年初到2016年11月，澳大利亚咖啡馆、餐馆和外卖的总支出（十亿美元），我们将预测之后的24个月。

---
# 示例：澳大利亚月度人口外出就餐支出

.scroll-output[
- 随着 `$K$` 的增加，傅里叶项会增加并将表现出更加不确定的季节性模式，需要简单的 ARIMA 模型来捕捉动态信息。

- 当 `$K=5$` 时，AICc 值取得最小值。从 `$K=4$` 到 `$K=5$`，AICc 值存在明显的跳跃，因此我们采用 `$K=5$` 的模型进行预测。

``` r
cafe04 <- window(auscafe, start=2004)
plots <- list()
for (i in seq(6)) {
  fit <- auto.arima(cafe04, xreg = fourier(cafe04, K = i), 
    seasonal = FALSE, lambda = 0)
  plots[[i]] <- autoplot(forecast(fit,xreg=fourier(cafe04, K=i, h=24))) +
    xlab(paste("K=",i,"   AICC=",round(fit$aicc,2))) + 
    ylab("") + ylim(1.5,4.7) 
}

gridExtra::grid.arrange(plots[[1]],plots[[2]],plots[[3]],
                        plots[[4]],plots[[5]],plots[[6]], nrow=3)
```

<img src="figure/eatout-1.svg" width="720" style="display: block; margin: auto;" />
]

---
# 5分钟客服中心电话接入量

``` r
autoplot(calls)
```

---
# 5分钟客服中心电话接入量

.scroll-output[

``` r
xreg <- fourier(calls, K = c(10,0))
(fit <- auto.arima(calls, xreg=xreg, seasonal=FALSE, stationary=TRUE))
#> Series: calls 
#> Regression with ARIMA(5,0,1) errors 
#> 
#> Coefficients:
#>          ar1      ar2      ar3      ar4     ar5      ma1  intercept   S1-169    C1-169   S2-169    C2-169    S3-169
#>       1.0587  -0.0512  -0.0324  -0.0045  0.0192  -0.8077   192.0786  55.2977  -79.0709  13.7401  -32.3622  -13.6468
#> s.e.  0.0095   0.0090   0.0088   0.0088  0.0069   0.0075     1.7542   0.7074    0.7068   0.3806    0.3804    0.2722
#>        C3-169   S4-169   C4-169   S5-169  C5-169  S6-169  C6-169  S7-169  C7-169  S8-169   C8-169   S9-169   C9-169
#>       -9.3187  -9.5112  -2.7921  -2.2355  2.8966  0.1704  3.3077  0.8531  0.2959  0.8592  -1.3888  -0.9797  -0.3421
#> s.e.   0.2721   0.2213   0.2212   0.1931  0.1931  0.1760  0.1760  0.1649  0.1649  0.1574   0.1574   0.1521   0.1521
#>       S10-169  C10-169
#>       -1.1846   0.8039
#> s.e.   0.1483   0.1482
#> 
#> sigma^2 = 242.5:  log likelihood = -115406.9
#> AIC=230869.9   AICc=230869.9   BIC=231100.3
```
]

---
# 5分钟客服中心电话接入量

.scroll-output[

``` r
fc <- forecast(fit, xreg = fourier(calls, c(10,0), 1690))
autoplot(fc)
```

<img src="figure/unnamed-chunk-9-1.svg" width="1080" style="display: block; margin: auto;" />
]

---
# Your turn

1. 在[这里](https://otexts.com/fpp2/extrafiles/retail.xlsx)下载澳大利亚月度零售数据，这些数据记录了澳大利亚各州不同类别的零售额。

2. 采用傅里叶项（Fourier terms）处理季节性成分，构建一个合适的动态回归模型，使用 AIC准则（赤池信息准则） 选择模型中需包含的傅里叶项数量。

3. 将本模型的预测结果与此前通过其他模型（如季节性ARIMA）得到的预测进行对比。

---
background-image: url("assets/titlepage_buaa.png")
background-position: 100% 50%
background-size: 100% 100%

.font300.bold[
**Thanks!**
]
<br>
<br>
<br>
.content-box-gray.font160[
<svg aria-hidden="true" role="img" viewBox="0 0 496 512" style="height:1em;width:0.97em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:#1da1f2;overflow:visible;position:relative;"><path d="M165.9 397.4c0 2-2.3 3.6-5.2 3.6-3.3.3-5.6-1.3-5.6-3.6 0-2 2.3-3.6 5.2-3.6 3-.3 5.6 1.3 5.6 3.6zm-31.1-4.5c-.7 2 1.3 4.3 4.3 4.9 2.6 1 5.6 0 6.2-2s-1.3-4.3-4.3-5.2c-2.6-.7-5.5.3-6.2 2.3zm44.2-1.7c-2.9.7-4.9 2.6-4.6 4.9.3 2 2.9 3.3 5.9 2.6 2.9-.7 4.9-2.6 4.6-4.6-.3-1.9-3-3.2-5.9-2.9zM244.8 8C106.1 8 0 113.3 0 252c0 110.9 69.8 205.8 169.5 239.2 12.8 2.3 17.3-5.6 17.3-12.1 0-6.2-.3-40.4-.3-61.4 0 0-70 15-84.7-29.8 0 0-11.4-29.1-27.8-36.6 0 0-22.9-15.7 1.6-15.4 0 0 24.9 2 38.6 25.8 21.9 38.6 58.6 27.5 72.9 20.9 2.3-16 8.8-27.1 16-33.7-55.9-6.2-112.3-14.3-112.3-110.5 0-27.5 7.6-41.3 23.6-58.9-2.6-6.5-11.1-33.3 2.6-67.9 20.9-6.5 69 27 69 27 20-5.6 41.5-8.5 62.8-8.5s42.8 2.9 62.8 8.5c0 0 48.1-33.6 69-27 13.7 34.7 5.2 61.4 2.6 67.9 16 17.7 25.8 31.5 25.8 58.9 0 96.5-58.9 104.2-114.8 110.5 9.2 7.9 17 22.9 17 46.4 0 33.7-.3 75.4-.3 83.6 0 6.5 4.6 14.4 17.3 12.1C428.2 457.8 496 362.9 496 252 496 113.3 383.5 8 244.8 8zM97.2 352.9c-1.3 1-1 3.3.7 5.2 1.6 1.6 3.9 2.3 5.2 1 1.3-1 1-3.3-.7-5.2-1.6-1.6-3.9-2.3-5.2-1zm-10.8-8.1c-.7 1.3.3 2.9 2.3 3.9 1.6 1 3.6.7 4.3-.7.7-1.3-.3-2.9-2.3-3.9-2-.6-3.6-.3-4.3.7zm32.4 35.6c-1.6 1.3-1 4.3 1.3 6.2 2.3 2.3 5.2 2.6 6.5 1 1.3-1.3.7-4.3-1.3-6.2-2.2-2.3-5.2-2.6-6.5-1zm-11.4-14.7c-1.6 1-1.6 3.6 0 5.9 1.6 2.3 4.3 3.3 5.6 2.3 1.6-1.3 1.6-3.9 0-6.2-1.4-2.3-4-3.3-5.6-2z"/></svg> [github.com/kl-lab](https://github.com/kl-lab)

]