analyze-time-series-data.qmd

# 时序数据分析 {#sec-analyze-time-series-data}

预测是非常古老的话题，几乎人人都想拥有预测未来的能力，唐朝袁天罡和李淳风的故事至今还广为流传。事实上，古时候只有至高无上的皇帝才可以去问钦天监了解星辰大海和国运命脉。时间序列数据的分析，以及根据分析得到的一般规律进行预测是经久不衰的命题。预测既包含一般规律指向的确定性，又有无法预知的不确定性，且**同时**包含认知局限带来的不确定性，后者往往更大。无休止地渴求往往伴随着巨大的挑战，而更大的挑战则是预测效果常常不能满足期待。

```{r}
#| message: false
#| eval: false

library(quantmod)    # 获取数据
library(ggplot2)     # 可视化
library(ggfortify)   # 静态展示
library(lmtest)      # 格兰杰因果检验
library(dygraphs)    # 交互展示
```

本章主要从以下几个方面展开：数据获取、数据探索、平稳性诊断、时间序列分解、模型拟合和预测。

## 数据获取

Joshua M. Ulrich 开发维护的 [**quantmod**](https://github.com/joshuaulrich/quantmod) 包可以下载国内外股票市场的数据。本节主要以美团股价数据为例，美团自 2018-09-20 在香港挂牌上市，股票代码 3690.HK。首先用 **quantmod** 包 [@quantmod2022] 获取美团上市至 2023-11-24 每天的股价数据，包含 Open 开盘价、High 最高价、Low 最低价、Close 闭市价、Adjusted 调整价和 Volume 成交量数据。

```{r}
#| label: downlaod-yahoo-data
#| eval: false
#| echo: true

library(quantmod)
# 美团股票代码 3690
meituan <- getSymbols("3690.HK", auto.assign = FALSE, src = "yahoo")
```

```{r}
#| echo: false
#| label: load-meituan-data
#| message: false

library(quantmod)
# 美团上市至 2023-11-24
meituan <- readRDS(file = "data/meituan.rds")
```

先来看数据的类型，数据类型颇为复杂，是由 `xts` 和 `zoo` 两种类型复合而成，`xts` 类型是继承自 `zoo` 类型的。

```{r}
class(meituan)
str(meituan)
```

数据集 `meituan` 是一个 `xts` 类型的时间序列数据对象，时间范围是 2018-09-20 至 2023-11-24，包含 4 个成分，分别如下

-   Data 部分显示为 906 行 6 列的双精度浮点存储的数值。
-   Columns 部分显示列名，依次是 3690.HK.Open、3690.HK.High、 3690.HK.Low 和 3690.HK.Close 等，当列数很多时，显示时会省略。
-   Index 部分表示索引列，有序是时间序列数据的本质特点。示例中索引存储数据点产生的先后顺序，索引是用日期来表示的，日期所在的时区是 "UTC"。
-   xts 部分是数据类型的一些属性（元数据），说明数据集的来源，什么时候制作的数据。示例中数据是从雅虎财经下载的，下载时间是 2023-11-27 14:31:12。

与时间序列数据相关的数据类型有很多，比如 Base R 提供的 Date 和 POSIX 等，扩展包 **timeDate** 和 **chron** 也都有自己的一套数据类型及处理方法。[**xts**](https://github.com/joshuaulrich/xts) 包是处理时间序列数据的主要工具之一，xts 是 eXtensible Time Series 的缩写。为了进一步了解用法，下面举个例子，使用该 R 包的函数 `xts()` 构造时间序列对象。

``` r
xts(x = NULL,
    order.by = index(x),
    frequency = NULL,
    unique = TRUE,
    tzone = Sys.getenv("TZ"),
    ...)
```

-   参数 `x` 表示数据。
-   参数 `order.by` 表示索引数据。
-   参数 `frequency` 表示频率。
-   参数 `unique` 表示唯一。
-   参数 `tzone` 表示时区。

```{r}
#| message: false

library(zoo)
library(xts)
# 数据矩阵
x <- matrix(1:4, ncol = 2, nrow = 2)
# 日期索引
idx <- as.Date(c("2018-01-01", "2019-12-12"))
# xts = matrix + index
xts(x, order.by = idx)
```

## 数据探索

### zoo

**zoo** 包提供 S3 范型函数 `autoplot.zoo()` 专门可视化 `zoo` 类型的数据，它接受一个 `zoo` 类型的数据对象，返回一个 `ggplot2` 数据对象，然后用户可以添加自定义的绘图设置，更多详情见帮助文档 `?autoplot.zoo()` 。

```{r}
#| label: fig-meituan-ggplot2
#| fig-cap: "美团在香港上市以来的股价走势"
#| fig-showtext: true

# xts 包需要先加载，否则 Index 不是日期类型而是数值类型
library(ggplot2)
autoplot(meituan[, "3690.HK.Adjusted"]) +
  theme_classic() +
  labs(x = "日期", y = "股价")
```

**zoo** 包还提供另一个范型函数 `fortify()` 将 `zoo` 数据对象转化为 `data.frame` ，这可以方便使用 **ggplot2** 包来展示数据。参数 `melt = TRUE` 意味着重塑原数据集，将数据从宽格式转长格式。参数 `names = c(Index = "Date")` 表示将 Index 列重命名为 date 列。

```{r}
meituan_df <- fortify(
  meituan[, c("3690.HK.Adjusted", "3690.HK.High")],
  melt = TRUE, names = c(Index = "Date")
)
```

数据集 `meituan_df` 中的 Series 列是因子型的，将其标签 `3690.HK.Adjusted` 、`3690.HK.High` 调整为调整价、最高价。根据日期字段 `Date` 提取年份字段 `year` 和一年中的第几天的字段 `day_of_year`。

```{r}
meituan_df <- within(meituan_df, {
  # 调整 Series 的标签
  Series <- factor(Series, labels = c("调整价", "最高价"))
  # 日期字段 Date 获取年份
  year <- format(Date, "%Y")
  # 日期字段 Date 一年中的第几天
  day_of_year <- as.integer(format(Date, "%j"))
})
```

调用 **ggplot2** 包绘制分面、分组时间序列图，以 `day_of_year` 为横轴，股价 `Value` 为纵轴，按 `year` 分组，按 `Series` 分面。

```{r}
#| label: fig-meituan-by-year
#| fig-cap: "美团调整的股价逐年走势"
#| fig-showtext: true

ggplot(data = meituan_df, aes(x = day_of_year, y = Value)) +
  geom_line(aes(color = year)) +
  facet_wrap(~Series, ncol = 1) + 
  theme_classic() +
  labs(x = "一年中的第几天", y = "调整的股价", color = "年份")
```

2019 年底开始出现疫情，2020 年整年陆续有疫情，美团股价一路狂飙突进，因疫情，利好外卖业务，市场看好外卖业务。2021 年政府去杠杆，互联网监管趋严，又监又管，受外部大环境，逆全球化趋势影响，整年股价一路走低。进入 2022 年，股价在 200 附近徘徊。

### xts

```{r}
library(xts)
```

**xts** 包提供 S3 泛型函数 `plot.xts()` 专门用来可视化 `xts` 类型的时间序列数据

```{r}
#| label: fig-meituan-plot
#| fig-cap: "美团在香港上市以来的股价走势"
#| fig-showtext: true

plot(meituan[, "3690.HK.Adjusted"], main = "调整的股价")
```

还可以任意选择一个时间窗口，展示相关数据

```{r}
#| label: fig-meituan-plot-xts
#| fig-cap: "美团 2021 年的股价走势"
#| fig-showtext: true

plot(meituan[, "3690.HK.Adjusted"],
  subset = "2022-01-01/2022-12-31", main = "调整的股价"
)
```

元旦节三天不开市，所以假期没有数据。

### ggfortify

**ggfortify** [@Tang2016] 支持快速地可视化 `ts`、`timeSeries` 、`stl` 等多种类型的时序数据， **ggplot2** 做数据探索会有一些帮助。

```{r}
#| label: fig-meituan-ggfortify
#| fig-cap: 美团股价走势
#| fig-showtext: true

library(ggfortify)
autoplot(meituan[, "3690.HK.Adjusted"], ts.geom = "line") +
  scale_x_date(
    date_breaks = "1 year",
    date_minor_breaks = "6 months",
    date_labels = "%b\n%Y"
  ) +
  theme_classic()
```

### dygraphs

[**dygraphs**](https://github.com/rstudio/dygraphs) 包专门绘制交互式时间序列图形，它封装了时序数据可视化库 [dygraphs](https://github.com/danvk/dygraphs) ，更多情况见 <https://dygraphs.com/>。下面以美团股价为例，展示时间窗口筛选、坐标轴名称、刻度标签、注释、事件标注、缩放等功能。

```{r}
#| label: fig-meituan-dygraphs
#| fig-cap: 美团股价变化趋势
#| fig-width: 6
#| fig-height: 4
#| echo: true
#| eval: !expr knitr::is_html_output(excludes = 'epub')

library(dygraphs)
# 缩放
dyUnzoom <- function(dygraph) {
  dyPlugin(
    dygraph = dygraph,
    name = "Unzoom",
    path = system.file("plugins/unzoom.js", package = "dygraphs")
  )
}

# 年月
getYearMonth <- '
  function(d) {
    var monthNames = ["01", "02", "03", "04", "05", "06","07", "08", "09", "10", "11", "12"];
    date = new Date(d);
    return date.getFullYear() + "-" + monthNames[date.getMonth()]; 
  }'

# 绘图
dygraph(meituan[, "3690.HK.Adjusted"], main = "美团股价走势") |>
  dyRangeSelector(dateWindow = c("2023-01-01", "2023-11-24")) |>
  dyAxis(name = "x", axisLabelFormatter = getYearMonth) |>
  dyAxis("y", valueRange = c(0, 500), label = "美团股价") |>
  dyEvent("2020-01-23", "武汉封城", labelLoc = "bottom") |>
  dyShading(from = "2020-01-23", to = "2020-04-08", color = "#FFE6E6") |>
  dyAnnotation("2020-01-23", text = "武汉封城", tooltip = "武汉封城", width = 60) |>
  dyAnnotation("2020-04-08", text = "武汉解封", tooltip = "武汉解封", width = 60) |>
  dyHighlight(highlightSeriesOpts = list(strokeWidth = 2)) |>
  dySeries(label = "调整股价") |>
  dyLegend(show = "follow", hideOnMouseOut = FALSE) |>
  dyOptions(fillGraph = TRUE, drawGrid = FALSE, gridLineColor = "lightblue") |>
  dyUnzoom()
```

```{r}
#| label: fig-meituan-dygraphs-img
#| echo: false
#| fig-cap: 美团股价变化趋势
#| eval: !expr knitr::is_latex_output()

knitr::include_graphics(path = "screenshots/meituan-stocks.png")
```

上图默认展示 YTD 数据，在一个动态的时间窗口内显示数据，假如今天是 2023-07-15，则展示 2023-01-01 至 2023-07-15 的股价数据。在函数 `dyRangeSelector()` 中设定时间窗口参数 `dateWindow`，实现数据范围的筛选。

## 平稳性诊断

### 自相关图

```{r}
#| label: fig-airpassengers-acf
#| fig-cap: 乘客数量自相关图
#| fig-showtext: true
#| fig-width: 5
#| fig-height: 4

autoplot(acf(AirPassengers, plot = FALSE)) +
  theme_classic()
```

### 偏自相关图

```{r}
#| label: fig-airpassengers-pacf
#| fig-cap: 乘客数量偏自相关图
#| fig-showtext: true
#| fig-width: 5
#| fig-height: 4

autoplot(pacf(AirPassengers, plot = FALSE)) +
  theme_classic()
```

### 延迟算子 {#sec-lag-operator}

```{r}
# 原始序列
AirPassengers
# 延迟 1 期
lag(AirPassengers, k = 1)
```

### 差分算子 {#sec-diff-operator}

函数 `diff()` 实现差分算子，默认参数 `lag = 1` ，`differences = 1` 表示延迟期数为 1 的一阶差分。

```{r}
# 延迟 1 期 1 阶差分
diff(AirPassengers, lag = 1, differences = 1)
```

### 单位根检验 {#sec-unitroot-test}

### 格兰杰因果检验 {#sec-granger-causality-test}

1969 年 Clive Granger 提出格兰杰因果检验，R 语言中 **lmtest** 包的函数 `grangertest()` 可以检验序列中变量之间的时间落差的相关性。

## 指数平滑模型 {#sec-exponential-smoothing}

### 指数平滑

首先来回答何为指数平滑？用历史数据的线性组合预测下一个时期的值，线性组合的权重随距离变远而按指数衰减。不妨设观测序列数据为 $\{x_i\}$ ，预测序列数据为 $\{y_i\}$，用数学公式表达，如下：

$$
y_h(1) = wx_h + w^2x_{h-1} + \cdots = \sum_{j=1}^{\infty} w^j x_{h+1-j}
$$

其中，权重 $0 < w < 1$ ，权重越小表示距离远的历史数据对当前预测的贡献越小。线性组合的权重之和等于 1，所以

$$
\sum_{j=1}^{\infty} w^j = \frac{w}{1-w}
$$

则第 $j$ 个权重应为

$$
\frac{w^j}{\frac{w}{1-w}} = (1-w)w^{j-1},j=1,2,\ldots
$$

则根据历史的 $h$ 期数据预测未来的 1 期数据 $y_h(1)$ 如下：

$$
y_h(1) = (1-w)(x_h + wx_{h-1} + w^2x_{h-2} + \cdots) = (1-w)\sum_{j=0}^{\infty}w^j x_{h-j}
$$

以上就是指数平滑（exponential smoothing），在早期应用中，权重 $w$ 的选取主要依靠经验。适用于没有明显趋势性、季节性、周期性的时间序列数据。

### 函数 `filter()`

函数 `filter()` 实现一元时间序列的线性过滤，或者对多元时间序列的单个序列分别做线性变换，它只是根据既定的平滑模型变换数据，没有拟合数据。函数 `filter()` 实现递归过滤和卷积过滤两种数据变换方式，分别对应自回归和移动平均两种时间序列平滑模型。

-   递归过滤（自回归）

$$
y_{i} = x_{i} + f_1 y_{i-1} +\cdots+ f_p y_{i-p}
$$ {#eq-filter-recursive}

-   卷积过滤（移动平均）

$$
y_{i} = f_1 x_{i+o} + \cdots + f_p x_{i+o-(p-1)}
$$ {#eq-filter-convolution}

其中，$p$ 代表模型的阶数， $o$ 代表漂移项，O 表示英文单词 offset 的首字母。下面举个具体的例子来说明函数 `filter()` 的作用，设输入序列 $\{x_i\}$ 是从 1 至 10 的整数。首先考虑自回归的情况，代码如下：

```{r}
x <- 1:10
# 自回归
filter(x, filter = c(2 / 3, 1 / 6, 1 / 6), method = "recursive")
```

参数 `x` 指定输入的时间序列 $\{x_i\}$，参数 `method` 指定平滑的方法，`method = "recursive"` 表示使用自回归方法，参数 `filter` 表示自回归的系数，系数向量的长度代表模型 @eq-filter-recursive 中的 $p$ ，`filter = c(2 / 3, 1 / 6, 1 / 6)` 对应的模型如下：

$$
\begin{aligned}
y_1 &= x_1 \\
y_2 &= x_2 + \frac{2}{3} y_1 \\
y_3 &= x_3 + \frac{2}{3} y_2 + \frac{1}{6} y_1 \\
y_i &= x_i + \frac{2}{3} y_{i-1} + \frac{1}{6} y_{i - 2} + \frac{1}{6} y_{i - 3}, \quad i \geq 4 \\
\end{aligned}
$$

其中，序列 $\{y_i\}$ 表示函数 `filter()` 的输出结果，由上述方程不难看出自回归模型的递归的特点。为了理解自回归和递归的过程，下面依次计算 $y_1$ 至 $y_4$ 。

```{r}
# y1
1
# y2
2 + 2/3 * 1
# y3
3 + 2/3 * (2 + 2/3 * 1) + 1/6 * 1
# y4
4 + 2/3 * (3 + 2/3 * (2 + 2/3 * 1) + 1/6 * 1) + 1/6 *(2 + 2/3 * 1) + 1/6 * 1
```

接下来，考虑移动平均的情况，代码如下：

```{r}
# 移动平均
filter(x, filter = c(2 / 3, 1 / 6, 1 / 6), method = "convolution", sides = 1)
```

参数 `method = "convolution"` 表示使用移动平均。参数 `sides` 仅适用于卷积过滤，`sides = 1` 表示系数都是作用于过去的值。为了对比自回归和移动平均，不妨设移动平均的系数同自回归的系数，则移动平均模型如下：

$$
\begin{aligned}
y_1 &~~ \text{不存在}\\
y_2 &~~ \text{不存在}\\
y_3 &= \frac{2}{3} x_{3} + \frac{1}{6} x_2 + \frac{1}{6} x_1\\
y_i &= \frac{2}{3} x_{i} + \frac{1}{6} x_{i - 1} + \frac{1}{6} x_{i - 2}, \quad i \geq 3
\end{aligned}
$$

比照模型 @eq-filter-convolution ，漂移项参数 $o$ 为 0，也就是没有漂移，移动平均作用于过去的 3 期数据，也就是 $p = 3$ 。因输出序列 $\{y_i\}$ 中 $y_1,y_2$ 不存在，下面仅计算 $y_3,y_4$ 。

```{r}
# y3
2/3 * 3 + 1/6 * 2 + 1/6 * 1
# y4
2/3 * 4 + 1/6 * 3 + 1/6 * 2
```

**TTR** 包提供许多移动平均的计算函数，比如 `SMA()` ，下面计算过去 3 个观察值的算术平均。

```{r}
library(TTR)
SMA(x, n = 3)
```

### 简单指数平滑

当时间序列不含趋势和季节性成分的时候，可以用简单指数平滑模型来拟合和预测。简单指数平滑模型如下：

$$
\begin{aligned}
\hat{y}_{t+h} &= a_{t} + h \times b_{t} + s_{t - p + 1 + (h - 1) \mod p} \\
a_{t} &= \alpha (y_{t} - s_{t-p}) + (1-\alpha) (a_{t-1} + b_{t-1}) \\
b_{t} &= b_{t-1} \\
s_{t} &= s_{t-p}
\end{aligned}
$$

其中，周期 $p$

```{r}
air_passengers_exp <- HoltWinters(AirPassengers, gamma = FALSE, beta = FALSE)
air_passengers_exp
```

预测的残差平方和 SSE sum-of-squared-errors

```{r}
air_passengers_exp$SSE
```

```{r}
#| label: fig-airpassengers-exp-fitted
#| fig-cap: 简单指数平滑模型
#| fig-showtext: true
#| fig-width: 5
#| fig-height: 4

# plot(air_passengers_exp)
autoplot(air_passengers_exp) +
  theme_classic()
```

向前预测 5 期

```{r}
air_passengers_pred <- predict(air_passengers_exp, n.ahead = 10, prediction.interval = TRUE)
```

预测值及其预测区间

```{r}
#| label: fig-airpassengers-exp-pred
#| fig-cap: 简单指数平滑模型预测
#| fig-showtext: true
#| fig-width: 5
#| fig-height: 5

plot(air_passengers_exp, air_passengers_pred)
```

### Holt 指数平滑 {#sec-holt}

当时间序列不含季节性成分，可以用 Holt 指数平滑模型拟合和预测 [@Holt2004] 。

$$
\begin{aligned}
\hat{y}_{t+h} &= a_{t} + h \times b_{t} + s_{t - p + 1 + (h - 1) \mod p} \\
a_{t} &= \alpha (y_{t} - s_{t-p}) + (1-\alpha) (a_{t-1} + b_{t-1}) \\
b_{t} &= \beta (a_{t} - a_{t-1}) + (1-\beta) b_{t-1} \\
s_{t} &= s_{t-p}
\end{aligned}
$$

```{r}
air_passengers_holt <- HoltWinters(AirPassengers, gamma = FALSE)
air_passengers_holt
```

可知，$\alpha = 1,\beta = 0.0032$

```{r}
#| label: fig-airpassengers-holt-fitted
#| fig-cap: holt 指数平滑模型
#| fig-showtext: true
#| fig-width: 5
#| fig-height: 5

plot(air_passengers_holt)
```

### Holt-Winters 指数平滑 {#sec-holt-winters}

时间序列同时含有趋势成分、季节性成分、随机成分，可以用 Holt-Winters 平滑模型来拟合和预测。根据趋势和季节性的关系，Holt-Winters 平滑模型分为可加 Holt-Winters 平滑和可乘 Holt-Winters 平滑。R 提供函数 `HoltWinters()` 拟合 Holt-Winters 平滑模型[@Holt2004; @Winters1960]。

可加 Holt-Winters 平滑模型如下：

$$
\begin{aligned}
\hat{y}_{t+h} &= a_{t} + h \times b_{t} + s_{t - p + 1 + (h - 1) \mod p} \\
a_{t} &= \alpha (y_{t} - s_{t-p}) + (1-\alpha) (a_{t-1} + b_{t-1}) \\
b_{t} &= \beta (a_{t} - a_{t-1}) + (1-\beta) b_{t-1} \\
s_{t} &= \gamma (y_{t} - a_{t}) + (1-\gamma) s_{t-p}
\end{aligned}
$$

可乘 Holt-Winters 平滑模型如下：

$$
\begin{aligned}
\hat{y}_{t+h} &= (a_{t} + h \times b_{t}) \times s_{t - p + 1 + (h - 1) \mod p} \\
a_{t} &= \alpha (y_{t} / s_{t-p}) + (1-\alpha) (a_{t-1} + b_{t-1}) \\
b_{t} &= \beta (a_{t} - a_{t-1}) + (1-\beta) b_{t-1} \\
s_{t} &= \gamma (y_{t} / a_{t}) + (1-\gamma) s_{t-p}
\end{aligned}
$$

其中 $\alpha, \beta, \gamma$ 是参数，$p$ 为周期长度，$a_{t}, b_{t}, s_{t}$ 分别代表水平、趋势和季节性成分。

```{r}
air_passengers_add <- HoltWinters(AirPassengers, seasonal = "additive")
air_passengers_add
```

可知，$\alpha = 0.248,\beta = 0.0345,\gamma = 1$

```{r}
#| label: fig-airpassengers-add-fitted
#| fig-cap: 可加 Holt-Winters 平滑模型拟合
#| fig-showtext: true
#| fig-width: 5
#| fig-height: 4

autoplot(air_passengers_add) +
  theme_classic()
```

```{r}
air_passengers_mult <- HoltWinters(AirPassengers, seasonal = "mult")
```

```{r}
#| label: fig-airpassengers-mult-fitted
#| fig-cap: 可乘 Holt-Winters 平滑模型拟合
#| fig-showtext: true
#| fig-width: 5
#| fig-height: 4

autoplot(air_passengers_mult) +
  theme_classic()
```

做一个 Shiny 应用展示参数 $\alpha, \beta, \gamma$ 对 Holt-Winters 平滑预测的影响。

## 时间序列分解 {#sec-time-series-decomposition}

-   可加模型

$$
y_t = T_t + S_t + e_t
$$

-   可乘模型

$$
y_t = T_t \times S_t \times e_t
$$

对时间序列 $\{y_t\}$ 分解，趋势性成分 $T_t$、季节性成分 $S_t$、剩余成分 $e_t$

### 函数 `decompose()` {#sec-decompose}

函数 `decompose()` 分解

```{r}
air_decomp_add <- decompose(x = AirPassengers, type = "additive")
```

函数返回一个列表，包含 6 个元素，分别是 `x` 原始序列，`seasonal` 季节性成分，`figure` 估计的季节图，`trend` 趋势成分，`random` 剩余成分，`type` 分解方法。

```{r}
#| label: fig-airpassengers-decomp
#| fig-cap: 变化趋势的分解
#| fig-showtext: true
#| fig-height: 6

# plot(air_decomp_add)
autoplot(air_decomp_add) +
  theme_classic()
```

去掉季节性部分

```{r}
#| label: fig-airpassengers-adjusted
#| fig-cap: 季节性调整
#| fig-showtext: true
#| fig-width: 5
#| fig-height: 5

AirPassengers_adjusted <- AirPassengers - air_decomp_add$seasonal
plot(AirPassengers_adjusted)
```

### 函数 `stl()` {#sec-stl}

函数 `stl()` 将时间序列分解为趋势性成分、季节性成分（周期性）、剩余成分。

```{r}
air_stl <- stl(x = AirPassengers, s.window = 12)
```

```{r}
#| label: fig-airpassengers-stl
#| fig-cap: 变化趋势的分解
#| fig-showtext: true
#| fig-height: 6

autoplot(air_stl) +
  theme_classic()
```

剩余成分不是平稳序列，是异方差的。

**xts** 包的 `periodicity()` 函数可以检测时间序列数据的周期，但时序数据对象最好是在 **xts** 框架内。

```{r}
xts::periodicity(AirPassengers)
```

## 经典时间序列模型 {#sec-classic-time-series-models}

### 自回归模型 {#sec-autoregressive-models}

函数 `ar()` 拟合 AR 模型

```{r}
ar(AirPassengers, order.max = 3)
```

### 移动平均模型 {#sec-moving-average-models}

将自回归的阶设为 0，函数 `arima()` 也可以用来拟合 MA 模型。

```{r}
arima(AirPassengers, order = c(0, 1, 3))
```

### 自回归移动平均模型 {#sec-autoregressive-moving-average-models}

函数 `arima()` 拟合 ARIMA 模型

```{r}
arima(AirPassengers, order = c(1, 1, 3))
```

**forecast** 包提供函数 `auto.arima()` 自动选择合适的自回归、差分和移动平均的阶来拟合数据。

``` r
forecast::auto.arima(AirPassengers)
```

```         
Series: AirPassengers 
ARIMA(2,1,1)(0,1,0)[12] 

Coefficients:
         ar1     ar2      ma1
      0.5960  0.2143  -0.9819
s.e.  0.0888  0.0880   0.0292

sigma^2 = 132.3:  log likelihood = -504.92
AIC=1017.85   AICc=1018.17   BIC=1029.35
```

## 总结 {#sec-time-series-summary}

方法没有好坏，只有适合与否。Holt-Winter 适合预警任务，算法简单，可以及时出预测结果，仅需要一步预测，不需要给出多步预测，要求快，以便迅速作出反应。Prophet 实现的贝叶斯结构可加模型适合短期预测任务，只要在可容许的时间范围内出结果即可，可以迅速出结果当然更好，需要给出多步预测结果，且结果需要强解释性，以便提前做一些商家供给、平台资源的分配。商分模型常常需要比较强的可解释性，算法策略模型重在预测精准度，对可解释性要求不高。

在时间序列数据的可视化方面，除了 Base R 提供的绘图方法外，静态的时序图 **lattice** 和 **ggplot2** 都不错，而交互式图形推荐使用 **plotly** 和 **dygraphs**。

**PortfolioAnalytics** 包做投资组合优化，均值-方差，收益和风险权衡。 [Rmetrics](https://www.rmetrics.org/) 提供系列时间序列数据分析和建模的 R 包，包括投资组合优化 **fPortfolio**、多元分析 **fMultivar、**自回归条件异方差模型 **fGarch**、二元相依结构的 Copulae 分析 **fCopulae** 、市场和基础统计 **fBasics** 。

**fable** 一元到多元时间序列预测问题，提供 ETS、ARIMA、TSLM 等模型，并有书籍时间序列预测原则。值得一提， **forecast** 包开发者 Rob J Hyndman 称已不再开发新的功能，推荐大家使用 **fable** 包。**feasts** 包辅助特征抽取、序列分解、汇总统计和绘制图形等， 插件包 **fable.prophet** 接入 Prophet 的预测能力。[**timetk**](https://github.com/business-science/timetk) 时间序列数据处理、分析、预测和可视化工具箱，提供一致的操作方式，试图形成完成的解决方案。The Rmetrics Association 开发了一系列 R 包专门处理金融时间序列数据，比如 **fGarch** 包提供条件自回归异方差模型。

从时间序列中寻找规律，这样才是真的数据建模，从数据到模型，而不是相反 [Finding Patterns in Time Series](https://mason.gmu.edu/~jgentle/papers/FindingPatternsTimeSeriesDraft.pdf)，识别金融时间序列的模式和统计规律。