商业预测

layout: true

---
class: inverse, center, middle

background-image: url(figures/titlepage16-9.png)
background-size: cover

<br>

# 商业预测
## 第五讲：动态回归模型及应用

#### 康雁飞 | 北航经管学院 | 数量经济与商务统计系

---

# 核心思路

- 在 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 模型的误差为白噪声。

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

<div class="figure" style="text-align: center">
<img src="BF-L5-dynamic_files/figure-html/elecscatter-1.png" alt="2014年澳大利亚维多利亚州每日耗电量与每日最高温度的散点图。"  />
<p class="caption">2014年澳大利亚维多利亚州每日耗电量与每日最高温度的散点图。</p>
</div>

---

.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="BF-L5-dynamic_files/figure-html/electime-out-1.png" 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)
```

---
# 动态回归模型的误差

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

---

.pull-left[
# 预测

```r
fcast <- forecast(fit, 
  xreg = cbind(rep(26,14), rep(26^2,14), c(0,1,0,0,1,1,1,1,1,0,0,1,1,1)))
autoplot(fcast) + ylab("日度耗电量（千兆瓦）") +
  theme(plot.title = element_text(hjust = 0.5))
```
]

.pull-right[
<img src="BF-L5-dynamic_files/figure-html/elecdailyf-out-1.png" style="display: block; margin: auto;" />

]

---
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 误差可以轻松处理短期动态。

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

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

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

```r
cafe04 <- window(auscafe, start=2004)
autoplot(cafe04)
```

<img src="BF-L5-dynamic_files/figure-html/unnamed-chunk-2-1.png" style="display: block; margin: auto;" />
]

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

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

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

```r
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="BF-L5-dynamic_files/figure-html/eatout-1.png" style="display: block; margin: auto;" />
]

---
# 示例：周度美国成品汽油产品供应量

.scroll-output[

```r
harmonics <- fourier(gasoline, K = 13)
(fit <- auto.arima(gasoline, xreg = harmonics, seasonal = FALSE))
```

```
## Series: gasoline 
## Regression with ARIMA(0,1,2) errors 
## 
## Coefficients:
##           ma1     ma2   drift   S1-52    C1-52    S2-52    C2-52   S3-52
##       -0.9612  0.0936  0.0014  0.0315  -0.2555  -0.0522  -0.0175  0.0242
## s.e.   0.0275  0.0286  0.0008  0.0124   0.0124   0.0090   0.0089  0.0082
##         C3-52   S4-52    C4-52    S5-52    C5-52   S6-52    C6-52   S7-52
##       -0.0989  0.0321  -0.0257  -0.0011  -0.0472  0.0580  -0.0320  0.0283
## s.e.   0.0082  0.0079   0.0079   0.0078   0.0078  0.0078   0.0078  0.0079
##        C7-52   S8-52   C8-52    S9-52   C9-52   S10-52  C10-52  S11-52   C11-52
##       0.0368  0.0238  0.0139  -0.0172  0.0119  -0.0236  0.0230  0.0001  -0.0191
## s.e.  0.0079  0.0079  0.0079   0.0080  0.0080   0.0081  0.0081  0.0082   0.0082
##        S12-52   C12-52  S13-52   C13-52
##       -0.0288  -0.0177  0.0012  -0.0176
## s.e.   0.0083   0.0083  0.0084   0.0084
## 
## sigma^2 = 0.05603:  log likelihood = 43.66
## AIC=-27.33   AICc=-25.92   BIC=129
```
]

---
# 示例：周度美国成品汽油产品供应量

```r
checkresiduals(fit, test=FALSE)
```

---
# 示例：周度美国成品汽油产品供应量

.scroll-output[

```r
newharmonics <- fourier(gasoline, K = 13, h = 156)
fc <- forecast(fit, xreg = newharmonics)
autoplot(fc)
```

<img src="BF-L5-dynamic_files/figure-html/unnamed-chunk-5-1.png" 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
##       1.0587  -0.0512  -0.0324  -0.0045  0.0192  -0.8077   192.0786  55.2977
## s.e.  0.0095   0.0090   0.0088   0.0088  0.0069   0.0075     1.7542   0.7074
##         C1-169   S2-169    C2-169    S3-169   C3-169   S4-169   C4-169   S5-169
##       -79.0709  13.7401  -32.3622  -13.6468  -9.3187  -9.5112  -2.7921  -2.2355
## s.e.    0.7068   0.3806    0.3804    0.2722   0.2721   0.2213   0.2212   0.1931
##       C5-169  S6-169  C6-169  S7-169  C7-169  S8-169   C8-169   S9-169   C9-169
##       2.8966  0.1704  3.3077  0.8531  0.2959  0.8592  -1.3888  -0.9797  -0.3421
## s.e.  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="BF-L5-dynamic_files/figure-html/unnamed-chunk-8-1.png" style="display: block; margin: auto;" />
]

---
class: inverse, center, middle
# 滞后预测变量

---
# 滞后预测变量

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

在这些情况下，我们需要再模型汇总引入预测变量的滞后项。假如模型中只有一个预测变量，则假如该预测变量的滞后项之后，滞后模型可写为：
`$$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$`。

---
# 例：电视广告和投保数量

.pull-left[
美国一家保险公司希望通过在国家电视台投放广告来增加保险的销售量（且通常为新出的保险种类）。

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

.pull-right[
<img src="BF-L5-dynamic_files/figure-html/tvadvert-out-1.png" style="display: block; margin: auto;" />
]

---
# 滞后模型

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

```r
# 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值来选择预测变量的最优滞后阶数。

```r
c(fit1$aicc,fit2$aicc,fit3$aicc,fit4$aicc)
```

```
## [1] 68.49968 60.02357 62.83253 65.45747
```

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

```r
(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.4117  -0.9317  0.3591     2.0393  1.2564  0.1625
## s.e.  0.1698   0.2545  0.1592     0.9931  0.0667  0.0591
## 
## sigma^2 = 0.2165:  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 \eta_{t-1}-0.93 \eta_{t-2}+0.36 \eta_{t-3}+\varepsilon_{t}$$`
其中， `$\varepsilon_t$` 为白噪声。

---
# 保险销量预测

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

```r
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(plot.title = element_text(hjust = 0.5))
```
]

.pull-right[
<img src="BF-L5-dynamic_files/figure-html/tvadvertf8-out-1.png" style="display: block; margin: auto;" />
]

---
# 实验5

数据集`advert`中包含汽车零部件的月度销售额和广告支出。
1. 使用`autoplot`函数绘制上述数据。参数`facets=TRUE`可以起到什么样的作用？
2. 利用`tslm()`函数拟合标准回归模型 `$y_t = a + b x_t + \eta_t$`，其中 `$y_t$` 代表销售额， `$x_t$` 代表广告支出。
3. 证明残差具有显著的自相关性。
4. 使用`auto.arima()`函数重新拟合。误差模型的参数估计有什么差异？应选用哪种 ARIMA 模型拟合残差？
5. 检验拟合模型的残差。
6. 假设未来六个月的广告支出为每月十个单位，通过模型预测月度销售额，并生成预测区间。

---
# 总结

1. 动态回归模型

2. 滞后模型

3. 适用场景