21、平稳时间序列分析全解析

平稳时间序列分析全解析

1. 一维自相关结构

时间序列是按时间顺序排列的数字序列，在时间序列分析和建模中，它被视为随机过程的实现，即一系列随机变量。若时间序列的均值恒定，且任意两个时刻值的协方差仅与时间间隔有关，则称该时间序列为二阶平稳时间序列，可用公式表示为：
- (E (x_t) = μ)
- (cov (x_t, x_s) = γ (t - s))

1.1 自协方差/相关函数

平稳时间序列 (x_t) 的自协方差函数定义为：
- (γ (τ) = cov (x_{t + τ}, x_t) = E (x_{t + τ} - μ) (x_t - μ))
其方差 (σ^2 = γ (0))，自相关函数为：
- (ρ(τ) = \frac{γ (τ)}{σ^2})

自协方差函数具有以下性质：
- (|γ (τ)| ≤ γ (0) = σ^2)
- (γ (τ) = γ (-τ))
- 对于任意 (p) 个整数 (τ_1, τ_2, \cdots, τ_p) 和实数 (a_1, a_2, \cdots, a_p)，有 (\sum_{i,j = 1}^{p} γ (τ_i - τ_j)a_ia_j ≥ 0)，即自协方差函数是非负定或半正定的。

我们可以通过 (var (λx_t + x_{t + τ}) ≥ 0)（对于任意实数 (λ)）来推导这些性质，对于最后一个性质，可利用 (var(\sum_{i} a_ix_{τ_i}) ≥ 0)。

1.2 时间序列模型

白噪声是一系列独立同分布的随机变量，均值为零，方差为 (σ^2_ε)，其自协方差是狄拉克脉冲，即 (γ_ε(τ) = δ_τ)（(τ = 0) 时为 1，其他时刻为 0）。但白噪声在实际中并不存在，气候等时间序列通常具有自相关性，因此人们提出了自回归移动平均（ARMA）模型来解释这种自相关性。

以下是一些基本的时间序列操作符：
- 后移或延迟算子 (B) ：对于离散时间序列，(Bx_t = x_{t - 1})，一般地，(B^mx_t = x_{t - m})，其逆算子 (B^{-1}) 为前移算子。对于常数 (c)，(Bc = c)，且 (B^mB^n = B^{m + n})。当 (|α| < 1) 时，(\frac{1}{1 - αB} x_t = (1 + αB + α^2B^2 + \cdots)x_t = x_t + αx_{t - 1} + \cdots)。
- 差分算子 (\nabla = 1 - B) ：(\nabla x_t = (1 - B)x_t = x_t - x_{t - 1})，例如 (\nabla^2x_t = (1 - B)^2x_t = x_t - 2x_{t - 1} + x_{t - 2})。
- 季节性差分算子 (\nabla_k = 1 - B^k) ：常用于处理季节性问题。
- 增益算子 ：对时间序列进行简单线性乘法，即 (ax_t)，参数 (a) 称为增益。
- 连续时间序列的差分算子 (D) ：对于连续时间序列 ({y(t), a ≤ t ≤ b})，(Dy(t) = \frac{dy(t)}{dt})（若可微）。
- 连续时间序列的移位算子 (B_u) ：(B_uy(t) = y(t - u))，且 (B_u = e^{-uD} = e^{-u\frac{d}{dt}})，可通过对 (y(t - u)) 进行泰勒展开推导该式。

ARMA 模型包括以下几种：
- 自回归模型 (AR(p)) ：(x_t = φ_1x_{t - 1} + φ_2x_{t - 2} + \cdots + φ_px_{t - p} + ε_t = (\sum_{k = 1}^{p} φ_kB^k)x_t + ε_t)。当 (p = 1) 时，为一阶自回归（AR(1)）模型，也称为红噪声，常用于气候研究。
- 移动平均模型 (MA(q)) ：(x_t = ε_t + φ_1ε_{t - 1} + \cdots + φ_qε_{t - q} = (1 + \sum_{k = 1}^{q} φ_kB^k)ε_t)。
- 自回归移动平均模型 (ARMA(p, q)) ：((1 - \sum_{k = 1}^{p} φ_kB^k)x_t = (1 + \sum_{k = 1}^{q} θ_kB^k)ε_t)，可简写为 (φ(B)x_t = θ(B)ε_t)，其平稳性要求 (φ(z) = 0) 的根在单位圆外。

