商业预测

layout: true

---
class: inverse, center, middle

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

<br>

# 商业预测
## 第六讲：商业预测前沿

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

---

# 商业预测前沿方法

- 分层或分组时间序列预测

- 向量自回归模型

- 神经网络模型

- 预测组合

- Bagging

- 预测现实

---
class: inverse, center, middle
# 分层或分组时间序列预测

---
# 分层时间序列

- 时间序列通常可以按照各种感兴趣的属性自然地进行分类。
    - 自行车制造商销售的自行车总数可以按产品类型分类，例如公路自行车、山地自行车、儿童自行车和混合动力车。
    - 其中每一个还可以进行更精细的类别。例如，混合动力自行车可分为城市、通勤、舒适和徒步自行车等等。
    - 这些类别嵌套在较大的组类别中，同样地，时间序列的集合遵循一个分层聚合的结构。
- 由于地理划分，经常出现分层时间序列。例如，自行车的总销售额可以按国家分类，然后在国家按州分类，每个州按地区分类，依此类推，直至出口水平。
- 通常可以根据分类的时间序列生成分类预测，我们一般要求预测是以和数据相同的方式进行累加。例如，应该将区域销售的预测加起来给出对州销售的预测，而各州的预测又应该加起来给出全国销售的预测。

---
# 例：澳大利亚旅游层次

- 澳大利亚分为八个地理区域（一些称为州，另一些称为地区），每个区域都有自己的政府和一些经济和行政自治机构。这些中的每一个可以按兴趣进一步细分为更小的部分，称为区。

- 商业规划者和旅游部门对整个澳大利亚，各州和地区以及各区的预测很感兴趣。在这个例子中，我们主要关注国内旅游的季度需求，以澳大利亚人出门在外的过夜数来衡量。

- 有六个州：新南威尔士州（NSW），昆士兰州（QLD），南澳大利亚州（SAU），维多利亚州（VIC），西澳大利亚州（WAU）和其他州（OTH）。对于其中的每一个州，我们考虑以下这些区的游客过夜数。

---
# 例：澳大利亚旅游层次

|州     |区                                   |
|:------|:---------------------------------------------------------------------|
|NSW    |大都会（NSWMetro），北海岸（NSWNthCo），南海岸（NSWSthCo），南部内陆（NSWSthIn），北部内陆（NSWNthIn）|
|QLD    |大都会（QLDMetro），中部（QLDCntrl），北海岸（QLDNthCo）|
|SAU    |大都会（SAUMetro），海岸（SAUCoast），内陆（SAUInner）|
|VIC    |大都会（VICMetro），西海岸（VICWstCo），东海岸（VICEstCo），内陆（VICInner）|
|WAU    |大都会（WAUMetro），海岸（WAUCoast），内陆（WAUInner）|
|OTH    |大都会（OTHMetro），非大都会（OTHNoMet）|

---
## `hts()`

.pull-left[

```r
library(hts)
tourism.hts <- hts(visnights, characters = c(3, 5))
tourism.hts %>% aggts(levels=0:1) %>%
  autoplot(facet=TRUE) +
  xlab("年份") + ylab("百万") +
  ggtitle("游客过夜数")+
  theme(plot.title = element_text(hjust = 0.5))
```
]

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

---
# 分组时间序列

- 分组时间序列比分层时间序列有更通用的聚合结构。

- 对于分组的时间序列，结构不会以独特的分层方式自然地分类，并且通常分类的因素都是嵌套和交叉的。

- 例如，我们可以通过旅行目的（例如度假，公务等。）进一步对澳大利亚旅游数据的所有地理级别再分类。

- 因此，我们可以考虑按照旅行目的来对整个澳大利亚，每个州以及每个区来划分游客过夜数。然后我们将结构描述为与地理层次“交叉”的旅行目的。

---
# 例：澳大利亚监狱人口数

<div class="figure" style="text-align: center">
<img src="BF-L6-advanced_files/figure-html/prison-1.png" alt="澳大利亚成人监狱总人口季度数据，按州、法律地位、性别分类结果。"  />
<p class="caption">澳大利亚成人监狱总人口季度数据，按州、法律地位、性别分类结果。</p>
</div>

