时间序列分析与预测

第五讲：探索性数据分析

---

黄嘉平

深圳大学中国经济特区研究中心
粤海校区汇文楼1510
https://huangjp.com/

数据的平稳化

为什么需要平稳时间序列？

对于时间序列数据的预测，最重要的是测量序列值之间的依赖关系，例如自相关系数。

• 非平稳序列的自相关系数 \rho(s,t) 依赖时间 s 和 t

• 平稳序列的自相关系数 \rho(h) 仅依赖时间之差 h=s-t

现实中的数据多是非平稳序列，比如例 1.1 中的强生公司的股票收益

非平稳序列之：趋势平稳模型

在众多非平稳序列中，最简单也是最容易处理的一类序列是趋势平稳模型，即序列在其趋势周围是平稳的（或去掉趋势之后是平稳的）。这类模型可以表达为

x_t = \mu_t + y_t

其中 x_t 是目标序列，\mu_t 代表趋势，y_t 是一个平稳序列。

我们可以尝试将 y_t 从 x_t 中分离出来，即消除趋势。例如，如果我们能获得趋势项的估计值 \widehat{\mu}_t ，就可以计算残差

\widehat{y}_t = x_t - \widehat{\mu}_t

然后针对 \widehat{y}_t 进行建模分析。

挪威三文鱼价格：线性趋势

假设线性趋势 \mu_t = \beta_0 + \beta_1 t。此模型的最小二乘估计结果是 \widehat{\mu}_t = -503 + 0.25t，因此消除趋势后的序列是

\widehat{y}_t = x_t + 503 - 0.25 t

挪威三文鱼价格：线性趋势

从 ACF 可以看出，三文鱼价格似乎存在一个3年左右的变化周期。

差分序列

差分是两个相邻时间点上变量的差，例如 x_t - x_{t-1} ，可以写成 \Delta x_t 或 \nabla x_t 。对于很多非平稳的时间序列数据，其差分序列 \{\nabla x_t\} 往往是平稳的。

例如，如果 x_t = \mu_t + y_t 中的趋势项 \mu_t 遵循带漂移项的随机游走模型，即

\mu_t = \delta + \mu_{t-1} + w_t ,

其中 w_t 是白噪声且与 y_t 独立。此时，x_t 的差分序列是

\begin{align*} x_t - x_{t-1} &= (\mu_t - \mu_{t-1}) + (y_t - y_{t-1}) \\ &= \delta + w_t + (y_t - y_{t-1}) \end{align*}

因为 y_t 是平稳的，可以证明 y_t - y_{t-1} 也是平稳的，因此差分序列 \{\nabla x_t\} 是平稳序列。

挪威三文鱼价格：差分

注意：差分法无法估计 y_t，只能产生平稳序列 \nabla x_t。ACF 显示 \nabla x_t 表现出较强的年周期性。

后移算子（backshift operator）

我们也可以针对差分序列 \nabla x_t 再取一次差分，也就是二阶差分。一阶差分可以移除线性趋势，那么二阶差分就可以移除二次趋势（即趋势项是时间的二次函数）。以此类推，我们也可以计算更高阶的差分。为了简化符号，我们定义后移算子。

后移算子（backshift operator）

后移算子 B 定义为

B x_t = x_{t-1}

因此，B^2 x_t = B(Bx_t) = Bx_{t-1} = x_{t-2}。以此类推，可得 B^k x_t = x_{t-k}。B^{-1} 称为前移算子。

1 阶差分 \nabla x_t = x_t - Bx_t = (1-B)x_t。因此，1-B 可以作为差分算子。d 阶差分可以表达为

\nabla^d x_t = (1-B)^d x_t

例如，二阶差分是 \nabla^2 x_t = (1-B)^2 x_t = (1-2B+B^2)x_t = x_t -2x_{t-1} + x_{t-2}。

全球地表温度数据（例 1.2）

全球地表温度数据的差分

含有滞后项的回归

滞后散点图（lagplot）

滞后散点图是序列 x_t 和自身的滞后项 x_{t-h} 或另一序列的滞后项 y_{t-h} 之间的散点图，通过观察两者间的关系（或添加拟合的回归曲线），可以帮助我们判断应把哪一项加入回归模型的解释变量。

• 图中的曲线是局部加权平滑法（locally weighted scatter plot smoothing, LOWESS）的拟合结果
• 右上角的数字是样本自相关系数
• 自相关系数测量的是两个变量间的线性相关关系，如果 LOWESS 拟合结果明显不是直线，则说明样本自相关系数意义不大

南方涛动指数 SOI 和新鱼数量间的滞后散点图

用 SOI 预测新鱼数量：含有滞后项的回归

样本 CCF 显示，SOI 指数的滞后项 S_{t-5}, S_{t-6}, S_{t-7}, S_{t-8} 和新鱼数量 R_t 间存在较强的负相关。通过 LOWESS 拟合则可以看出，在 SOI 取正值和负值时，相关性的强弱是不同的。因此，我们可以假设下面的回归模型：