确定时间序列模型的方法有多种，例如可根据 ARMA 模型的自相关函数（阻尼指数和/或正弦波）来选择模型，也可使用偏自相关函数。还可以利用信息理论，通过拟合一系列模型，计算残差估计值 (\hat{ε}) 及其方差 (\hat{σ}^2)，然后计算赤池信息准则（AIC）：
- (AIC = log(\hat{σ}^2) + \frac{2}{n}(P + 1))
其中 (P) 是待估计的参数数量，AIC 最小的模型为最佳模型。

2. 功率谱

假设平稳时间序列 (x_t) 的自协方差函数可求和，即 (\sum_{k} γ (k) < ∞)，其功率谱定义为：
- (f (ω) = \frac{1}{2π}\sum_{k = -∞}^{∞} γ (k)e^{-ikω})
利用自协方差函数的对称性，功率谱可表示为：
- (f (ω) = \frac{σ^2}{2π}(1 + 2\sum_{k = 1}^{∞} ρ(k)coskω))

功率谱具有以下性质：
- (f (-ω) = f (ω))
- (f (ω) ≥ 0)，对于所有 (ω \in [-π, π])
- (γ (τ) = \int_{-π}^{π} e^{iωτ}f (ω)dω = \int_{-π}^{π} cosτωf (ω)dω)，即自协方差函数是功率谱的逆傅里叶变换，且 (σ^2 = \int_{-π}^{π} f (ω)dω)，功率谱用于分配方差。

以下是一些功率谱的例子：
- 白噪声的功率谱是常数，即 (f (ω) = \frac{σ^2}{2π})
- 对于红噪声时间序列 (x_t = αx_{t - 1} + ε_t)，自相关函数为 (ρ(τ) = α^{|τ|})，功率谱为 (f (ω) = \frac{σ^2}{2π}(1 - 2αcosω + α^2)^{-1})

计算 ARMA 过程功率谱的直接方法是利用线性滤波的结果。对于延迟操作 (y_t = Bx_t)，其傅里叶变换关系为 (y(ω) = e^{iω}x(ω))；当 (y_t = αx_t + βBx_t) 时，(y(ω) = (α + βe^{iω})x(ω))。对于 ARMA 时间序列模型，其功率谱为 (f_x(ω) = \frac{σ^2_ε}{|\frac{θ(e^{iω})}{φ(e^{iω})}|^2})，其中要求 (φ(z)) 的根在单位圆外（平稳性），(θ(z)) 也满足类似条件（可逆性）。

3. 多元情况

多元时间序列 (x_t) 中的每个元素 (x_t = [x_{t1}, x_{t2}, \cdots, x_{tp}]) 是 (p) 维的，假设其均值为零，协方差矩阵为 (Σ_0)。

3.1 自协方差结构

滞后交叉或自协方差矩阵 (Σ(τ)) 定义为 (Σ(τ) = E[x_{t + τ}x^T_t])，其元素 ([Σ(τ)] {ij} = E[x {t + τ, i}x_{t, j}])。对角元素是构成 (x_t) 的各个一维时间序列的自协方差，非对角元素是滞后交叉协方差。滞后协方差矩阵具有以下性质：
- (Σ(-τ) = [Σ(τ)]^T)
- (Σ(0)) 是 (x_t) 的协方差矩阵 (Σ_0)
- (Σ(τ)) 是半正定的，即对于任意整数 (m > 0) 和实向量 (a_1, \cdots, a_m)，有 (\sum_{i,j = 1}^{m} a^T_i Σ(i - j)a_j ≥ 0)

滞后交叉相关矩阵 (Υ(τ) = Σ^{-1/2} 0 Σ(τ)Σ^{-1/2}_0)，其元素 (ρ {ij}(τ) = \frac{γ_{ij}(τ)}{\sqrt{γ_{ii}(0)γ_{jj}(0)}}) 具有类似性质，且 (|ρ_{ij}(τ)| ≤ 1)。

3.2 交叉谱

多元时间序列 (x_t) 的谱密度矩阵 (F(ω)) 定义为自协方差矩阵的傅里叶变换：
- (F(ω) = \frac{1}{2π}\sum_{τ = -∞}^{∞} e^{-iτω}Σ(τ))
当 (\sum_{τ} |Σ(τ)| < ∞) 时，(F(ω)) 存在。与一元情况不同，谱密度矩阵可能是复数，其对角元素是构成 (x_t) 的各个一维时间序列的功率谱，实部是同谱矩阵，虚部是正交谱矩阵。

