时间序列分析与预测
第二讲：时间序列模型简介
黄嘉平
深圳大学中国经济特区研究中心
粤海校区汇文楼1510
https://huangjp.com/
什么是时间序列模型？
为了利用已知的时间序列数据预测未来的变化趋势，我们需要了解数据是如何生成的。然而世界是复杂的，而可获得的数据是有限的，因此我们无法完全理解数据的生成机制。
科学的发展告诉我们，通过简单的模型描述复杂的现象在很多时候是有效的。
描述时间序列数据 \{x_t\} 的模型是一个能够产生 \{x_t\} 的随机变量序列 \{X_t\} 1。
白噪声 white noise
最简单的时间序列模型是白噪声模型。它是一组相互不相关的随机变量 w_t, 且每个随机变量都具有零均值和有限方差 \sigma_w^2 。
随机过程 \{w_t\}，其中 E(w_t) =0，\mathit{Var}\,(w_t) = \sigma_w^2 < \infty。 也可以写成 w_t \sim \mathit{wn}\,(0, \sigma_w^2)。我们常常假设独立正态分布，即 w_t \sim \mathrm{iid}N(0, \sigma_w^2)。
library(astsa)
w <- rnorm(500) # 生成 500 个服从 N(0,1) 的随机数
tsplot(w, col=4, main="white noise")
移动平均、平滑、滤波
白噪声过程看起来充满了震荡，这意味着它在局部的变化幅度非常大。很多时候我们希望看到更加平滑（smooth）的时间序列数据。最常用的平滑方法是取相邻观测值的平均值，称为移动平均（moving average）。
令 v_t = \tfrac{1}{3} (w_{t-1} + w_t + w_{t+1}) ，则 \{v_t\} 是 \{w_t\} 的一个移动平均过程。
从时间序列数据中剔除噪声的方法称为滤波。比较常见的线性滤波指时间序列数据的线性结合。因此，移动平均过程是线性滤波的一种。
白噪声过程与移动平均滤波
w <- rnorm(500)
v <- filter(w, sides = 2, filter = rep(1/3, 3)) # 用 filter 命令生成 3 期移动平均滤波
s <- filter(w, sides = 2, filter = rep(1/15, 15)) # 15 期移动平均滤波
tsplot(w, col = "tomato3", lwd = 1)
lines(v, col = "blue4", lwd = 3)
lines(s, col = "chartreuse2", lwd = 6)
![]()
自回归模型 autoregression models
考虑下面的方程：
x_t = 1.5 x_{t-1} - 0.75 x_{t-2} + w_t, \quad w_t \sim N(0,1)
此时，时间序列 \{x_t\} 是一个自回归序列。
自回归序列 autoregressive series
给定初始值 x_0 = x_{-1} = 0 时，\{x_t\} 可以用递归形式表达：
\begin{align*}
x_1 &= 1.5 x_{0} - 0.75 x_{-1} + w_1 = w_1 \\
x_2 &= 1.5 x_{1} - 0.75 x_{0} + w_2 = 1.5w_1 + w_2 \\
x_3 &= 1.5 x_{2} - 0.75 x_{1} + w_3 = 1.5w_1 + 1.5w_2 + w_3 \\
&\ \ \vdots
\end{align*}
自回归模型 autoregression models
set.seed(90210) # 指定随机数生成器的种子
w <- rnorm(250 + 50) # 多生成 50 个噪声项
x <- filter(w, filter=c(1.5,-0.75), method="recursive")[-(1:50)]
# method="recursive" 对应自回归模型 [-(1:50)] 代表去掉前 50 个值
tsplot(x, main="autoregression", col=4)
![]()
带漂移项的随机游走模型
白噪声的累积称为随机游走模型，即
x_t = \sum_{j=1}^t w_j \quad \Leftrightarrow \quad x_t = x_{t-1} + w_t
带漂移项的随机游走模型 random walk with drift
在随机游走模型中加入固定趋势即为带漂移项的随机游走模型，即
x_t = \delta t + \sum_{j=1}^t w_j \quad \Leftrightarrow \quad x_t = \delta + x_{t-1} + w_t
\delta 称为漂移项 drift。
带漂移项的随机游走模型
set.seed(314159265)
w <- rnorm(200); x <- cumsum(w) # 在命令后加分号 ; 代表一个命令结束
wd <- w +.3; xd <- cumsum(wd) # cumsum() 为取第一项至当前项之和
tsplot(xd, ylim=c(-2,80), main="random walk", ylab="", col=4)
abline(a=0, b=.3, lty=2, col=4) # 画出漂移项（趋势）
lines(x, col="darkred"); abline(h=0, col="darkred", lty=2)
![]()
信号 + 噪声
在信号处理领域，含有信息的量称为信号（signal），而不含信息的量称为噪声。现实中，任意一个物理观测量都可以看作信号，例如气温、气压、浓度、电压、电流、速度等。非物理观测量也可以是信号，例如股价等。
时间序列数据可以看作在不同时间点上观测的信号或噪声。
周期性信号常用三角函数表达，例如下面模型中的第一项：
x_t = 2\cos \Big( 2\pi \frac{t+15}{50}\Big) + w_t
余弦波可以表达为 A\cos (2\pi \omega t + \varphi) ，其中 A 为振幅（amplitude），\omega 为频率，\varphi 为初始相位。
信号 + 噪声
t <- 1:500
cs <- 2*cos(2*pi*(t+15)/50) # 信号 signal
w <- rnorm(500) # 噪声 noise
par(mfrow=c(3,1))
tsplot(cs, col=4, main=expression(2*cos(2*pi*(t+15)/50)))
tsplot(cs+w, col=4, main=expression(2*cos(2*pi*(t+15)/50)+N(0,1))) # 白噪声方差为 1
tsplot(cs+5*w,col=4, main=expression(2*cos(2*pi*(t+15)/50)+N(0,5^2))) # 白噪声方差为 25
![]()
思考题
随机游走模型
x_t = x_{t-1} + w_t, \quad w_t \sim N(0,1)
R 命令为
set.seed(135)
w <- rnorm(500)
x <- cumsum(w)
tsplot(x)
基于模型生成数据
时间序列分析的目的是通过样本数据判断背后的生成机制，即从数据推断模型。
那么，想象一下你在非洲大草原上或加拿大的针叶林中看到一种动物，你怎么判断它是什么物种呢？
你需要事先掌握物种的知识，例如去动物园参观，或利用百科全书学习。
同理，我们需要学习不同的模型，掌握不同模型生成数据的特点。而学习模型，除了数学推导，最好的方法就是用计算机模拟生成数据。
用计算机生成伪随机数
人们把无法预测或无法解释的现象称为“随机”。因此，所有用计算机生成的数值都不是随机的。但是，我们可以生成很难预测的近似随机的数值，称为伪随机数（pseudo-random number）。
线性同余法 linear congruential generator
最古老的伪随机数生成器。根据下面的递归方程生成伪随机数列 \{x_n\}
x_{n+1} = (a x_n + c) \mod m
其中，0 < a < m ， 0 < c < m ，0 < x_0 < m。
- 例如，(a, c, m, x_0) = (3, 4, 15, 1) 时，\{x_n\} = (1, 7, 10, 4, 1, 7, 10, 4, 1, \dots)，周期为 4。
- 当 (a, c, m) 满足特定条件时，\{x_n\} 的周期为 m 且不受 x_0 取值的影响。
- Borland C/C++ 中 (a, c, m) = (22695477, 1, 2^{31}) 。
- Microsoft Visual C++ 中 (a, c, m) = (214013, 2531011, 2^{31}) 。
R 中的伪随机数生成器
x1 = runif(1000); x2 = rnorm(1000); x3 = rexp(1000, rate=1)
par(mfrow = c(1,3))
hist(x1, main="Uniform"); hist(x2, main="Standard Normal"); hist(x3, main="Exponential")
Monte Carlo Simulation
带有随机成分的模拟方法统称为蒙特卡洛法，该名称源自摩纳哥的著名赌场所在地 Monte-Carlo。
当模型中含有随机变量时，我们感兴趣的量就不再是某一个标量，而是一个概率分布。例如
Y = X + 3Z, \quad X \sim N(1,4^2), \ Z \sim N(5, 1) , \ X \perp Z
我们想知道 Y 的分布。如果你有足够的统计学知识，你会知道 Y \sim N(16, 5^2) 。
蒙特卡洛法则是基于 X 和 Z 的信息，利用伪随机数生成多个 Y 的样本点，从而模拟 Y 的分布。当模型非常复杂时，推导目标变量的分布变得困难甚至不可能，此时蒙特卡洛法就能发挥重要作用。
Monte Carlo Simulation
x <- rnorm(10000, mean=1, sd=4)
z <- rnorm(10000, mean=5, sd=1)
y <- x + 3*z
par(mar = c(5, 4, 1, 2) + 0.1) # 调整图周围的空白空间
plot(density(x), col=2, lwd=4, xlim=c(-15, 35), ylim=c(0,0.42), xlab="", main="")
lines(density(z), col=3, lwd=4); lines(density(y), col=4, lwd=4)
legend("topright", inset=c(0.02,0.05), legend=c("X", "Z", "Y"), col=2:4, lwd=4)
![]()
模拟自回归过程
模拟自回归过程 \quad x_t = 1.5 x_{t-1} - 0.75 x_{t-2} + w_t, \ w_t \sim N(0,1)
w <- rnorm(100) # 准备白噪声序列
x <- vector("numeric", 100) # 在内存中为 x_t 留出空间
x[1] <- w[1] # 定义初期值 x_1
x[2] <- 1.5*x[1] + w[2] # 定义初期值 x_2
for (i in 3:100) { # 将序号 i 从 3 至 100 进行迭代
x[i] <- 1.5*x[i-1] - 0.75*x[i-2] + w[i] # 定义 x_i
}
这是任何编程语言都可以使用的方法，也和数学模型的描述最接近。
w <- rnorm(100)
x <- filter(w, filter=c(1.5,-0.75), method="recursive")
这个方法利用了 stats 包里的命令 filter()。
利用 R 中的独特函数可以大幅减少编程时间，也不需要考虑算法优化、内存等问题，但前提是对 R 语言和各种软件包非常熟悉。
思考题
模拟序列 \{y_t: t = 1, 2, \dots, 200\} ：
\begin{align*}
y_t &= 8 + 0.7 y_{t-1} + 1.4 x_{t-1} - 0.5 x_{t-2} + 0.1t + w_t, \quad w_t \sim N(0, 1) \\
x_t &= x_{t-1} + v_t , \quad v_t \sim N(0, 2^2), \quad v_t \perp w_t \\
\end{align*}
w <- rnorm(200); v <- rnorm(200, sd=____); x <- cumsum(v)
y <- vector("numeric", 200)
y[1] <- 8 + 0.1 + w[1]
y[2] <- 8 + 0.7*y[1] + 1.4*x[1] + 0.2 + w[2]
for (i in ____) {
y[i] <- 8 + 0.7*y[i-1] + 1.4*x[i-1] - 0.5*x[i-2] + ____ + w[i]
}
w <- rnorm(200); v <- rnorm(200, sd=2); x <- cumsum(v)
y <- vector("numeric", 200)
y[1] <- 8 + 0.1 + w[1]
y[2] <- 8 + 0.7*y[1] + 1.4*x[1] + 0.2 + w[2]
for (i in 3:200) {
y[i] <- 8 + 0.7*y[i-1] + 1.4*x[i-1] - 0.5*x[i-2] + 0.1*i + w[i]
}