---
# 可视化

绘制主要组别的一种方法如下所示。

```r
prison.gts <- gts(prison/1e3, characters = c(3,1,9),
  gnames = c("State", "Gender", "Legal",
             "State*Gender", "State*Legal",
             "State*Gender*Legal"))
prison.gts %>% aggts(level=0:3) %>% autoplot()
```

---
# 自下而上预测

```r
forecast(tourism.hts, method="bu", fmethod="arima")
```

```
## Hierarchical Time Series 
## 3 Levels 
## Number of nodes at each level: 1 6 20 
## Total number of series: 27 
## Number of observations in each historical series: 76 
## Number of forecasts per series: 8 
## Top level series of forecasts: 
##          Qtr1     Qtr2     Qtr3     Qtr4
## 2017 90.76257 72.40181 76.00051 78.42748
## 2018 90.16914 72.78868 75.65169 78.33982
```

---
# 自上而下预测

```r
forecast(tourism.hts, method="tdfp", fmethod="arima")
```

```
## Hierarchical Time Series 
## 3 Levels 
## Number of nodes at each level: 1 6 20 
## Total number of series: 27 
## Number of observations in each historical series: 76 
## Number of forecasts per series: 8 
## Top level series of forecasts: 
##          Qtr1     Qtr2     Qtr3     Qtr4
## 2017 95.65193 78.02317 82.30592 83.67373
## 2018 97.22648 79.59772 83.88048 85.24828
```

---
# 中间突破法预测

```r
forecast(tourism.hts, method="mo", fmethod="arima", level=1)
```

```
## Hierarchical Time Series 
## 3 Levels 
## Number of nodes at each level: 1 6 20 
## Total number of series: 27 
## Number of observations in each historical series: 76 
## Number of forecasts per series: 8 
## Top level series of forecasts: 
##          Qtr1     Qtr2     Qtr3     Qtr4
## 2017 93.32627 75.70450 79.51271 81.15298
## 2018 93.90874 76.57744 80.37253 82.02744
```

---
# 最优调和法预测

```r
forecast(tourism.hts, method="comb", fmethod="arima")
```

```
## Hierarchical Time Series 
## 3 Levels 
## Number of nodes at each level: 1 6 20 
## Total number of series: 27 
## Number of observations in each historical series: 76 
## Number of forecasts per series: 8 
## Top level series of forecasts: 
##          Qtr1     Qtr2     Qtr3     Qtr4
## 2017 92.55729 74.59201 78.37077 80.33918
## 2018 92.75248 75.35934 78.79227 80.88062
```

---
# 预测效果评价

```r
train <- window(tourism.hts, end = c(2014, 4))
test <- window(tourism.hts, start = 2015)
fcsts.opt <- forecast(train, h = 8, method = "comb", fmethod = "arima")
fcsts.bu <- forecast(train, h = 8, method = "bu", fmethod = "arima")
cbind(accuracy(fcsts.opt, test, levels = 0), accuracy(fcsts.bu, test, levels = 0))
```

```
##          Total     Total
## ME    8.338903  8.521813
## RMSE  8.708361  8.894777
## MAE   8.338903  8.521813
## MAPE 10.232566 10.453860
## MPE  10.232566 10.453860
## MASE  2.552022  2.607999
```

---
# 预测效果比较

<table>
<caption></caption>
 <thead>