谱密度矩阵具有以下性质：
- (F(-ω) = [F(ω)]^{ T})（Hermitian 性质）
- (Σ(τ) = \int_{-π}^{π} F(ω)e^{iτω}dω)
- (Σ_0 = \int_{-π}^{π} F(ω)dω)，且 (2πF(0) = \sum_{k} Σ(k))
- (F(ω)) 是半正定的（Hermitian），即对于任意整数 (m > 0) 和复数 (c_1, c_2, \cdots, c_p)，有 (c^{ T} F(ω)c = \sum_{i,j = 1}^{p} c^{*} i F {ij}(ω)c_j ≥ 0)

对于 (i \neq j)，(x_{t, i}) 和 (x_{t, j}) 之间的相干性和相位分别为：
- (c_{ij}(ω) = \frac{|F_{ij}(ω)|^2}{F_{ii}(ω)F_{jj}(ω)})
- (φ_{ij}(ω) = Atan(\frac{Im(F_{ij}(ω))}{Re(F_{ij}(ω))}))
相干性衡量了两个时间序列在频域中的相关系数的平方，相位衡量了它们之间的时间滞后。

4. 样本空间中的自相关结构

4.1 自协方差/自相关估计

对于有限时间序列样本 (x_t)，常用的自协方差函数估计方法有：
- (\hat{γ} 1(τ) = \frac{1}{n}\sum {t = 1}^{n - τ} (x_t - \bar{x}) (x_{t + τ} - \bar{x}))
- (\hat{γ} 2(τ) = \frac{1}{n - τ}\sum {t = 1}^{n - τ} (x_t - \bar{x}) (x_{t + τ} - \bar{x}))

假设样本均值为零，(\hat{γ}_1(τ)) 有轻微偏差，但渐近无偏，且是一致估计（样本量趋于无穷时方差趋于零）；(\hat{γ}_2(τ)) 无偏，但不一致（样本量增大时方差趋于无穷）。对于固定滞后，两者渐近无偏，且近似方差满足 (var[\hat{γ}_1(τ)] \approx O(\frac{1}{n})) 和 (var[\hat{γ}_2(τ)] \approx O(\frac{1}{n - k}))。

自相关函数可通过 (\hat{ρ}(τ) = \frac{\hat{γ} (τ)}{\hat{σ}^2}) 估计，其中 (\hat{γ} ()) 是自协方差函数的估计，(\hat{σ}^2 = \frac{1}{n - 1}\sum_{t = 1}^{n} (x_t - \bar{x})^2) 是样本方差。样本估计 (\hat{ρ}_1(τ)) 是半正定的，而 (\hat{ρ}_2(.)) 一般不满足。样本自相关函数与滞后的关系图称为相关图，可基于渐近正态性和白噪声进行显著性检验，即 (E[\hat{ρ}(τ)] \approx 0)（(τ \neq 0)），(var[\hat{ρ}(τ)] \approx \frac{1}{n})（(τ \neq 0)），利用这些近似值可构建样本自相关函数的置信区间。

4.2 周期图

• 原始周期图 ：对于中心化的时间序列样本 (x_t)，在谱估计中通常考虑傅里叶频率 (\omega_k = \frac{2πk}{n})，其中 (k = -[\frac{n - 1}{2}], \cdots, [\frac{n}{2}])，([x]) 是 (x) 的整数部分。奈奎斯特频率 (\frac{2π}{2\Delta t}) 是可分辨的最高频率，功率谱只能估计低于该频率的值。
  时间序列的周期图定义为傅里叶系数幅值的平方，即 (I_n(\omega_k) = \frac{1}{n}|\sum_{t = 1}^{n} x_te^{-it\omega_k}|^2)。由周期图的定义可得 ((n - 1)\hat{σ}^2 = \sum_{t = 1}^{n} x^2_t = \sum_{k = -[\frac{n - 1}{2}]}^{[\frac{n}{2}]} |α_k|^2 = \sum_{k = -[\frac{n - 1}{2}]}^{[\frac{n}{2}]} I_n(\omega_k))，说明周期图分配了样本方差。通过展开周期图的定义式可得 (I_n(\omega_p) = \sum_{k = -(n - 1)}^{n - 1} \hat{γ} (k) cos(\omega_pk))，因此 (\frac{1}{2π}I_n(\omega_p)) 是功率谱 (f (ω_p)) 的一个候选估计，且 (E[I_n(\omega_p)] \approx 2πf (ω_p))，即周期图是功率谱的渐近无偏估计，但它不是一致估计，因为其方差恒定，且具有较大的采样波动，需要进行平滑处理。
