时间序列是对某个度量在规律时间间隔上的测量。典型的时间序列分析示例包括金融、股票价格、气象数据、公用事业研究等。
时间序列建模可从以下方面入手：
理解生成观测数据的潜在机制与结构。
开始拟合模型，并进行预测、监控，甚至反馈与前馈控制。
在本教程中，您将获得对平稳与非平稳时间序列模型的概览。我们将展示如何通过计算 ACF 和 PACF 来识别时间序列。这些函数的图形有助于判断序列的平稳性。通过差分，我们可以将非平稳序列转化为平稳序列。了解序列的性质后，便可以从其所遵循的模型中更容易地预测未来值。教程还将包含 R 的 TSA 包中的一组真实数据示例。
平稳过程
设 $X$ 为以时间（通常记作 $t$）为索引的随机变量，观测 $\left\{x_t,\,t\in \textbf{N}\right\}$ 称为时间序列。此处 $\textbf{N}$ 为整数集，作为时间索引集；$N$ 也可以是时间戳。平稳性是时间序列模型中的关键假设，意味着序列具有同质性，即无论时间如何变化，序列的行为方式相似，其统计性质不会随时间改变。平稳性有强与弱两种形式。
定义 1
若平稳过程 $\left\{x_t,\,t\in \textbf{N}\right\}$ 在时间尺度平移后其统计分布保持不变，则称其为严格（或强）平稳。由于随机过程的分布由有限维分布函数刻画，我们可以给出严格平稳的等价定义。若对任意 $n$，任意时刻选择 $t_1, t_2, \ldots, t_n\in \textbf{N}$，以及任意满足 $t_{i+k}\in\textbf{N}$ 的滞后 $k$，$n$ 维随机向量 $(X_{t_1+k}, X_{t_2+k}, \ldots, X_{t_n+k})$ 与向量 $(X_{t_1}, X_{t_2}, \ldots, X_{t_n})$ 同分布，则该过程为严格平稳。即对 $h$ 和 $x_i$ 有
\begin{align*} P(X_{t_1}\leq x_1, X_{t_2}\leq x_2, \ldots, X_{t_k}\leq x_k) &= F(x_{t_1}, x_{t_2}, \ldots, x_{t_k})\\ &= F(x_{h+t_1}, x_{h+t_2}, \ldots, x_{h+t_k})\\ &= P(X_{h+t-1}\leq x_1, X_{h+t_2}\leq x_2, \ldots, X_{h+t_k}\leq x_k) \end{align*}
对任意时间平移 $h$ 与观测 $x_j$ 均成立。若 $\left\{X-t,\,t\in\textbf{N}\right\}$ 严格平稳，则 $X_t$ 的边缘分布与 $t$ 无关。此外，二维分布 $(X_{t_1}, X_{t_2})$ 与 $t_1$ 与 $t_2$ 的绝对位置无关，只与间隔 $t_1-t_2$ 有关。因此，均值函数 $E(X)$ 为常数，协方差 $Cov(X_t,X_{t-k})$ 仅为 $k$ 的函数，而与 $k$ 与 $t$ 的绝对位置无关。在更高阶矩（如三阶矩）下，若对 $s, t, u$ 同加一常数时间平移，$E[X_uX_tX_v]$ 亦保持不变。
定义 2
若单变量时间序列 $X_t$ 的均值、方差与协方差均与时间无关，则称其为平稳。于是，若时间序列（或按时间排序的随机变量）$X_t$ 定义在 $t=1,2,3,\ldots, n$ 以及 $t=0, -1, -2, -3, \ldots$ 上，则当且仅当满足
\begin{align} (i) & E[X_t] = \mu\\ (ii) & E\left[(X-\mu)^2\right] = Var(X_t) = \gamma(0) = \sigma^2\\ (iii) & E\left[(X_t-\mu)(X_{t-k}-\mu)\right] = Cov(X_t, X_{t-k}) = \gamma(k) \end{align}
时，$X_t$ 在数学上为弱平稳。前两条（(i) 与 (ii)）分别要求过程具有恒定均值与恒定方差，而 (iii) 要求任意两值之间的协方差（形式上称为协方差函数）仅取决于二者之间的时间间隔 $k$，而与具体时点 $t$ 无关。
若过程为高斯且具有有限二阶矩，则弱平稳等价于强平稳。严格平稳仅在必要矩存在时才蕴含弱平稳。强平稳还需要分布方面的假设。强形式通常被认为过于严格，因此我们主要关注弱平稳，也称协方差平稳、广义平稳或二阶平稳。
平稳性与随机时间序列（随机过程）
在一个时间序列中，若观测围绕恒定均值波动、方差恒定且相互独立，则称为随机时间序列。此类序列不呈现任何模式：
观测不呈上升或下降趋势
方差不随时间增大或减小
观测不会在某些时期显著偏大
一个平稳随机模型的例子可写为
$$X_t = \mu+\varepsilon_t$$ 其中 $\mu$ 为常数均值，满足 $E[X_t] = \mu$；$\varepsilon_t$ 为噪声项，假设其均值为零、方差恒定且相互独立（亦称白噪声）。
随机时间序列的模拟
eps <- rnorm( 100 , mean = 0 , sd = 1 )
mu <- 2
X_t <- mu + eps
ts.plot( X_t, main = "Example of (random) stationary time series" , ylab = expression( X[ t] ) ) 技术支持 该模拟过程围绕常数均值 $\mu = 2$ 波动。
自协方差函数
平稳随机过程的理论自协方差函数（ACF）是评估时间序列性质的重要工具。设 $X_t$ 为均值 $\mu$、方差 $\sigma^2$ 的平稳随机过程。滞后 $k$ 处的 ACF，$\gamma(k)$，为 $$\gamma(k) = \frac{\gamma(k)}{\gamma(0)} = \frac{\gamma(k)}{\sigma^2}$$
ACF 是对自协方差的归一化度量，并具有若干性质。
ACF 的性质
$\rho(0) = 1$
ACF 关于滞后是偶函数，即 $\rho(k) = \rho(-k)$
$|\rho(-k)| \leq 1$
上述过程的 ACF 如下图所示。
acf( X_t, main = "Auto-covariance function of X" ) 技术支持 注 1：ACF 的一个特征是非唯一性。即便某个随机过程具有唯一的协方差结构，反之通常不成立：可能存在多个随机过程共享同一 ACF。这导致的模型设定问题在 [@jenkinsd] 中有所说明。
注 2：平稳过程的自相关函数可构造出一种特殊矩阵，称为 Toeplitz 矩阵。它是一种阶数为 $m = t-k$（滞后，包含至 $m-1$ 阶自相关）的方差-协方差矩阵，因而是对角常数、对称且正定的。
\begin{equation} \begin{pmatrix} 1 & \rho(1) & \rho(2) & \cdots & \rho(m-1)\\ \rho(1) & 1 & \rho(1) & \cdots & \rho(m-2)\\ \cdots &\cdots &\cdots &\cdots &\cdots\\ \rho(m-1) & \rho(m-2) & \cdots & \rho(1) & 1 \end{pmatrix} \end{equation}
纯随机过程或白噪声
若离散过程 $\left\{Z_t\right\}$ 中的随机变量 $Z_t$ 构成一组相互独立且同分布（i.i.d.）的序列，则称其为纯随机过程。该定义意味着过程具有恒定均值与方差。
$$\gamma(k) = cov(Z_t, Z_{t+k}) = 0,\quad\forall k\in -3,-2,-1,0,1,2,3, \ldots$$ 由于均值与自协方差函数（acvf）不依赖时间，该过程为二阶平稳。
随机游走过程
若过程 $\left\{X_t\right\}$ 满足 $X_t = X_{t-1}+Z_t$，其中 $Z_t$ 为均值 $\mu$、方差 $\sigma^2_Z$ 的纯随机过程，则称其为随机游走过程。过程通常从 $t = 0$ 开始，且有 $X_1 = Z_0$，即 $X_0 = 0$。于是
\begin{align*} X_1 &= X_0 + Z_1, \quad\text{当 } t = 1\\ X_2 &= X_1 + Z_2 = X_0 + Z_1+ Z_2, \quad\text{当 } t = 2\\ X_3 &= X_2 + Z_3 = X_0 + Z_1+ Z_2 + Z_3, \quad\text{当 } t = 3\\ &\cdots\\ X_t &= X_0 + \sum_{i = 1}^tZ_i \end{align*}
该过程的一阶矩（期望）为 $$E[X_t] = X_0 +\sum_{i = 1}^tE[Z_i] = X_0 + t\mu_z = t\mu_z$$ 方差为 $$Var(X_t) = t\sigma^2_Z$$ 注意均值与方差随时间变化，因此该过程为\textbf{非平稳}。随机游走行为的时间序列示例包括股票价格。
随机游走过程模拟
X <- 0
Z <- rnorm( 100 , mean = 0.5 , sd = 1.5 )
for ( i in 2 : length( Z) ) {
X[ i] <- X[ i- 1 ] + Z[ i]
}
ts.plot( X, main = "Random walk process" ) 技术支持 差分
差分是使时间序列平稳的最常用方法。这是一种特殊的滤波方式，特别适合去除趋势。对季节性数据，一阶差分通常足以使均值达到平稳。设 $X_t = \left\{X_1, X_2,\ldots, X_n\right\}$ 为非平稳时间序列。其对应的平稳序列可表示为
$$\Delta X_{t+1} = X_{t+1}-X_t\quad\text{或 }\,\Delta X_{t} =X_t-X_{t-1}$$ 即所谓一阶差分。若需要二阶差分，可使用算子 $\Delta^2$，即对一阶差分再作差分：
$$\Delta^2X_{t+2} = \Delta X_{t+2} - \Delta X_{t+1}$$
ts.plot( diff( X) ) 技术支持 由此得到的一阶差分围绕常数均值 0 波动。因为从数学上有
$$\Delta X_{t+1} = X_{t+1}-X_t = Z_t$$ 这是一种平稳过程，因为它是均值与方差恒定的纯随机过程。
q 阶移动平均：MA($q$)
设 $\left\{Z_t\right\}$ 为均值为零、方差为 $\sigma^2_Z$ 的纯随机过程。若过程满足
$$X_t = \beta_0Z_t-\beta_1Z_{t-1}-\cdots-\beta_q Z_{t-q}$$ 则称其为 q 阶移动平均过程，其中 $\beta_i,\, i = 1,2,\ldots,q$ 为常数。随机变量 $Z_t, \,t\in\textbf{N}$ 通常缩放使得 $\beta_0=1$。
一阶 MA 的模拟
Z <- rnorm( 100 , mean = 0 , sd = 1.5 )
X <- c( )
for ( i in 2 : length( Z) ) {
X[ i] <- Z[ i] - 0.45 * Z[ i- 1 ]
}
ts.plot( X, main = "Moving Average or order 1 process" ) 技术支持 对 MA(1) 过程，可如下验证三条条件：
\begin{align} E[X_t] &= 0\\ Var(X_t) &= E[X_t^2] - 0 = E\left[Z_t^2-2\beta Z_tZ_{t-1}+\beta^2Z_{t-1}^2\right]= \sigma^2_Z+\beta^2\sigma^2_Z\\ \gamma(k) &= cov(X_t, X_{t+k}) = E[X_tX_{t+k}] - E[X_t]E[X_{t+k}] = E[X_tX_{t+k}]\\ &= E\left[(Z_t-\beta Z_{t-1})(Z_{t+k}-\beta Z_{t-1+k})\right]\\ &= E\left[Z_tZ_{t+k} - \beta Z_tZ_{t-1+k} -\beta Z_{t-1}Z_{t+k} + \beta^2 Z_{t-1}Z_{t-1+k}\right]\\ \end{align}
当 $k = 0$，$\, \gamma(0) = Var(X_t) = \sigma^2_Z(1+\beta^2)$
当 $k = 1$，$\, \gamma(1) = E\left[Z_tZ_{t+1} - \beta Z_t^2 -\beta Z_{t-1}Z_{t+1} + \beta^2 Z_{t-1}Z_{t}\right] = -\beta\sigma^2_Z$
当 $k>1$，$\, \gamma(k) = 0$。
因此，MA(1) 过程在位移超过一个周期时协方差为零。也就是说，它仅具有一个周期的记忆。
MA(1) 的 ACF 为
$$\rho(k) = \frac{\gamma(k)}{\gamma(0)} = \begin{cases}1, &\text{k }=0\\\frac{-\beta}{1+\beta^2}, & \text{k }\pm 1\end{cases}$$
问题： 将上述推导扩展至 MA(2)。
可以看出，样本自相关函数在确定移动平均过程的阶数上很有用：MA(q) 过程的自相关函数 $\rho(k)$ 在前 q 个滞后处为非零（显著偏离零），在 $k > q$ 时为零。MA 的平稳性对 $\left\{\beta_i\right\}$ 无额外限制。但为保证可逆性，$\left\{\beta_i\right\}$ 需要满足一定约束。
p 阶自回归：AR($p$)
设 $\left\{Z_t\right\}$ 为均值为零、方差为 $\sigma_Z^2$ 的纯随机过程。若过程 $\left\{X_t\right\}$ 满足 $$X_t = \alpha_1X_{t-1}+\alpha_2X_{t-2}+\cdots+\alpha_pX_{t-p}+Z_t$$ 则称其为 p 阶自回归过程。在 p 阶自回归过程中，当前观测 $X_t$（例如今日收益）由回溯至第 $p$ 期的过去观测的加权和生成，并叠加一个与过去值独立的随机变量，因此称为自回归。此类过程由 [@greenwood1920inquiry] 引入。上述 AR($p$) 可加入常数项写作
$$X_t = \delta+\alpha_1X_{t-1}+\alpha_2X_{t-2}+\cdots+\alpha_pX_{t-p}+Z_t$$ 其中 $\mu$ 为与时间序列均值相关的常数项，$\alpha_1, \alpha_2,\ldots, \alpha_p$ 可正可负。若 $E[X_t] = E[X_{t-1}] = \cdots = E[X_p] = \mu$，则 AR($p$) 为平稳。于是，
\begin{align*} E[X_t] = \mu &= \delta + \alpha_1\mu+\alpha_2\mu+\cdots+\alpha_p\mu+ 0\\ \Rightarrow \mu &= \frac{\delta}{1-\alpha_1-\alpha_2-\cdots-\alpha_p} \end{align*} 为使上述公式给出常数，需满足 $\alpha_1+\alpha_2+\cdots+\alpha_p<1$。
</1>
一阶 AR 过程：AR(1)
当 $p=1$ 时，$$X_t = \alpha X_{t-1}+Z_t$$ 即为一阶自回归 AR(1)，亦称马尔可夫过程 。利用后移算子 $BX_t = X_{t-1}$，可将 AR(1) 表示为无限阶 MA 过程。由于 $$(1 - \alpha B)X_t = Z_t$$ 因而 \begin{align*} X_t &= \frac{Z_t}{1-\alpha B}\\ &= (1+\alpha B +\alpha^2B^2+ \cdots)Z_t\\ &=Z_t+\alpha Z_{t-1}+\alpha^2Z_{t-2}+\cdots\\ &=Z_t+\beta_1 Z_{t-1}+\beta_2 Z_{t-2}+\cdots\\ \end{align*}
于是 $E[X_t] = 0$，且 $Var(X_t) = \sigma^2_Z(1+\alpha^2+\alpha^4+\cdots)$。当 $|\alpha|<1$ 时，级数收敛。
问题 ：给定 AR(1) 过程 $X_t = \alpha_1X_1+Z_t$ 为均值为零、方差为 $\sigma^2_Z$ 的纯随机过程，且 $\alpha$ 为满足必要条件的常数，推导 $X_t$ 的方差与自协方差函数。
AR(1) 的模拟
alpha = 0.5
Z <- rnorm( 100 , mean = 0 , sd = 1.5 )
X <- rnorm( 1 )
for ( i in 2 : length( Z) ) {
X[ i] <- 0.7 * X[ i- 1 ] + Z[ i]
}
ts.plot( X) 技术支持 问题： 给出 AR(2) 模型关于参数取值的平稳条件。即证明下列条件：
有限阶的 AR 过程（如 p 阶）也可以表示为无限阶 MA 过程。可通过逐次代入或使用后移算子实现。
自回归移动平均过程：ARMA($p,q$)
在建模中，可能需要在同一模型中同时引入 AR 与 MA 项，这便得到混合的自回归–移动平均（ARMA）过程。包含 $p$ 个自回归项与 $q$ 个移动平均项的 ARMA 过程称为 $(p, q)$ 阶，其形式为
$$X_t = \alpha_1X_{t-1}+\alpha_2X_{t-2}+\cdots+\alpha_pX_{t-p}+Z_t-\beta_1Z_{t-1}-\beta_2Z_{t-2}-\cdots-\beta_qZ_{t-q}$$ 使用后移算子 B，可写为
$$\alpha_p(B)X_t = \beta_q(B)Z_t$$ 其中 $\alpha_p(B)$ 与 $\beta_q(B)$ 分别为 $p$ 与 $q$ 阶多项式，满足
\begin{align*} \alpha_p(B) &= \left(1-\alpha_1 B-\cdots-\alpha_p B^p\right)\\ \beta_q(B) &= \left(1-\beta_1 B-\cdots-\beta_p B^p\right) \end{align*}
为保证可逆，$\beta_q(B)$ 的根必须位于单位圆之外。为保证平稳，需 $\alpha_p(B)=0$ 的根位于单位圆之外。并假设 $\alpha_p(B)=0$ 与 $\beta_q(B) = 0$ 无公共根。
ARMA(1,1) 过程
该过程为 $(1,1)$ 阶，其方程为 $$X_t = \alpha_1X_{t-1}+Z_t-\beta_1Z_{t-1}$$
为证明平稳与可逆，需满足 $|\alpha_1|<1$ 与 $|\beta_1|<1$。当 $\alpha < 0$ 时，ARMA(1,1) 退化为 MA(1)；当 $\beta_1<0$ 时，则为 AR(1)。ARMA(1,1) 可使用后移算子转化为纯自回归表示：
$$(1-\alpha_1B)X_t = (1-\beta_1B)Z_t$$
于是
\begin{align*} \frac{1-\alpha_1B}{1-\beta_1B}X_t &= Z_t\\ \pi(B)X_t &= Z_t \end{align*}
其中
\begin{align*} \pi(B) = 1-\pi_1B-\pi_2B^2-\ldots &= \frac{1-\alpha_1B}{1-\beta_1B}\\ (1-\beta_1B)(1-\pi_1B-\pi_2B^2-\ldots) &= 1-\alpha_1B\\ 1-\left[(\pi_1+\beta_1)B - (\pi_2+\beta_2)B^2 - (\pi_3+\beta_3)B^3\right] & = 1-\alpha_1B \end{align*}
比较等式两边 B 的系数，可得未知项 $$\pi_j = \beta_1^{j-1}(\alpha_1-\beta_1),\,\text{当 }\, j\geq1$$
ARMA (1, 1) 的自相关函数
$X_t$ 的期望为 $E[X_t] = \alpha_1E[X_{t-1}]$，在平稳性假设下，$E[X_t] = E[X_{t-1}] = \mu = 0$。这有助于计算自协方差函数，步骤如下：
$X_tX_{t-k} = \alpha_1X_{t-k}X_{t-1}+X_{t-k}Z_t-\beta_1X_{t-k}Z_{t-1}$
两边取期望，得：\begin{align*} E\left[X_tX_{t-k}\right] &= E\left[\alpha_1X_{t-k}X_{t-1}\right]+E\left[X_{t-k}Z_t\right]-\beta_1E\left[X_{t-k}Z_{t-1}\right]\\ \gamma(k)&=\alpha_1\gamma(k-1)+E\left[X_{t-k}Z_t\right]-\beta_1E\left[X_{t-k}Z_{t-1}\right] \end{align*}
当 $k = 0$：
\begin{align*} \gamma(1)&=\alpha_1\gamma(0)+E\left[X_{t-1}Z_t\right]-\beta_1E\left[X_{t-1}Z_{t-1}\right]\\ &=\alpha_1\left(\alpha_1\gamma(1)+ \sigma^2_Z - \beta_1(\alpha_1-\beta_1)\sigma^2_Z\right)-\beta_1\sigma^2_Z\\ &= \alpha^2\gamma(1)+(\alpha_1-\beta_1)(1-\alpha_1\beta_1)\sigma^2_Z\\ \Rightarrow \gamma(1)&=\frac{(\alpha_1-\beta_1)(1-\alpha_1\beta_1)\sigma^2_Z}{1-\alpha^2}\\ \Rightarrow \gamma(0)&=\alpha_1\frac{(\alpha_1-\beta_1)(1-\alpha_1\beta_1)\sigma^2_Z}{1-\alpha^2}+ \sigma^2_Z - \beta_1(\alpha_1-\beta_1)\sigma^2_Z\\ &=\frac{\left(1+\beta_1^2-2\alpha_1\beta_1\right)\sigma_Z^2}{1-\alpha^2} \end{align*}
当 $k \geq 2$：
\begin{align*} \gamma(2)&=\alpha_1\gamma(k-1) \end{align*}
因此，ARMA(1,1) 的 ACF 为：
$$\rho(k) = \begin{cases} 1, &\text{当 }\,k=0\\ \frac{(\alpha_1-\beta_1)(1-\alpha_1\beta_1)}{1+\beta_1^2-2\alpha_1\beta_1},&\text{当 }\,k=1\\ \alpha_1\rho(k-1), &\text{当 }\,k\geq2\\ \end{cases} $$ ARMA(1,1) 的 ACF 同时体现 AR(1) 与 MA(1) 的特征。$\rho(1)$ 同时包含 MA(1) 与 AR(1) 的参数。在 $\rho(1)$ 之后，ARIMA(1,1) 模型的 ACF 与 AR(1) 的 ACF 具有相同的衰减模式。
问题： 求 ARMA(1,1) 过程的偏自相关函数（PACF）。
注意： 时间序列过程的特征常以其 ACF 与 PACF 给出。时间序列分析中最关键的一步，是基于已有数据识别并建立模型，此时 ACF 与 PACF 尚未知。
ARMA(1,1) 过程模拟
Z <- rnorm( 100 , mean = 0 , sd = 1.5 )
X <- rnorm( 1 )
for ( i in 2 : length( Z) ) {
X[ i] <- 0.35 * X[ i- 1 ] + Z[ i] + 0.4 * Z[ i- 1 ]
}
ts.plot( X, main = "ARMA(1,1) process" )
par( mfrow = c( 1 , 2 ) )
acf( X) ; pacf( X) 技术支持 ARIMA($p,d,q$) 过程
自回归差分移动平均模型由以下方程定义：
SARIMA(p,d,q)(P,D,Q) 过程
预测
本节将使用真实时间序列来拟合其最优模型。为此，您将使用 forecast 包。函数 auto.arima 可基于数据拟合并选择最优模型，forecast 函数可实现向前预测 h 个时期。
library( forecast)
library( TSA) 技术支持 示例 1：
data( "co2" )
data( "boardings" )
fit <- auto.arima( co2)
plot( fc <- forecast( fit, h = 15 ) ) 技术支持 示例 2：
data( "boardings" )
fit2 <- auto.arima( boardings[ , "log.price" ] )
plot( fc2 <- forecast( fit2, h = 15 ) ) 技术支持 参考文献
Greenwood, Major, and G Udny Yule. 1920. “An Inquiry into the Nature of Frequency Distributions Representative of Multiple Happenings with Particular Reference to the Occurrence of Multiple Attacks of Disease or of Repeated Accidents.” Journal of the Royal Statistical Society 83 (2). JSTOR: 255–79.
Jenkins, GM. n.d. “D. G. Watts (1968) Spectral Analysis and Its Applications.” San Francisco.
结语
在本教程中，您系统学习了 R 中的时间序列。内容涵盖平稳过程、随机变量模拟、随机时间序列模拟、随机游走过程等；还涉及 p 阶自回归 AR(p)、SARIMA(p,d,q)(P, D, Q) 过程，以及预测。
若您希望进一步学习 R 中的时间序列，请参考以下 DataCamp 课程：
%22%3E%3Cpath%20d%3D%22m-369%200h408v79h-408z%22%20fill%3D%22%237933ff%22%2F%3E%3Cpath%20d%3D%22m-148%200h306v45h-306z%22%20fill%3D%22%23974dff%22%2F%3E%3C%2Fg%3E%3C%2Fsvg%3E)
%22%2F%3E%3C%2FclipPath%3E%3Cg%20clip-path%3D%22url(%23a)%22%3E%3Cpath%20d%3D%22m76%200h484v36h-484z%22%20fill%3D%22%23974dff%22%2F%3E%3Cpath%20d%3D%22m0%200h384v24h-384z%22%20fill%3D%22%2306bdfc%22%20transform%3D%22matrix(1%200%200%20-1%200%2060)%22%2F%3E%3C%2Fg%3E%3C%2Fsvg%3E)
时间序列是对某个度量在规律时间间隔上的测量。典型的时间序列分析示例包括金融、股票价格、气象数据、公用事业研究等。
时间序列建模可从以下方面入手：
在本教程中，您将获得对平稳与非平稳时间序列模型的概览。我们将展示如何通过计算 ACF 和 PACF 来识别时间序列。这些函数的图形有助于判断序列的平稳性。通过差分，我们可以将非平稳序列转化为平稳序列。了解序列的性质后，便可以从其所遵循的模型中更容易地预测未来值。教程还将包含 R 的 TSA 包中的一组真实数据示例。
平稳过程
设 $X$ 为以时间（通常记作 $t$）为索引的随机变量，观测 $\left\{x_t,\,t\in \textbf{N}\right\}$ 称为时间序列。此处 $\textbf{N}$ 为整数集，作为时间索引集；$N$ 也可以是时间戳。平稳性是时间序列模型中的关键假设，意味着序列具有同质性，即无论时间如何变化，序列的行为方式相似，其统计性质不会随时间改变。平稳性有强与弱两种形式。
定义 1
若平稳过程 $\left\{x_t,\,t\in \textbf{N}\right\}$ 在时间尺度平移后其统计分布保持不变，则称其为严格（或强）平稳。由于随机过程的分布由有限维分布函数刻画，我们可以给出严格平稳的等价定义。若对任意 $n$，任意时刻选择 $t_1, t_2, \ldots, t_n\in \textbf{N}$，以及任意满足 $t_{i+k}\in\textbf{N}$ 的滞后 $k$，$n$ 维随机向量 $(X_{t_1+k}, X_{t_2+k}, \ldots, X_{t_n+k})$ 与向量 $(X_{t_1}, X_{t_2}, \ldots, X_{t_n})$ 同分布，则该过程为严格平稳。即对 $h$ 和 $x_i$ 有
\begin{align*} P(X_{t_1}\leq x_1, X_{t_2}\leq x_2, \ldots, X_{t_k}\leq x_k) &= F(x_{t_1}, x_{t_2}, \ldots, x_{t_k})\\ &= F(x_{h+t_1}, x_{h+t_2}, \ldots, x_{h+t_k})\\ &= P(X_{h+t-1}\leq x_1, X_{h+t_2}\leq x_2, \ldots, X_{h+t_k}\leq x_k) \end{align*}
对任意时间平移 $h$ 与观测 $x_j$ 均成立。若 $\left\{X-t,\,t\in\textbf{N}\right\}$ 严格平稳，则 $X_t$ 的边缘分布与 $t$ 无关。此外，二维分布 $(X_{t_1}, X_{t_2})$ 与 $t_1$ 与 $t_2$ 的绝对位置无关，只与间隔 $t_1-t_2$ 有关。因此，均值函数 $E(X)$ 为常数，协方差 $Cov(X_t,X_{t-k})$ 仅为 $k$ 的函数，而与 $k$ 与 $t$ 的绝对位置无关。在更高阶矩（如三阶矩）下，若对 $s, t, u$ 同加一常数时间平移，$E[X_uX_tX_v]$ 亦保持不变。
定义 2
若单变量时间序列 $X_t$ 的均值、方差与协方差均与时间无关，则称其为平稳。于是，若时间序列（或按时间排序的随机变量）$X_t$ 定义在 $t=1,2,3,\ldots, n$ 以及 $t=0, -1, -2, -3, \ldots$ 上，则当且仅当满足
\begin{align} (i) & E[X_t] = \mu\\ (ii) & E\left[(X-\mu)^2\right] = Var(X_t) = \gamma(0) = \sigma^2\\ (iii) & E\left[(X_t-\mu)(X_{t-k}-\mu)\right] = Cov(X_t, X_{t-k}) = \gamma(k) \end{align}
时，$X_t$ 在数学上为弱平稳。前两条（(i) 与 (ii)）分别要求过程具有恒定均值与恒定方差，而 (iii) 要求任意两值之间的协方差（形式上称为协方差函数）仅取决于二者之间的时间间隔 $k$，而与具体时点 $t$ 无关。
若过程为高斯且具有有限二阶矩，则弱平稳等价于强平稳。严格平稳仅在必要矩存在时才蕴含弱平稳。强平稳还需要分布方面的假设。强形式通常被认为过于严格，因此我们主要关注弱平稳，也称协方差平稳、广义平稳或二阶平稳。
平稳性与随机时间序列（随机过程）
在一个时间序列中，若观测围绕恒定均值波动、方差恒定且相互独立，则称为
随机时间序列。此类序列不呈现任何模式：一个平稳随机模型的例子可写为
$$X_t = \mu+\varepsilon_t$$ 其中 $\mu$ 为常数均值，满足 $E[X_t] = \mu$；$\varepsilon_t$ 为噪声项，假设其均值为零、方差恒定且相互独立（亦称白噪声）。
随机时间序列的模拟
该模拟过程围绕常数均值 $\mu = 2$ 波动。
自协方差函数
平稳随机过程的理论自协方差函数（ACF）是评估时间序列性质的重要工具。设 $X_t$ 为均值 $\mu$、方差 $\sigma^2$ 的平稳随机过程。滞后 $k$ 处的 ACF，$\gamma(k)$，为 $$\gamma(k) = \frac{\gamma(k)}{\gamma(0)} = \frac{\gamma(k)}{\sigma^2}$$
ACF 是对自协方差的归一化度量，并具有若干性质。
ACF 的性质
上述过程的 ACF 如下图所示。
注 1：ACF 的一个特征是非唯一性。即便某个随机过程具有唯一的协方差结构，反之通常不成立：可能存在多个随机过程共享同一 ACF。这导致的模型设定问题在 [@jenkinsd] 中有所说明。
注 2：平稳过程的自相关函数可构造出一种特殊矩阵，称为 Toeplitz 矩阵。它是一种阶数为 $m = t-k$（滞后，包含至 $m-1$ 阶自相关）的方差-协方差矩阵，因而是对角常数、对称且正定的。
\begin{equation} \begin{pmatrix} 1 & \rho(1) & \rho(2) & \cdots & \rho(m-1)\\ \rho(1) & 1 & \rho(1) & \cdots & \rho(m-2)\\ \cdots &\cdots &\cdots &\cdots &\cdots\\ \rho(m-1) & \rho(m-2) & \cdots & \rho(1) & 1 \end{pmatrix} \end{equation}
纯随机过程或白噪声
若离散过程 $\left\{Z_t\right\}$ 中的随机变量 $Z_t$ 构成一组相互独立且同分布（i.i.d.）的序列，则称其为纯随机过程。该定义意味着过程具有恒定均值与方差。
$$\gamma(k) = cov(Z_t, Z_{t+k}) = 0,\quad\forall k\in -3,-2,-1,0,1,2,3, \ldots$$ 由于均值与自协方差函数（acvf）不依赖时间，该过程为二阶平稳。
随机游走过程
若过程 $\left\{X_t\right\}$ 满足 $X_t = X_{t-1}+Z_t$，其中 $Z_t$ 为均值 $\mu$、方差 $\sigma^2_Z$ 的纯随机过程，则称其为随机游走过程。过程通常从 $t = 0$ 开始，且有 $X_1 = Z_0$，即 $X_0 = 0$。于是
\begin{align*} X_1 &= X_0 + Z_1, \quad\text{当 } t = 1\\ X_2 &= X_1 + Z_2 = X_0 + Z_1+ Z_2, \quad\text{当 } t = 2\\ X_3 &= X_2 + Z_3 = X_0 + Z_1+ Z_2 + Z_3, \quad\text{当 } t = 3\\ &\cdots\\ X_t &= X_0 + \sum_{i = 1}^tZ_i \end{align*}
该过程的一阶矩（期望）为 $$E[X_t] = X_0 +\sum_{i = 1}^tE[Z_i] = X_0 + t\mu_z = t\mu_z$$ 方差为 $$Var(X_t) = t\sigma^2_Z$$ 注意均值与方差随时间变化，因此该过程为\textbf{非平稳}。随机游走行为的时间序列示例包括股票价格。
随机游走过程模拟
差分
差分是使时间序列平稳的最常用方法。这是一种特殊的滤波方式，特别适合去除趋势。对季节性数据，一阶差分通常足以使均值达到平稳。设 $X_t = \left\{X_1, X_2,\ldots, X_n\right\}$ 为非平稳时间序列。其对应的平稳序列可表示为
$$\Delta X_{t+1} = X_{t+1}-X_t\quad\text{或 }\,\Delta X_{t} =X_t-X_{t-1}$$ 即所谓一阶差分。若需要二阶差分，可使用算子 $\Delta^2$，即对一阶差分再作差分：
$$\Delta^2X_{t+2} = \Delta X_{t+2} - \Delta X_{t+1}$$
由此得到的一阶差分围绕常数均值 0 波动。因为从数学上有
$$\Delta X_{t+1} = X_{t+1}-X_t = Z_t$$ 这是一种平稳过程，因为它是均值与方差恒定的纯随机过程。
q 阶移动平均：MA($q$)
设 $\left\{Z_t\right\}$ 为均值为零、方差为 $\sigma^2_Z$ 的纯随机过程。若过程满足
$$X_t = \beta_0Z_t-\beta_1Z_{t-1}-\cdots-\beta_q Z_{t-q}$$ 则称其为 q 阶移动平均过程，其中 $\beta_i,\, i = 1,2,\ldots,q$ 为常数。随机变量 $Z_t, \,t\in\textbf{N}$ 通常缩放使得 $\beta_0=1$。
一阶 MA 的模拟
对 MA(1) 过程，可如下验证三条条件：
\begin{align} E[X_t] &= 0\\ Var(X_t) &= E[X_t^2] - 0 = E\left[Z_t^2-2\beta Z_tZ_{t-1}+\beta^2Z_{t-1}^2\right]= \sigma^2_Z+\beta^2\sigma^2_Z\\ \gamma(k) &= cov(X_t, X_{t+k}) = E[X_tX_{t+k}] - E[X_t]E[X_{t+k}] = E[X_tX_{t+k}]\\ &= E\left[(Z_t-\beta Z_{t-1})(Z_{t+k}-\beta Z_{t-1+k})\right]\\ &= E\left[Z_tZ_{t+k} - \beta Z_tZ_{t-1+k} -\beta Z_{t-1}Z_{t+k} + \beta^2 Z_{t-1}Z_{t-1+k}\right]\\ \end{align}
当 $k = 0$，$\, \gamma(0) = Var(X_t) = \sigma^2_Z(1+\beta^2)$
当 $k = 1$，$\, \gamma(1) = E\left[Z_tZ_{t+1} - \beta Z_t^2 -\beta Z_{t-1}Z_{t+1} + \beta^2 Z_{t-1}Z_{t}\right] = -\beta\sigma^2_Z$
当 $k>1$，$\, \gamma(k) = 0$。
因此，MA(1) 过程在位移超过一个周期时协方差为零。也就是说，它仅具有一个周期的记忆。
MA(1) 的 ACF 为
$$\rho(k) = \frac{\gamma(k)}{\gamma(0)} = \begin{cases}1, &\text{k }=0\\\frac{-\beta}{1+\beta^2}, & \text{k }\pm 1\end{cases}$$
问题： 将上述推导扩展至 MA(2)。
可以看出，样本自相关函数在确定移动平均过程的阶数上很有用：MA(q) 过程的自相关函数 $\rho(k)$ 在前 q 个滞后处为非零（显著偏离零），在 $k > q$ 时为零。MA 的平稳性对 $\left\{\beta_i\right\}$ 无额外限制。但为保证
可逆性，$\left\{\beta_i\right\}$ 需要满足一定约束。p 阶自回归：AR($p$)
设 $\left\{Z_t\right\}$ 为均值为零、方差为 $\sigma_Z^2$ 的纯随机过程。若过程 $\left\{X_t\right\}$ 满足 $$X_t = \alpha_1X_{t-1}+\alpha_2X_{t-2}+\cdots+\alpha_pX_{t-p}+Z_t$$ 则称其为 p 阶自回归过程。在 p 阶自回归过程中，当前观测 $X_t$（例如今日收益）由回溯至第 $p$ 期的过去观测的加权和生成，并叠加一个与过去值独立的随机变量，因此称为自回归。此类过程由 [@greenwood1920inquiry] 引入。上述 AR($p$) 可加入常数项写作
$$X_t = \delta+\alpha_1X_{t-1}+\alpha_2X_{t-2}+\cdots+\alpha_pX_{t-p}+Z_t$$ 其中 $\mu$ 为与时间序列均值相关的常数项，$\alpha_1, \alpha_2,\ldots, \alpha_p$ 可正可负。若 $E[X_t] = E[X_{t-1}] = \cdots = E[X_p] = \mu$，则 AR($p$) 为平稳。于是，
\begin{align*} E[X_t] = \mu &= \delta + \alpha_1\mu+\alpha_2\mu+\cdots+\alpha_p\mu+ 0\\ \Rightarrow \mu &= \frac{\delta}{1-\alpha_1-\alpha_2-\cdots-\alpha_p} \end{align*} 为使上述公式给出常数，需满足 $\alpha_1+\alpha_2+\cdots+\alpha_p<1$。
</1>
一阶 AR 过程：AR(1)
当 $p=1$ 时，$$X_t = \alpha X_{t-1}+Z_t$$ 即为一阶自回归 AR(1)，亦称马尔可夫过程。利用后移算子 $BX_t = X_{t-1}$，可将 AR(1) 表示为无限阶 MA 过程。由于 $$(1 - \alpha B)X_t = Z_t$$ 因而 \begin{align*} X_t &= \frac{Z_t}{1-\alpha B}\\ &= (1+\alpha B +\alpha^2B^2+ \cdots)Z_t\\ &=Z_t+\alpha Z_{t-1}+\alpha^2Z_{t-2}+\cdots\\ &=Z_t+\beta_1 Z_{t-1}+\beta_2 Z_{t-2}+\cdots\\ \end{align*}
于是 $E[X_t] = 0$，且 $Var(X_t) = \sigma^2_Z(1+\alpha^2+\alpha^4+\cdots)$。当 $|\alpha|<1$ 时，级数收敛。
问题：给定 AR(1) 过程 $X_t = \alpha_1X_1+Z_t$ 为均值为零、方差为 $\sigma^2_Z$ 的纯随机过程，且 $\alpha$ 为满足必要条件的常数，推导 $X_t$ 的方差与自协方差函数。
AR(1) 的模拟
问题： 给出 AR(2) 模型关于参数取值的平稳条件。即证明下列条件：
有限阶的 AR 过程（如 p 阶）也可以表示为无限阶 MA 过程。可通过逐次代入或使用后移算子实现。
自回归移动平均过程：ARMA($p,q$)
在建模中，可能需要在同一模型中同时引入 AR 与 MA 项，这便得到混合的自回归–移动平均（ARMA）过程。包含 $p$ 个自回归项与 $q$ 个移动平均项的 ARMA 过程称为 $(p, q)$ 阶，其形式为
$$X_t = \alpha_1X_{t-1}+\alpha_2X_{t-2}+\cdots+\alpha_pX_{t-p}+Z_t-\beta_1Z_{t-1}-\beta_2Z_{t-2}-\cdots-\beta_qZ_{t-q}$$ 使用后移算子 B，可写为
$$\alpha_p(B)X_t = \beta_q(B)Z_t$$ 其中 $\alpha_p(B)$ 与 $\beta_q(B)$ 分别为 $p$ 与 $q$ 阶多项式，满足
\begin{align*} \alpha_p(B) &= \left(1-\alpha_1 B-\cdots-\alpha_p B^p\right)\\ \beta_q(B) &= \left(1-\beta_1 B-\cdots-\beta_p B^p\right) \end{align*}
为保证可逆，$\beta_q(B)$ 的根必须位于单位圆之外。为保证平稳，需 $\alpha_p(B)=0$ 的根位于单位圆之外。并假设 $\alpha_p(B)=0$ 与 $\beta_q(B) = 0$ 无公共根。
ARMA(1,1) 过程
该过程为 $(1,1)$ 阶，其方程为 $$X_t = \alpha_1X_{t-1}+Z_t-\beta_1Z_{t-1}$$
为证明平稳与可逆，需满足 $|\alpha_1|<1$ 与 $|\beta_1|<1$。当 $\alpha < 0$ 时，ARMA(1,1) 退化为 MA(1)；当 $\beta_1<0$ 时，则为 AR(1)。ARMA(1,1) 可使用后移算子转化为纯自回归表示：
$$(1-\alpha_1B)X_t = (1-\beta_1B)Z_t$$
于是
\begin{align*} \frac{1-\alpha_1B}{1-\beta_1B}X_t &= Z_t\\ \pi(B)X_t &= Z_t \end{align*}
其中
\begin{align*} \pi(B) = 1-\pi_1B-\pi_2B^2-\ldots &= \frac{1-\alpha_1B}{1-\beta_1B}\\ (1-\beta_1B)(1-\pi_1B-\pi_2B^2-\ldots) &= 1-\alpha_1B\\ 1-\left[(\pi_1+\beta_1)B - (\pi_2+\beta_2)B^2 - (\pi_3+\beta_3)B^3\right] & = 1-\alpha_1B \end{align*}
比较等式两边 B 的系数，可得未知项 $$\pi_j = \beta_1^{j-1}(\alpha_1-\beta_1),\,\text{当 }\, j\geq1$$
ARMA (1, 1) 的自相关函数
$X_t$ 的期望为 $E[X_t] = \alpha_1E[X_{t-1}]$，在平稳性假设下，$E[X_t] = E[X_{t-1}] = \mu = 0$。这有助于计算自协方差函数，步骤如下：
当 $k = 0$：
\begin{align*} \gamma(1)&=\alpha_1\gamma(0)+E\left[X_{t-1}Z_t\right]-\beta_1E\left[X_{t-1}Z_{t-1}\right]\\ &=\alpha_1\left(\alpha_1\gamma(1)+ \sigma^2_Z - \beta_1(\alpha_1-\beta_1)\sigma^2_Z\right)-\beta_1\sigma^2_Z\\ &= \alpha^2\gamma(1)+(\alpha_1-\beta_1)(1-\alpha_1\beta_1)\sigma^2_Z\\ \Rightarrow \gamma(1)&=\frac{(\alpha_1-\beta_1)(1-\alpha_1\beta_1)\sigma^2_Z}{1-\alpha^2}\\ \Rightarrow \gamma(0)&=\alpha_1\frac{(\alpha_1-\beta_1)(1-\alpha_1\beta_1)\sigma^2_Z}{1-\alpha^2}+ \sigma^2_Z - \beta_1(\alpha_1-\beta_1)\sigma^2_Z\\ &=\frac{\left(1+\beta_1^2-2\alpha_1\beta_1\right)\sigma_Z^2}{1-\alpha^2} \end{align*}
当 $k \geq 2$：
\begin{align*} \gamma(2)&=\alpha_1\gamma(k-1) \end{align*}
因此，ARMA(1,1) 的 ACF 为：
$$\rho(k) = \begin{cases} 1, &\text{当 }\,k=0\\ \frac{(\alpha_1-\beta_1)(1-\alpha_1\beta_1)}{1+\beta_1^2-2\alpha_1\beta_1},&\text{当 }\,k=1\\ \alpha_1\rho(k-1), &\text{当 }\,k\geq2\\ \end{cases} $$ ARMA(1,1) 的 ACF 同时体现 AR(1) 与 MA(1) 的特征。$\rho(1)$ 同时包含 MA(1) 与 AR(1) 的参数。在 $\rho(1)$ 之后，ARIMA(1,1) 模型的 ACF 与 AR(1) 的 ACF 具有相同的衰减模式。
问题： 求 ARMA(1,1) 过程的偏自相关函数（PACF）。
注意： 时间序列过程的特征常以其 ACF 与 PACF 给出。时间序列分析中最关键的一步，是基于已有数据识别并建立模型，此时 ACF 与 PACF 尚未知。
ARMA(1,1) 过程模拟
ARIMA($p,d,q$) 过程
自回归差分移动平均模型由以下方程定义：
SARIMA(p,d,q)(P,D,Q) 过程
预测
本节将使用真实时间序列来拟合其最优模型。为此，您将使用
forecast包。函数auto.arima可基于数据拟合并选择最优模型，forecast函数可实现向前预测h个时期。示例 1：
示例 2：
参考文献
Greenwood, Major, and G Udny Yule. 1920. “An Inquiry into the Nature of Frequency Distributions Representative of Multiple Happenings with Particular Reference to the Occurrence of Multiple Attacks of Disease or of Repeated Accidents.” Journal of the Royal Statistical Society 83 (2). JSTOR: 255–79.
Jenkins, GM. n.d. “D. G. Watts (1968) Spectral Analysis and Its Applications.” San Francisco.
结语
在本教程中，您系统学习了 R 中的时间序列。内容涵盖平稳过程、随机变量模拟、随机时间序列模拟、随机游走过程等；还涉及 p 阶自回归 AR(p)、SARIMA(p,d,q)(P, D, Q) 过程，以及预测。
若您希望进一步学习 R 中的时间序列，请参考以下 DataCamp 课程：