R_t = \beta_0 + \beta_1 S_{t-6} + \beta_2 D_{t-6} + \beta_3 S_{t-6} \times D_{t-6} + w_t , \quad D_{t} = \begin{cases} 1 & \text{if }\ S_t \geq 0 \\ 0 & \text{if }\ S_t < 0 \end{cases}

这里我们只用了 S_{t-6} 作为独立变量，但是加入了虚拟变量 D_{t-6} 以及交互项，这样就可以在 S_{t-6} 取正值和负值时获得不同的线性拟合结果。这个模型的另一种写法是

R_t = \begin{cases} \beta_0 + \beta_1 S_{t-6} + w_t & \text{if }\ S_{t-6} \geq 0 \\ (\beta_0 + \beta_2) + (\beta_1 + \beta_3) S_{t-6} + w_t & \text{if }\ S_{t-6} < 0 \end{cases}

用 SOI 预测新鱼数量：含有滞后项的回归

library(zoo) # zoo allows easy use of the variable names
dummy = ifelse(soi<0, 0, 1)
fish = as.zoo(ts.intersect(rec, soiL6=lag(soi,-6), dL6=lag(dummy,-6))) 
summary(fit <- lm(rec~ soiL6*dL6, data=fish, na.action=NULL))

Call:
lm(formula = rec ~ soiL6 * dL6, data = fish, na.action = NULL)

Residuals:
    Min      1Q  Median      3Q     Max 
-63.291 -15.821   2.224  15.791  61.788 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   74.479      2.865  25.998  < 2e-16 ***
soiL6        -15.358      7.401  -2.075   0.0386 *  
dL6           -1.139      3.711  -0.307   0.7590    
soiL6:dL6    -51.244      9.523  -5.381  1.2e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.84 on 443 degrees of freedom
Multiple R-squared:  0.4024,    Adjusted R-squared:  0.3984 
F-statistic: 99.43 on 3 and 443 DF,  p-value: < 2.2e-16

用 SOI 预测新鱼数量：含有滞后项的回归

plot(fish$soiL6, fish$rec, panel.first=Grid(), col="dodgerblue3") 
points(fish$soiL6, fitted(fit), pch=3, col=6) 
lines(lowess(fish$soiL6, fish$rec), col=4, lwd=2) 

用 SOI 预测新鱼数量：含有滞后项的回归

注意：残差依然不像是白噪声，后面我们会尝试寻找更合适的模型。

平滑方法

移动平均

以序列 \{x_t\} 定义的对称移动平均序列定义为

m_t = \sum_{j= -k}^k a_j x_{t-j}

其中 a_j = a_{-j} \geq 0, \ \sum_{j=-k}^k a_j = 1。

核平滑（kernel smoothing）

以核函数（kernel function）定义权重的移动平均过程称为核平滑（kernel smoothing），即

m_t = \sum_{i=1}^n w_i(t) x_{i}

其中 w_i(t) = K\big(\frac{t-t_i}{b}\big)\big/\sum_{k=1}^n K\big(\frac{t-t_k}{b}\big) 是权重，K(\cdot) 是核函数。通常使用正态核（normal kermel）K(z)=\exp(-z^2/2)。

LOWESS

LOWESS 的方法比较复杂，但其核心思想接近最近邻回归（nearest-neighbor regression）。

k 最近邻回归是用 \{x_{t-k/2}, \dots, x_{t}, \dots, x_{t+k/2}\} 中的观测值通过回归预测 x_t，并令 m_t = \hat{x}_t。LOWESS 在这个基础上会对观测值进行加权，距离 x_t 越近的观测值权重越大。

R 中有 lowess() 和 loess() 两个函数。基本的 LOWESS 拟合使用 lowess() 即可，而 loess() 会丢弃数据的时间序列属性，但可以指定回归函数并输出置信区间。

经典结构模型（decomposition）

时间序列分析的一个传统方法是将时间序列数据分解成趋势（T_t），季节性（S_t）和噪声（N_t）三部分，即

x_t = T_t + S_t + N_t

当数据存在循环周期时，也可加入周期项（C_t）。

需要注意的是，这种分解没有唯一解。根据不同的方法或不同的参数设定可以获得不同的结果。最简单的分解方法是：

1. 估计趋势项 T_t；
2. 针对去趋势的序列 x_t - \hat{T}_t，用季节虚拟变量进行回归以估计 S_t；
3. 获得残差 \hat{N}_t = x_t - \hat{T}_t - \hat{S}_t。

R 中可用的命令包括 decompose() 和 stl() 。前者利用了移动平均平滑法，后者则使用 loess 估计趋势项。

夏威夷酒店入住率：使用 decompose()

夏威夷酒店入住率：使用 stl()