商业预测

layout: true

---
class: inverse, center, middle

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

<br>

# 商业预测
## 第四讲：ARIMA 模型及应用

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

---

# 案例回顾：耗电量预测

---

.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="BF-L4-arima_files/figure-html/eledaily-lm1-out-1.png" 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="BF-L4-arima_files/figure-html/eledaily-lm2-out-1.png" 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="BF-L4-arima_files/figure-html/eledaily-lm3-out-1.png" style="display: block; margin: auto;" />
]

---
# 模型残差

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

---
# 自相关性

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

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

---
# AR (Autoregression)

- 自回归 `$$\mathrm{Demand}_t = a_0 + a_1 \times  \mathrm{Demand}_{t-1} + \mathrm{Error}_t.$$`

- 时间序列中的趋势通常在我们估计AR模型之前被过滤掉，方法是首先对数据进行一阶差分。

---
# Integration

- 差分 `$$\mathrm{Demand}_{t} - \mathrm{Demand}_{t-1} = a_{0} + a_{1} \times \left(\mathrm{Demand}_{t-1} - \mathrm{Demand}_{t-2}\right) + \mathrm{Error}_{t}.$$`
    
- 也就是 `$$\Delta \mathrm{Demand}_{t} = a_{0} + a_{1} \times \Delta \mathrm{Demand}_{t-1} + \mathrm{Error}_{t}.$$`

、

---
# 平稳性和差分

- 平稳的时间序列的性质不随观测时间的变化而变化。

- 因此具有趋势或季节性的时间序列是平稳时间序列吗？

- 它们是平稳序列吗？

---
# 差分

- 我们注意到（a）中谷歌股价数并不平稳，但（b）中谷歌股价每天的变化量则是平稳的。这向我们展示了一种让非平稳时间序列变平稳的方法——计算相邻观测值之间的差值，这种方法被称为**差分**。

- 差分则可以通过去除时间序列中的一些变化特征来平稳化它的均值，并因此消除（或减小）时间序列的趋势和季节性。

- 有时差分后的数据仍然不平稳，所以可能需要再一次对数据进行差分来得到一个平稳的序列。

---
# 季节性差分

时间序列中的季节通常在我们估计AR模型之前被过滤掉，方法是首先对数据进行季节差分。季节性差分是对一个观测值和相对应的前一年的观测值之间进行差分。因此有：
`$$y_t' = y_t - y_{t-m}，$$`
其中 `$m=$` 一年中的季节数量。

---
# 季节性差分

.pull-left[

```r
cbind("销售量 ($百万)" = a10,
      "每月销量对数" = log(a10),
      "每年销量变化对数" = diff(log(a10),12)) %>%
  autoplot(facets=TRUE) +
    xlab("年份") + ylab("") +
    ggtitle("抗糖尿病药剂销量")+
  theme(plot.title = element_text(hjust = 0.5))
```
]

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

---
# MA (Moving average)

- 未来需求不取决于前一需求，而是取决于前一误差，我们可以将MA(1)模型表示为以下形式 `$$\mathrm{Demand}_{t} = a_{0} + a_{1} \times \mathrm{Error}_{t-1} + \mathrm{Error}_{t}.$$`

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

---
# 自回归模型

- 在多元线性回归模型中，我们通过对多个预测变量的线性组合预测了目标变量。
- 在自回归模型中，我们则是基于目标变量历史数据的组合对目标变量进行预测。
- 一个 `$p$` 阶的自回归模型可以表示如下：
`$$y_{t} = c + \phi_{1}y_{t-1} + \phi_{2}y_{t-2} + \dots + \phi_{p}y_{t-p} + \varepsilon_{t},$$`
这里的 `$\varepsilon_t$` 是白噪声。这就相当于将预测变量替换为*目标变量的历史值*的多元回归。我们将这个模型称为 **AR( `$p$` ) 模型**—— `$p$` 阶自回归模型。

---
# 自回归模型

- 自回归模型在处理拥有复杂特征的时间序列上十分灵活。
- 在自回归模型中，系数 `$\phi_1,\dots,\phi_p$` 的变化将使得时间序列拥有不同的特征。
- 而误差项 `$\varepsilon_t$`的方差则只会改变序列的数值范围，而不会改变它的特征。

---
class: inverse, center, middle
# 移动平均模型

---
# 移动平均模型

- 不同于使用预测变量的历史值来进行回归，移动平均模型（moving average model）使用历史预测误差来建立一个类似回归的模型。 
`$$y_{t} = c + \varepsilon_t + \theta_{1}\varepsilon_{t-1} + \theta_{2}\varepsilon_{t-2} + \dots + \theta_{q}\varepsilon_{t-q},$$`
- 上式中的 `$\varepsilon_t$` 是白噪声。
- 我们将这个模型称之为 MA(q) 模型。

---
# 移动平均模型示例

- `$y_t$` 的每一个值都可以被认为是一个历史预测误差的加权移动平均值.
- 改变 `$\theta_1,\dots,\theta_q$` 这些系数将会使数据显示出不同的时间序列特征。
- 和自回归模型一样，误差项 `$\varepsilon_t$` 的方差之后改变序列的数值范围，而不会改变它的特征。

---
class: inverse, center, middle
# 非季节性ARIMA模型

---
# 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) 模型**。

---
# ARIMA 中的参数定义

|   |                 |
|--:|:----------------|
|  p|自回归模型阶数   |
|  d|差分阶数         |
|  q|移动平均模型阶数 |

---
# 例：美国女性被谋杀率

```r
autoplot(wmurders) +
  xlab("年份") + ylab("谋杀率（每十万）")+
  theme(plot.title = element_text(hjust = 0.5))
```

---
# 自动选择模型

```r
fit <- auto.arima(wmurders)
```

这是一个ARIMA(1,2,1)模型：
`$$y_{t}''=-0.2434 y_{t-1}''-0.8261 \varepsilon_{t-1}+\varepsilon_{t}.$$`

---
# 预测

```r
fit %>% forecast(h=10) %>% autoplot()
```

---
class: inverse, center, middle
# 季节性ARIMA模型

---
# 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()
```

---
# 实验4

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