• 周期图平滑 ：构建谱密度函数一致估计的常用方法是平滑处理。平滑后的周期图通过将原始周期图与“谱窗” (W()) 卷积得到：
  (\hat{f} (ω) = \frac{1}{2π}\sum_{k = -[\frac{n - 1}{2}]}^{[\frac{n}{2}]} W(ω - \omega_k)I_n(\omega_k))
  谱窗是对称的核函数，积分值为 1，且在大值处衰减。这种平滑等价于使用（时域）滞后窗 (\lambda(.)) 对加权自协方差估计进行离散傅里叶变换：
  (\hat{f} (ω) = \frac{1}{2π}\sum_{k = -(n - 1)}^{n - 1} \lambda(k) \hat{γ} (k) cos(\omega_pk))
  该求和通常在滞后窗的截断点处截断。谱窗 (W()) 是滞后窗的傅里叶变换，其目的是忽略样本自协方差函数中来自大滞后的贡献，这意味着时间上的局部化与频域上的宽泛性相关，反之亦然。常见的滞后/谱窗有巴特利特窗和帕曾窗：
  • 巴特利特窗 ：滞后窗定义为 (\lambda(τ) = \begin{cases} 1 - \frac{τ}{M}, & |τ| < M \ 0, & \text{otherwise} \end{cases})，对应的谱窗为 (W(ω) = \frac{M}{n}(\frac{sin(πMω)}{πMω})^2)
  • 帕曾窗 ：滞后窗定义为 (\lambda(τ) = \begin{cases} 1 - 6(\frac{τ}{M})^2 + 6(\frac{τ}{M})^3, & |τ| ≤ \frac{M}{2} \ 2(1 - \frac{τ}{M})^3, & \frac{M}{2} < |τ| < M \ 0, & \text{otherwise} \end{cases})，谱窗为 (W(ω) = \frac{6}{πM^3}(\frac{sin(Mω/4)}{sinω/2})^4)

随着参数 (M) 的增大，滞后窗变窄，由于 (M) 可视为时间分辨率，因此方差随 (M) 增大而增大，反之亦然。除了平滑周期图法，还可以使用最大熵方法（MEM）来估计功率谱，即通过拟合自回归模型到时间序列，然后利用模型参数计算功率谱。交叉协方差和交叉谱的估计方法与样本协方差函数和样本谱类似，例如两个零均值时间序列样本 (x_t) 和 (y_t) 的交叉协方差可通过 (\hat{γ} {12}(τ) = \frac{1}{n}\sum {t = 1}^{n - τ} x_ty_{t + τ}) 估计（(τ = 0, 1, \cdots, n - 1)），并通过对称性补充 (\hat{γ} {12}(-τ) = \hat{γ} {21}(τ))，交叉谱可通过 (\hat{f} {12}(ω) = \frac{1}{2π}\sum {k = -M}^{M} \lambda(k) \hat{γ}_{12}(k)e^{iωk}) 估计。

综上所述，平稳时间序列分析在多个领域都有重要应用，掌握这些理论和方法对于处理和分析时间序列数据具有重要意义。通过对自协方差、自相关、功率谱等概念的理解和运用，我们可以更好地建模和预测时间序列的行为。同时，不同的估计方法和操作符为我们提供了丰富的工具，以适应不同的实际需求。

平稳时间序列分析全解析

5. 应用与案例分析

平稳时间序列分析在众多领域都有着广泛的应用，下面通过几个具体案例来展示其实际应用价值。

5.1 气候数据预测

在气候研究中，许多气候变量的时间序列具有一定的自相关性。以气温时间序列为例，我们可以使用 ARMA 模型进行建模和预测。

步骤如下 ：
1. 数据收集 ：收集某地区多年的每日气温数据，形成时间序列 (x_t)。
2. 数据预处理 ：检查数据的平稳性，如果不平稳，可进行差分处理使其平稳。
3. 模型选择 ：计算样本自相关函数和偏自相关函数，根据其特征选择合适的 ARMA(p, q) 模型阶数 (p) 和 (q)。也可以使用 AIC 准则进行模型选择，选择 AIC 最小的模型。
4. 模型拟合 ：使用选定的模型对数据进行拟合，估计模型参数 (φ_1, φ_2, \cdots, φ_p) 和 (θ_1, θ_2, \cdots, θ_q)。
5. 预测 ：使用拟合好的模型对未来的气温进行预测。