<tr>
<th style="empty-cells: hide;border-bottom:hidden;" colspan="1"></th>
<th style="border-bottom:hidden;padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="2"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px; ">自下而上法</div></th>
<th style="border-bottom:hidden;padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="2"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px; ">最优调和法</div></th>
</tr>
  <tr>
   <th style="text-align:left;">   </th>
   <th style="text-align:right;"> MAPE </th>
   <th style="text-align:right;"> MASE </th>
   <th style="text-align:right;"> MAPE </th>
   <th style="text-align:right;"> MASE </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> 总计 </td>
   <td style="text-align:right;"> 10.45 </td>
   <td style="text-align:right;"> 2.61 </td>
   <td style="text-align:right;"> 10.23 </td>
   <td style="text-align:right;"> 2.55 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 州 </td>
   <td style="text-align:right;"> 12.01 </td>
   <td style="text-align:right;"> 1.71 </td>
   <td style="text-align:right;"> 11.93 </td>
   <td style="text-align:right;"> 1.73 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 底层 </td>
   <td style="text-align:right;"> 13.92 </td>
   <td style="text-align:right;"> 1.29 </td>
   <td style="text-align:right;"> 13.86 </td>
   <td style="text-align:right;"> 1.29 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 所有序列 </td>
   <td style="text-align:right;"> 13.37 </td>
   <td style="text-align:right;"> 1.43 </td>
   <td style="text-align:right;"> 13.29 </td>
   <td style="text-align:right;"> 1.44 </td>
  </tr>
</tbody>
</table>

---
# 最优调和法预测
.pull-left[
### 例：预测澳大利亚监狱人口

```r
prisonfc <- forecast(prison.gts)
```

要获得每个聚合级别的预测值，我们可以使用 `aggts()` 函数。

```r
fcsts <- aggts(prisonfc, levels=0:3)
```

使用下列代码画简单的图：

```r
groups <- aggts(prison.gts, levels=0:3)
autoplot(fcsts) + autolayer(groups)
```
]

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

]

---
# 预测效果评价

<table>
<caption>澳大利亚成年监狱人口按所有属性的交互的分组的一致预测。</caption>
 <thead>
<tr>
<th style="empty-cells: hide;border-bottom:hidden;" colspan="1"></th>
<th style="border-bottom:hidden;padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="2"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px; ">自下而上法</div></th>
<th style="border-bottom:hidden;padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="2"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px; ">最优调和法</div></th>
</tr>
  <tr>
   <th style="text-align:left;">   </th>
   <th style="text-align:right;"> MAPE </th>
   <th style="text-align:right;"> MASE </th>
   <th style="text-align:right;"> MAPE </th>
   <th style="text-align:right;"> MASE </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> 总计 </td>
   <td style="text-align:right;"> 5.32 </td>
   <td style="text-align:right;"> 1.84 </td>
   <td style="text-align:right;"> 3.08 </td>
   <td style="text-align:right;"> 1.06 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 州 </td>
   <td style="text-align:right;"> 7.59 </td>
   <td style="text-align:right;"> 1.88 </td>
   <td style="text-align:right;"> 7.62 </td>
   <td style="text-align:right;"> 1.85 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 法律地位 </td>
   <td style="text-align:right;"> 6.40 </td>
   <td style="text-align:right;"> 1.76 </td>
   <td style="text-align:right;"> 4.32 </td>
   <td style="text-align:right;"> 1.14 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 性别 </td>
   <td style="text-align:right;"> 8.62 </td>
   <td style="text-align:right;"> 2.68 </td>
   <td style="text-align:right;"> 8.72 </td>
   <td style="text-align:right;"> 2.74 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 底层 </td>
   <td style="text-align:right;"> 15.82 </td>
   <td style="text-align:right;"> 2.23 </td>
   <td style="text-align:right;"> 15.25 </td>
   <td style="text-align:right;"> 2.16 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 所有序列 </td>
   <td style="text-align:right;"> 12.41 </td>
   <td style="text-align:right;"> 2.16 </td>
   <td style="text-align:right;"> 12.02 </td>
   <td style="text-align:right;"> 2.08 </td>
  </tr>
</tbody>
</table>

---
class: inverse, center, middle
# 向量自回归模型

---
# 向量自回归模型

- 到目前为止，我们考虑的所有模型都存在一个强加单向关系的局限性 —— 预测变量受预测因子的影响，但反之不成立。

- 个人消费支出的变化（ `$C_t$` ）是基于个人可支配收入的变化（ `$I_t$` ）进行预测的，但在该例中，双向影响关系可能会更合适。