例如，若选择 AR(1) 模型 (x_t = φ_1x_{t - 1} + ε_t)，通过数据拟合得到 (φ_1 = 0.8)，则可以根据当前的气温 (x_t) 预测下一天的气温 (x_{t + 1} = 0.8x_t + ε_{t + 1})。

5.2 金融市场分析

在金融市场中，股票价格、汇率等时间序列数据通常也具有一定的自相关性。以股票价格为例，我们可以使用功率谱分析来研究其波动特征。

步骤如下 ：
1. 数据收集 ：收集某股票的历史价格数据，形成时间序列 (x_t)。
2. 计算功率谱 ：根据功率谱的定义 (f (ω) = \frac{1}{2π}\sum_{k = -∞}^{∞} γ (k)e^{-ikω})，计算股票价格时间序列的功率谱。可以使用样本自协方差函数 (\hat{γ} (k)) 来估计功率谱。
3. 分析功率谱特征 ：观察功率谱的形状和峰值，分析股票价格波动的主要频率成分。例如，如果功率谱在某个频率处有明显的峰值，说明该股票价格在该频率下有较强的波动。
4. 风险评估 ：根据功率谱分析结果，评估股票价格的波动风险。例如，功率谱的高频成分较多，说明股票价格波动较为剧烈，风险较高。

6. 总结与展望

平稳时间序列分析是一种强大的数据分析工具，它通过对时间序列的自相关结构、功率谱等进行分析，能够揭示时间序列的内在规律和特征。本文介绍了一维和多元时间序列的自相关结构、时间序列模型、功率谱等基本概念和方法，以及样本空间中的自相关结构估计和周期图分析。

在实际应用中，平稳时间序列分析可以用于预测、风险评估、信号处理等多个领域。然而，实际的时间序列数据往往具有复杂性和不确定性，可能存在非平稳性、非线性等问题。因此，未来的研究方向可以包括：
1. 非平稳时间序列分析 ：研究如何处理非平稳时间序列，如趋势分析、季节性调整等。
2. 非线性时间序列模型 ：开发更复杂的非线性时间序列模型，以更好地描述实际数据的特征。
3. 多变量时间序列分析 ：进一步研究多变量时间序列之间的关系和相互作用，如格兰杰因果检验、协整分析等。
4. 大数据和机器学习应用 ：结合大数据和机器学习技术，提高时间序列分析的效率和准确性。

总之，平稳时间序列分析在不断发展和完善，未来将在更多领域发挥重要作用。通过不断探索和创新，我们可以更好地理解和利用时间序列数据，为决策提供更有力的支持。

7. 常见问题解答

在平稳时间序列分析过程中，可能会遇到一些常见问题，下面对这些问题进行解答。

问题	解答
如何判断时间序列是否平稳？	可以通过观察时间序列的均值和方差是否随时间变化来初步判断。也可以使用单位根检验等方法进行严格检验，如 ADF 检验、PP 检验等。
如何选择 ARMA 模型的阶数 (p) 和 (q)？	可以根据样本自相关函数和偏自相关函数的特征来初步确定阶数。也可以使用 AIC、BIC 等信息准则进行模型选择，选择信息准则值最小的模型。
周期图为什么需要平滑处理？	周期图是功率谱的渐近无偏估计，但它具有较大的采样波动，方差恒定，不是一致估计。平滑处理可以减少采样波动，提高估计的稳定性和准确性。
不同的滞后/谱窗有什么区别？	不同的滞后/谱窗具有不同的形状和性质，对功率谱估计的效果也不同。例如，巴特利特窗具有较好的平滑效果，但分辨率较低；帕曾窗在分辨率和平滑效果之间有较好的平衡。

8. 流程图总结

下面是一个平稳时间序列分析的基本流程图：

graph TD;
    A[数据收集] --> B[数据预处理];
    B --> C[平稳性检验];
    C -- 平稳 --> D[模型选择];
    C -- 非平稳 --> E[差分处理];
    E --> C;
    D --> F[模型拟合];
    F --> G[模型评估];
    G -- 合格 --> H[预测/分析];
    G -- 不合格 --> D;

这个流程图展示了平稳时间序列分析的基本步骤，从数据收集开始，经过数据预处理、平稳性检验、模型选择、拟合和评估，最终进行预测或分析。如果模型评估不合格，则需要重新选择模型。通过这个流程图，可以清晰地了解平稳时间序列分析的整体流程。