- 在2008年至2009年全球金融危机期间，澳大利亚就出现了这种情况。澳大利亚政府在2008年12月推出了包括现金支付在内的一揽子刺激计划，正好赶上圣诞节支出。结果，零售商报告显示销量增加，经济也得到刺激，结果是收入也增加了。

---
# 向量自回归模型

- 向量自回归（ VAR ）框架允许这类双向反馈关系。
- 在VAR框架中，所有的变量均被平等对待。
- 在建模过程中，变量之间平等地相互影响。
- 在更正式的术语中，所有变量都被视为“内生的”。
- 为了表示这一点，我们改变记号并将所有变量写成 `$y$` ： `$y_{1,t}$` 表示变量 `$y_1$` 的第 `$t$` 个观察值， `$y_{2,t}$` 表示变量 `$y_2$` 的第 `$t$` 个观察值，依此类推。
- 我们考虑滞后一期的二变量 VAR 模型，于是得到一个二维的 VAR(1) 模型
`\begin{align}
  y_{1,t} &= c_1+\phi _{11,1}y_{1,t-1}+\phi _{12,1}y_{2,t-1}+e_{1,t} \\
  y_{2,t} &= c_2+\phi _{21,1}y_{1,t-1}+\phi _{22,1}y_{2,t-1}+e_{2,t}, 
\end{align}`
其中 `$e_{1,t}$` 和 `$e_{2,t}$` 都是白噪声过程，且可能存在同期相关。系数 `$\phi_{ii,\ell}$` 表示变量 `$y_i$` 的 `$\ell$` 期滞后项对其自身的影响，系数 `$\phi_{ij,\ell}$` 表示变量 `$y_j$` 的 `$\ell$` 期滞后项对 `$y_i$` 的影响。

---
# 例：预测美国消费的 VAR 模型

```r
library(vars)
VARselect(uschange[,1:2], lag.max=8, type="const")[["selection"]]
```

```
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      5      1      1      5
```

---
# 例：预测美国消费的 VAR 模型

R 输出结果显示了在 **vars** 包中根据每个信息准则选择得到的滞后期长度。由 AIC 选择的 VAR(5) 与 BIC 选择的 VAR(1) 之间存在较大的差异。这种情况并不罕见。因此，我们首先拟合由 BIC 选择的 VAR(1) 。

```r
var1 <- VAR(uschange[,1:2], p=1, type="const")
serial.test(var1, lags.pt=10, type="PT.asymptotic")
var2 <- VAR(uschange[,1:2], p=2, type="const")
serial.test(var2, lags.pt=10, type="PT.asymptotic")
```

---
# 例：预测美国消费的 VAR 模型

与单变量 ARIMA 方法类似，我们使用混合检验对残差的相关性进行检验。VAR(1) 和 VAR(2) 的残差均存在序列相关，因此我们拟合 VAR(3)

```r
var3 <- VAR(uschange[,1:2], p=3, type="const")
serial.test(var3, lags.pt=10, type="PT.asymptotic")
```

```
## 
## 	Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object var3
## Chi-squared = 33.617, df = 28, p-value = 0.2138
```

该模型的残差通过了序列相关的检验。

---
# 预测

```r
forecast(var3) %>%
  autoplot() + xlab("年") 
```

<div class="figure" style="text-align: center">
<img src="BF-L6-advanced_files/figure-html/VAR3-1.png" alt=" VAR(3) 对美国消费和收入的预测。"  />
<p class="caption"> VAR(3) 对美国消费和收入的预测。</p>
</div>

---
class: inverse, center, middle
# 神经网络模型

---
# 神经网络架构

.pull-left[
- 神经网络可以被理解为是一种分层的“神经元”网络结构。

- 预测变量（或输入）构成底层，响应变量（或输出）构成顶层。还可能存在包含“隐藏神经元”的中间层。

- 最简单的网络不包含中间的隐藏层，等价于线性回归。
]

.pull-right[
<div class="figure" style="text-align: center">
<img src="figures/nnet1.png" alt="与线性回归等价的简单神经网络。" width="75%" />
<p class="caption">与线性回归等价的简单神经网络。</p>
</div>
]

---
# 神经网络架构

如果我们增加一个包含隐藏神经元的中间层，神经网络就成为非线性形式。

<div class="figure" style="text-align: center">
<img src="figures/nnet2.png" alt="包含4个输入及1个隐藏层的神经网络，其中隐藏层中包含3个隐藏神经元。" width="80%" />
<p class="caption">包含4个输入及1个隐藏层的神经网络，其中隐藏层中包含3个隐藏神经元。</p>
</div>

---
# 神经网络自回归

- 类似于我们在线性自回归模型中使用滞后值作为自变量一样，对于时间序列数据，时间序列的滞后值也可以作为神经网络的输入。我们称之为神经网络自回归或 NNAR 模型。
- 我们只考虑包含一个隐藏层的前馈网络的情况，使用符号 NNAR($p,k$) 来表示隐藏层中有 `$p$` 期滞后输入和 `$k$` 个节点。
- 例如， NNAR(9,5) 模型是一个使用滞后9期观测值 `$(y_{t-1},y_{t-2},\dots,y_{t-9})$` 作为输入，预测 `$y_t$` 的神经网络，且该神经网络中隐藏层有5个神经元。 NNAR( `$p,0$` ) 模型相当于 ARIMA( `$p,0,0$`) 模型，但这里不包含为确保平稳性而对参数的限制。
- 对于季节性数据，可以考虑增加同一季节的最后观测值作为输入。例如， NNAR(3,1,2) `$_{12}$` 模型的输入为 `$y_{t-1}$` ,  `$y_{t-2}$` ,  `$y_{t-3}$` ,  `$y_{t-12}$`  ，隐藏层中有两个神经元。更一般地， NNAR( `$p,P,k$`) `$_m$` 模型的输入为  `$(y_{t-1},y_{t-2},\dots,y_{t-p},y_{t-m},y_{t-2m},y_{t-Pm})$` ，隐藏层有 `$k$` 个神经元。 NNAR( `$p,P,0$`) `$_m$` 模型相当于 ARIMA( `$p,0,0$`)( `$P$`,0,0) `$_m$` 模型，但不包含为确保平稳性而对参数的限制。

---
# 神经网络自回归

- 函数 `nnetar()`可用于拟合 NNAR( `$p,P,k$`) `$_m$` 模型。

- 使用神经网络进行预测的时候，计算会迭代进行。比如向前预测一步时，我们只需要将历史数据作为输入。向前预测两步时，需要将历史数据和向前一步预测值作为神经网络的输入。这样的迭代会继续向后进行，直到所有需要的预测值都被计算出来。

---
# 例：太阳黑子

.pull-left[
太阳的表面含有磁区，这些磁区看起来像黑色的点。它们会影响无线电波的传播，因此电信公司喜欢预测太阳黑子的活动，以便为将来的困难做好准备。太阳黑子的一个周期为9到14年。下图中展示了 NNAR(10,6) 模型未来30年的预测值。我们进行了`lambda=0`的Box-Cox 变换，以确保预测值始终为正数值。

```r
fit <- nnetar(sunspotarea, lambda=0)
autoplot(forecast(fit,h=30)) +
  xlab("时间") + ylab("太阳黑子活动") 
```
]
.pull-right[
<img src="BF-L6-advanced_files/figure-html/sunspotnnetar-out-1.png" style="display: block; margin: auto;" />
]

---
class: inverse, center, middle
# 预测组合

---
# 模型组合

把不同的模型组合在一起进行预测通常会得到比最优预测更优的结果。

- 简单组合： `$f^{c}=\frac{1}{P} \sum_{i=1}^{P} f_{i}$`。

- 基于回归的组合： `$f^{c}= \sum_{i=1}^{P} \beta_i f_{i}$`。

- 基于特征的组合： `$f^{c}= \sum_{i=1}^{P} \omega_i f_{i}$`。

---
# 模型组合

---
# 示例：澳大利亚外出就餐的月度支出数据

我们使用以下模型进行预测: ETS ， ARIMA ， STL-ETS ， NNAR 和 TBATS ；我们对过去5年(60个月)的观察结果进行比较。

.scroll-output[

```r
train <- window(auscafe, end=c(2012,9))
h <- length(auscafe) - length(train)
ETS <- forecast(ets(train), h=h)
ARIMA <- forecast(auto.arima(train, lambda=0, biasadj=TRUE), h=h)
STL <- stlf(train, lambda=0, h=h, biasadj=TRUE)
NNAR <- forecast(nnetar(train), h=h)
TBATS <- forecast(tbats(train, biasadj=TRUE), h=h)
Combination <- (ETS[["mean"]] + ARIMA[["mean"]] +
  STL[["mean"]] + NNAR[["mean"]] + TBATS[["mean"]])/5
autoplot(auscafe) +
  autolayer(ETS, series="ETS", PI=FALSE) +
  autolayer(ARIMA, series="ARIMA", PI=FALSE) +
  autolayer(STL, series="STL", PI=FALSE) +
  autolayer(NNAR, series="NNAR", PI=FALSE) +
  autolayer(TBATS, series="TBATS", PI=FALSE) +
  autolayer(Combination, series="Combination") +
  xlab("年") + ylab("十亿美元") +
  ggtitle("澳大利亚外出就餐月支出")
```

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

---
# 预测评价

```r
c(ETS=accuracy(ETS, auscafe)["Test set","RMSE"],
  ARIMA=accuracy(ARIMA, auscafe)["Test set","RMSE"],
  `STL-ETS`=accuracy(STL, auscafe)["Test set","RMSE"],
  NNAR=accuracy(NNAR, auscafe)["Test set","RMSE"],
  TBATS=accuracy(TBATS, auscafe)["Test set","RMSE"],
  Combination=accuracy(Combination, auscafe)["Test set","RMSE"])
```

```
##         ETS       ARIMA     STL-ETS        NNAR       TBATS Combination 
##  0.13699696  0.15919858  0.19310148  0.33548889  0.09406039  0.07217230
```

---
# 参考

参考预测组合综述论文:

Xiaoqian Wang, Rob Hyndman, Feng Li, Yanfei Kang (2022). Forecast combinations: an over 50-year review. [Working paper](https://arxiv.org/pdf/2205.04216.pdf).

---
class: inverse, center, middle
# Bagging

---
# Bagging
###  **B**ootstrap **agg**regat**ing**

![](https://otexts.com/fppcn/fpp_files/figure-html/stlbookts-1.png)

---
# Bagging

.pull-left[
- Bootstrap 得到的时间序列可以用来提高预测的准确性。

- 如果我们对每一个 Bootstrap 时间序列进行预测，并对预测结果进行平均，那么我们会得到相比直接对原始时间序列进行预测更好的预测结果。

- 这个过程称为"bagging"。
]

.pull-right[
![](https://otexts.com/fppcn/fpp_files/figure-html/baggedets-1.png)
]

---
class: inverse, center, middle
# 预测现实

---
class: inverse, middle

> "In theory, there is no difference between theory
and practice. But, in practice, there is."

> -- Benjamin Brewster (1882).

---
# 预测现实 1：周数据、日数据及日以下数据

周数据很难进行预测是因为季节周期（一年中的周数）非常大且非整数。一年中的平均周数是 52.18 。我们考虑的大多数方法都是基于整数季节周期进行的。即使我们将周数据的季节周期近似为52，大多数方法也难以有效地处理这么大的季节周期数据。

1.  STL 分解，并将非季节性方法应用于季节调整数据

2. 动态谐波回归模型

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

.scroll-output[

```r
gasoline %>% stlf() %>% autoplot() + 
  xlab("时间")
```

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

---
# 日数据及日以下数据

- 日数据和日以下数据难以预测的原因与周数据相比则有所不同 —— 它们通常涉及多个季节模式，所以我们需要使用一种能处理这种复杂季节性数据的方法。

- STL、动态谐波回归或者 TBATS。

- 适用于常规季节性的时间序列。捕捉像复活节、农历新年这些移动事件的季节性是十分困难的。即使有月度数据也很棘手，因为这些节日可能在三月或四月（复活节），一月或二月（农历新年），也可能在一年中的任何时候。

- 处理移动节日的最佳方法是使用虚拟变量作为协变量，此时可以使用动态回归模型/其他形式的回归模型，在此模型中，自变量中可以包含描述节日效应的虚拟变量（也可以是傅里叶项）。

---
# 预测现实 2：计数时间序列

.pull-left[

```r
productC %>% croston() %>% 
  autoplot() +
   xlab("时间") +
  ggtitle('某石油公司润滑剂销量')
```

参考 [`tsintermittent`](https://cran.r-project.org/package=tsintermittent) 包。

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

---
# 预测现实 3：历史数据有限

- 如何使用很少的数据点来拟合时间序列模型？

- 取决于要估计的模型参数的数量和数据的随机性。

- 对于 ARIMA 建模通常给出最小样本量为30。

- 对于数据量有限的数据，推荐简单模型。

---
# 预测现实 4：历史数据太长

- 大多数时间序列模型对长期时间序列都不适用。
    1. 真实数据并非来自我们使用的模型。当观测值数量不太大（比如大约200个）时，模型通常可以很好的近似生成数据的过程。但最终当你得到足够多的数据时，真正的数据生成过程与模型之间的差异开始变得更加明显。
    2. 由于观测样本量增加，模型参数的优化变得更加耗时。

- 如何处理这些问题取决于模型的目的。可以使用更加灵活和复杂的模型，但这仍然需要假定模型结构适用于整个数据周期。更好的方法通常是允许模型本身随时间变化。

- 如果你只对预测接下来的少数观测值感兴趣，一个简单的方法就是抛弃最早的观测值，只用最近的观测值拟合模型。

---
# 预测现实 5：缺失和异常处理

- 当缺失值引发报错时，至少有两种方法可以解决问题。首先，假如有足够长的观测序列可以用来进行有意义的预测时，我们可以只选取最后一个缺失值之后的数据部分。或者，我们可以用估计值替换缺失值。`na.interp()`函数可以实现这个过程。

- 函数 `tsoutliers()` 函数用于识别离群值，并提供可能的替换值。

- 另一个有用的函数是 `tsclean()`，它能够标识并替换离群值，也能够替换缺失值。

---
# 示例：每天早上的黄金价格

.scroll-output[

```r
gold2 <- na.interp(gold)
autoplot(gold2, series="Interpolated") +
  autolayer(gold, series="Original") +
  scale_color_manual(values=c(`Interpolated`="red",`Original`="gray")) +
  xlab("时间")
```

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

---
# 示例：每天早上的黄金价格

```r
tsoutliers(gold)
```

```
## $index
## [1] 770
## 
## $replacements
## [1] 494.9
```

---
# 示例：每天早上的黄金价格

.scroll-output[

```r
gold %>%
  tsclean() %>%
  ets() %>%
  forecast(h=50) %>%
  autoplot() + 
   xlab("时间")
```

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

---
# 预测现实 6：数据量大且相关

- 全局模型 （M5 竞赛)

- 预测组合 （M 竞赛）

- No-Free-Lunch（M4 竞赛）

---
# 基于特征的预测组合
 
<img src="figures/forecastdiag.jpg" width="90%" style="display: block; margin: auto;" />

---
# 基于特征的组合
 
<img src="figures/fforma.png" style="display: block; margin: auto;" />

---
# 什么是特征？

- 信息熵
- 自相关性
- 异常性
- 等等
- 详见 **tsfeatures** 包。

---
# 怎样学习权重 `$\omega_i$` ？

- 优化理论
- 机器学习
- 例： `xgboost`
- 详见 [M4metalearning 包](https://github.com/robjhyndman/M4metalearning/blob/master/docs/metalearning_example.md)

---
# 预测现实 6：领域相关

参考预测百科论文: https://yanfei.site/docs/paper/ftp.pdf

---
# 课程报告

任选感兴趣的时间序列数据，根据其特点，选择适当的预测方法完成预测。