阿尔伯塔时间序列分析笔记-全-

阿尔伯塔时间序列分析笔记（全）

001：噪声过程

在本节课中，我们将要学习时间序列分析的基础，特别是不同类型的噪声过程。这些过程是构成所有时间序列数据的基石。我们将从最简单的白噪声开始，逐步介绍自回归过程和移动平均过程，并探讨它们的基本性质。

概述：什么是时间序列分析？

我们的世界在不断变化，作为统计学家，我们的工作是理解这些变化。统计学479课程——时间序列分析，就是为此而生。与线性回归不同，时间序列分析处理的是具有时间依赖性的数据。例如，预测明天的气温，最好的依据可能是今天、昨天甚至更早的气温。这引出了核心挑战：在时间序列中，观测值之间通常不独立，而是存在时间上的依赖关系（或称时间因果性），即过去会影响未来。

在深入分析真实数据之前，我们需要先理解一些数学模型。这引出了本讲的第一部分：噪声的类型。时间序列的核心就是噪声，即随机波动过程。

白噪声：时间序列的基石

我们首先讨论白噪声。这是所有时间序列模型的基础构建模块。

定义：令 W_t 为一个以时间 t 为索引的随机过程，其中 t 在区间 [0, T] 内。白噪声满足以下假设：

零均值：E[W_t] = 0。
恒定方差：Var(W_t) = σ²。
零协方差：对于任意 t ≠ s，Cov(W_t, W_s) = 0。这意味着不同时间点的观测值不相关。

注意：不相关仅意味着没有线性关系，但仍可能存在非线性依赖。

以上定义被称为弱白噪声。我们可以在此基础上加强条件，得到更严格的白噪声：

独立同分布白噪声：在弱白噪声的基础上，将“不相关”替换为“独立”。这为模型提供了更强的基础。
高斯白噪声：在独立同分布白噪声的基础上，增加条件：W_t 服从正态分布，即 W_t ~ N(0, σ²)。正态分布因其数学上的便利性和在实践中的普遍性而备受青睐。

为什么叫“白”噪声？
这个术语源于信号处理。白光包含了所有频率的光。类似地，白噪声在信号意义上包含了所有频率，其图形看起来是完全混乱的波动，这与只包含单一频率的平滑正弦波形成对比。

自回归过程：用过去预测现在

上一节我们介绍了作为基础的白噪声，本节中我们来看看如何用它构建更复杂的过程。第一个是自回归过程。

“回归”一词提示我们，这与线性回归有相似之处。在线性回归中，我们用自变量的线性组合来预测因变量。在时间序列中，我们用变量自身的过去值来预测其当前值。

其一般公式如下：
X_t = θ_1 * X_{t-1} + θ_2 * X_{t-2} + ... + θ_p * X_{t-p} + W_t
其中 θ_i 是固定系数，W_t 是白噪声。

这个公式表明，当前值 X_t 是过去 p 个值的线性组合，再加上当前时刻的噪声。

以下是两个重要的特例：

随机游走：这是自回归过程的一个特例。
X_t = X_{t-1} + W_t
这意味着当前值完全是前一个值加上一个随机扰动。它没有特定的趋势，只是上下随机波动。
带漂移的随机游走：
X_t = A + X_{t-1} + W_t
其中 A 是一个固定常数。如果 A > 0，过程在随机波动的基础上有一个向上的趋势；如果 A < 0，则趋势向下。这可以用来建模具有长期趋势的现象，如全球气温上升或股票价格的整体走势。

通过模拟，我们可以观察这些过程。例如，白噪声图看起来完全随机；随机游走图则显示出相邻时间点值相近，但长期可能变化很大的特点；高阶自回归过程（如AR2、AR3）的图形可能更复杂。仅通过看图，通常难以直接判断时间序列属于哪种自回归过程，这就需要我们后续使用统计方法进行识别。

移动平均过程：平滑噪声

自回归过程构成了本课程的一大基石，另一大基石是移动平均过程。

这是一种平滑过程，在现实世界中很常见。例如，在报告COVID-19病例时，常使用“7日移动平均”来平滑每日波动较大的数据，以更好地观察趋势。

从数学上讲，移动平均过程是对白噪声的平滑。其公式如下：
X_t = φ_1 * W_{t-1} + φ_2 * W_{t-2} + ... + φ_q * W_{t-q} + W_t
这被称为 q 阶移动平均过程。

这个公式表明，当前值 X_t 是过去 q 个噪声项的加权和，再加上当前时刻的噪声。

一个直观的理解是：将 X_t 视为股票价格，W_t 视为市场在时刻 t 受到的随机冲击。那么，移动平均模型意味着过去的市场冲击（“震荡”）会以衰减的方式影响当前价格。如果 φ 为正，一个正向冲击会在未来一段时间内产生持续的正面影响；如果为负，则会产生类似“回弹”的相反效果。

移动平均过程还有另一种表示形式，即对以当前时刻为中心的一个时间窗口内的噪声进行加权平均：
X_t = Σ_{j=-Q/2}^{Q/2} φ_j * W_{t+j}
通常，权重 φ_j 之和为1。当窗口长度（Q值）增大时，得到的图形会越来越平滑。例如，一个MA(21)过程可能看起来与某些自回归过程或随机游走非常相似。这再次说明，仅凭图形难以区分不同类型的时间序列过程。

重要概念与定义

在继续之前，我们需要明确一些贯穿课程的重要概念和定义。

时间索引：

理论上，时间可以是连续的（t ∈ [0, T]）。
在本课程中，我们主要处理离散且等间隔的时间序列数据，即 t = 0, 1, 2, ..., T。这简化了分析，避免了不规则采样带来的复杂性。

以下是两个关键的概率概念：

马尔可夫过程：
一个过程 X_t 是马尔可夫过程，如果给定全部历史信息时，其当前时刻的期望值只依赖于最近一期的观测值。即：
E[X_t | X_1, ..., X_{t-1}] = E[X_t | X_{t-1}]
例如，AR(1)过程就是马尔可夫的，因为 X_t 只直接依赖于 X_{t-1}。

鞅：
鞅可以看作广义的随机游走，有时被称为“公平游戏”。其定义为：
E[X_t | X_1, ..., X_{t-1}] = X_{t-1}
这意味着，基于过去信息，对下一期值的最佳预测就是当前观测值。股票价格等金融数据常被假设具有鞅性质。

高斯过程：
高斯过程是多元正态分布在连续时间上的推广。对于任意有限个时间点 t_1, t_2, ..., t_k，向量 (X_{t_1}, ..., X_{t_k}) 服从多元正态分布。

一个高斯过程完全由其均值函数和协方差函数刻画。
正态分布的性质非常良好：多元正态分布的任意线性组合仍是正态分布。
因此，如果基础白噪声 W_t 是高斯白噪声，那么由它构建出的任何AR或MA过程也都是高斯过程。我们只需要关注这些过程的均值和协方差结构即可。

线性过程：一个统一的框架

最后，我们介绍一个更一般化的概念——线性过程，它涵盖了本课程将讨论的大多数时间序列模型。

定义：如果 W_t 是白噪声（索引 t 从负无穷到正无穷），那么由下式定义的过程 X_t 称为线性过程：
X_t = μ + Σ_{j=-∞}^{∞} θ_j * W_{t-j}
其中 μ 是常数均值项。

这个定义看起来像是移动平均过程的无限阶版本。但首先需要回答一个数学问题：这个无穷和有意义吗？
答案是肯定的，只要系数平方和收敛：Σ_{j=-∞}^{∞} θ_j² < ∞。
这是因为我们需要确保 X_t 的方差是有限的（计算过程涉及方差公式与无穷级数）。

注意：上述定义是非因果的，因为 X_t 依赖于未来（j < 0 时）的噪声。在时间序列中，我们通常关注因果性，即过去影响未来。

因此，我们通常使用因果线性过程的定义：
X_t = μ + Σ_{j=0}^{∞} θ_j * W_{t-j}
这里 j 从0到无穷，W_{t-j} 随着 j 增大而不断回溯到更早的过去，不涉及未来。

示例：将AR(1)过程表示为线性过程

让我们以AR(1)过程为例，看看它如何纳入线性过程的框架。
AR(1)过程定义为：X_t = θ * X_{t-1} + W_t。

我们可以递归地展开这个等式：

X_t = θ * X_{t-1} + W_t
     = θ * (θ * X_{t-2} + W_{t-1}) + W_t
     = θ² * X_{t-2} + θ * W_{t-1} + W_t
     = ...

假设 |θ| < 1，当递归进行到无限过去时，包含初始项 X_{t-k} 的部分会趋于零（因为 θ^k 趋于零）。最终，我们得到：
X_t = Σ_{j=0}^{∞} θ^j * W_{t-j}
这正是因果线性过程的形式，其中 θ_j = θ^j。

现在，我们严谨地考虑这个无穷和。定义部分和：
S_N(θ) = Σ_{j=0}^{N} θ^j * W_{t-j}

其期望为：E[S_N(θ)] = 0。
其方差为：Var(S_N(θ)) = σ² * Σ_{j=0}^{N} θ^{2j} = σ² * (1 - θ^{2(N+1)}) / (1 - θ²)。

当 |θ| < 1 且 N → ∞ 时，方差收敛到：
lim_{N→∞} Var(S_N(θ)) = σ² / (1 - θ²)
此外，如果 W_t 是高斯白噪声，那么 S_N(θ) 的极限分布就是正态分布 N(0, σ²/(1-θ²))。

边界情况：

如果 |θ| > 1，方差级数发散，过程不稳定。
如果 θ = 1，我们就得到了随机游走：X_t = X_{t-1} + W_t。此时，线性过程的表示不收敛，但根据中心极限定理，适当标准化后的部分和会收敛到正态分布。随机游走正好处于平稳与非平稳过程的边界上。

这种对参数 θ 取值的依赖性和边界行为，在时间序列分析中非常重要。

总结

本节课中我们一起学习了时间序列分析的入门知识：

白噪声：介绍了弱白噪声、独立同分布白噪声和高斯白噪声，它们是构建更复杂模型的基础。
自回归过程：学习了如何用变量自身的过去值来建模当前值，包括随机游走和带漂移的随机游走等特例。
移动平均过程：学习了如何通过对过去的噪声项进行加权平均来获得平滑的序列。
重要概念：明确了离散时间索引、马尔可夫性、鞅、高斯过程等关键思想。
线性过程：介绍了一个统一框架，并展示了AR(1)过程如何表示为线性过程，同时讨论了其收敛条件。

这些基本的噪声过程和概念是我们理解、建模和分析真实世界时间序列数据的工具箱。在接下来的课程中，我们将深入探讨时间序列的更多性质，如自协方差、平稳性等。

002：自协方差与平稳性

在本节课中，我们将深入学习时间序列的核心属性：自协方差、自相关和平稳性。理解这些概念是分析任何时间序列数据的基础。

上一讲我们介绍了不同类型的随机过程，如自回归、移动平均和白噪声过程。本节中，我们将进一步探讨时间序列的统计特性，并学习如何在实际数据中估计这些特性。

自协方差与自相关

首先，我们需要回顾协方差的概念。在概率论中，对于两个随机变量 X 和 Y，其协方差定义为：
Cov(X, Y) = E[(X - E[X])(Y - E[Y])]
协方差衡量了两个变量之间的线性关系。一个常见的误解是协方差为零意味着独立，但事实并非如此。协方差为零不意味着独立，但独立一定意味着协方差为零。只有在联合正态分布的情况下，协方差为零才意味着独立。

那么，什么是自协方差呢？对于一个时间序列过程 {X_t}，其在时间点 s 和 t 的自协方差定义为该过程在两个不同时间点的协方差：
K_X(s, t) = Cov(X_s, X_t)
自协方差函数具有两个重要性质：

对称性：K_X(s, t) = K_X(t, s)。
正定性：对于任意有限个时间点集合 {t_1, ..., t_k}，由 K_X(t_i, t_j) 构成的矩阵是正定（或半正定）的。

自相关函数是自协方差的标准化版本，它消除了量纲的影响，其值在 -1 到 1 之间。自相关函数定义为：
ρ_X(s, t) = K_X(s, t) / sqrt(K_X(s, s) * K_X(t, t))
在实际分析中，我们经常使用自相关函数（或其估计值）来判断时间序列中是否存在依赖关系，并初步识别过程的类型。

示例：带漂移项的 AR(1) 过程

为了加深理解，让我们计算一个具体例子：带常数项（漂移项）的 AR(1) 过程。
X_t = a + θ * X_{t-1} + W_t
其中，{W_t} 是方差为 σ² 的弱白噪声。通过递归展开，我们可以将其写成线性过程的形式，进而计算其均值和自协方差。

首先，其均值为常数：
E[X_t] = a / (1 - θ)
其次，其自协方差函数仅依赖于时间差 |s-t|，具体形式为：
K_X(s, t) = (σ² / (1 - θ²)) * θ^{|s-t|}
特别地，其方差（即 s = t 时的自协方差）为：
Var(X_t) = σ² / (1 - θ²)
请注意，这个方差不依赖于时间 t。这是一个非常重要的性质，它将我们引向下一个核心概念。

平稳性

平稳性是时间序列分析中一个至关重要的概念。粗略地说，它意味着时间序列的统计特性不随时间平移而改变。

我们主要关注两种平稳性：

弱平稳性：要求均值恒定，且自协方差仅依赖于时间差，即：
- E[X_t] = μ （对所有 t 成立）
- Cov(X_t, X_{t+τ}) = K_X(τ) （仅依赖于时滞 τ，而非具体时间 t）
严平稳性：要求时间序列任意有限维联合分布都不随时间平移而改变。这比弱平稳性的条件更强。

为什么平稳性如此重要？在统计学中，我们经常需要估计参数（如均值）。如果均值随时间变化，那么用整个时间序列的样本均值去估计就会得到错误的结果。平稳性保证了我们估计的目标是一个恒定值。对于非平稳序列，我们通常需要先通过差分、去趋势等方法将其转化为平稳序列，然后再进行分析。

对于弱平稳过程，其自协方差函数还有两个有用性质：

对称性：K_X(τ) = K_X(-τ)
有界性：|K_X(τ)| ≤ K_X(0) （由柯西-施瓦茨不等式可得）

联合平稳性与样本均值的方差

对于两个时间序列 {X_t} 和 {Y_t}，如果它们各自弱平稳，且其互协方差 Cov(X_s, Y_t) 也仅依赖于时间差 s-t，则称它们是联合平稳的。

现在，让我们看看在时间序列背景下，一个最简单的统计量——样本均值——的估计有何不同。假设我们有一个弱平稳时间序列 {X_t}, t=1,...,T，其真实均值为 μ。我们使用样本均值进行估计：
μ̂ = (1/T) * Σ_{t=1}^{T} X_t
这个估计量是无偏的，即 E[μ̂] = μ。但其方差不再像独立同分布数据那样简单。

在经典统计（IID数据）中，样本均值的方差为：
Var(μ̂_IID) = σ² / T
然而，在时间序列数据存在自相关的情况下，样本均值的方差为：
Var(μ̂) = (1/T) * Σ_{τ=-(T-1)}^{T-1} (1 - |τ|/T) * K_X(τ)
这个公式可以展开为：
Var(μ̂) = (1/T) * K_X(0) + (2/T) * Σ_{τ=1}^{T-1} (1 - τ/T) * K_X(τ)
其中第一项 K_X(0)/T 就是 IID 情况下的方差。第二项是由自相关引起的附加项。如果自协方差 K_X(τ) 多为正，那么样本均值的实际方差将大于 IID 情况下的方差，这意味着我们的估计不如在独立数据中那么精确。理解这一点对于正确评估估计量的不确定性至关重要。

总结

本节课中我们一起学习了时间序列分析的核心基础概念：

自协方差和自相关函数，它们量化了时间序列内部不同时间点之间的线性依赖关系。
平稳性，特别是弱平稳性，它要求均值恒定且自协方差仅依赖于时间差。平稳性是进行许多统计推断的前提。
在平稳 AR(1) 过程示例中，我们看到了理论均值和自协方差的计算方法。
我们探讨了时间序列数据中样本均值估计量的方差会如何受到自相关性的影响，这与独立数据的情形有显著不同。

下一讲，我们将学习如何从实际数据中估计自协方差和自相关函数，并利用这些工具来初步判断一个时间序列是否类似于白噪声过程，这将是正式进行时间序列建模的第一步。

003：估计与线性模型

在本节课中，我们将学习自协方差和自相关函数，如何估计它们，以及它们能告诉我们关于时间序列数据的哪些信息。课程最后，我们将使用R语言分析一个关于20世纪80年代至90年代美国处方药价格的真实数据集。

回顾平稳过程

上一节我们介绍了平稳过程的概念，本节中我们来看看如何估计其参数。对于一个弱平稳过程，其均值不随时间变化，且自协方差函数仅依赖于时间点之间的间隔，而非具体位置。这意味着我们可以使用所有数据点来估计一个全局均值。

对于平稳过程，总体均值的估计量就是通常的样本均值，记作 \(\bar{x}\)：

\[\bar{x} = \frac{1}{T} \sum_{t=1}^{T} X_t \]

然而，由于时间序列数据点之间存在依赖性，该估计量的方差与独立数据的情况有所不同。

估计自协方差与自相关函数

在时间序列分析中，我们真正感兴趣的是估计自协方差和自相关函数，它们是探索性分析的关键工具。

对于一个平稳过程 \(X_t\)，其自协方差函数可以简化为单参数形式 \(K_X(h)\)，其中 \(h\) 是滞后阶数。当 \(h=0\) 时，两个时间点重合；\(h=1, 2, 3...\) 时，它们逐渐远离。

对于滞后阶数 \(h\)（从0到 \(T-1\)），自协方差函数的估计量 \(\hat{K}(h)\) 定义如下：

\[\hat{K}(h) = \frac{1}{T} \sum_{t=1}^{T-h} (X_{t+h} - \bar{x})(X_t - \bar{x}) \]

以下是关于此估计量的几点说明：

求和项的数量为 \(T-h\) 项。这意味着随着 \(h\) 增大，用于估计自协方差的项数减少，估计的准确性会下降。
在实际分析中，我们通常只关注较小的 \(h\) 值（如前10或15阶），因此这个问题通常不构成严重障碍。

接下来，我们可以估计自相关函数 \(\hat{\rho}(h)\)。它只是估计的自协方差除以估计的方差（即 \(h=0\) 时的自协方差）：

\[\hat{\rho}(h) = \frac{\hat{K}(h)}{\hat{K}(0)} \]

这类似于计算相关系数：用协方差除以各自的标准差。由于过程平稳，方差是常数。

此外，我们还可以估计互协方差和互相关函数。对于两个联合平稳的过程 \(X\) 和 \(Y\)，在滞后 \(h\) 处的互协方差估计为：

\[\hat{K}_{XY}(h) = \frac{1}{T} \sum_{t=1}^{T-h} (X_{t+h} - \bar{x})(Y_t - \bar{y}) \]

相应的互相关估计为：

\[\hat{\rho}_{XY}(h) = \frac{\hat{K}_{XY}(h)}{\sqrt{\hat{K}_X(0) \hat{K}_Y(0)}} \]

关于自协方差估计量公式的一个技术细节是：我们在求和时除以 \(T\) 而非 \(T-h\)。这是为了确保估计的协方差矩阵具有正定性，这是协方差矩阵的一个重要数学性质。

在R中分析自相关函数

现在，让我们在R中查看一些估计的自相关函数示例。我们可以使用 acf() 函数来计算和绘制自相关图。该函数默认计算从滞后0到滞后15的自相关估计。

自相关图能直观地揭示时间序列的结构：

在滞后0处，自相关总是1。
对于白噪声过程，我们期望在除0以外的所有滞后处，自相关都接近于0。图中的蓝虚线（大约在 \(\pm 2/\sqrt{T}\) 处）提供了一个粗略的参考范围，如果估计的自相关值大部分落在此范围内，则表明白噪声的可能性很大。
对于随机游走这种非平稳过程，自相关图会显示所有正滞后处都有持续很高的正自相关，这是非平稳性的一个警示信号。
对于移动平均过程，自相关图可能显示出自相关值在某个滞后（如阶数q）之后突然截尾（变为接近0）。
对于自回归过程，自相关图通常显示出自相关值呈几何（指数）衰减至0，而非线性截尾。

通过观察自相关图，我们可以对数据的特性进行初步探索性分析，例如判断是否平稳、是否存在移动平均或自回归成分等。

检测白噪声

既然我们有了估计自相关函数的方法，就可以用它来检测一个过程是否是白噪声。对于白噪声过程，其理论自相关在 \(h \neq 0\) 时应为0。

在实践中，我们的估计值 \(\hat{\rho}(h)\) 不会恰好为0。我们需要判断它是否“足够小”。可以证明，对于白噪声过程，估计的自相关 \(\hat{\rho}(h) (h \neq 0)\) 近似服从均值为0、方差约为 \(1/T\) 的分布。

因此，一个常用的经验法则是检查自相关图上的蓝虚线，它们大约位于 \(\pm 2/\sqrt{T}\)。如果绝大多数估计的自相关值都落在这个区间内，则没有强烈证据反对白噪声的假设。需要强调的是，这只是一个粗略的定性指导，而非严格的统计假设检验。

时间序列的统计模型

现在，我们开始探讨时间序列的统计建模。首先介绍两个有用的算子定义，它们能简化后续的公式表达。

后移算子 \(B\) 定义为：\(B X_t = X_{t-1}\)。将其迭代应用：\(B^k X_t = X_{t-k}\)。
类似地，有前移算子 \(B^{-1}\)：\(B^{-1} X_t = X_{t+1}\)。

差分算子 \(\nabla\)（或 \(\Delta\)）用于计算序列的变化，类似于离散导数。一阶差分定义为：

\[\nabla X_t = X_t - X_{t-1} = (1 - B) X_t \]

k阶差分可以写为：

\[\nabla^k X_t = (1 - B)^k X_t \]

展开后涉及二项式系数，形式较复杂，但用算子表示非常简洁。

时间序列模型 vs. 线性回归模型

线性回归是理解时间序列模型的一个良好起点。在线性回归中，我们通常有：

\[X_t = \beta_1 Z_{t1} + ... + \beta_p Z_{tp} + \epsilon_t \]

其中 \(\epsilon_t\) 是独立同分布（或至少不相关）的误差项。根据高斯-马尔可夫定理，最小二乘估计量是最优线性无偏估计量。

在时间序列的背景下，我们可以考虑类似的模型，但误差项 \(\epsilon_t\) 被一个平稳过程 \(Y_t\) 所取代：

\[X_t = \beta_1 Z_{t1} + ... + \beta_p Z_{tp} + Y_t \]

这里，\(Y_t\) 不一定是白噪声，它可以是任何平稳过程（如AR或MA过程）。这就引入了时间依赖性。

如果我们忽略这种依赖性，仍然用最小二乘法估计 \(\beta\)，那么得到的残差 \(r_t = X_t - Z_t \hat{\beta}\) 将近似等于 \(Y_t\)，并且可能表现出自相关。关键区别在于：在时间序列分析中，我们将残差过程本身视为一个需要研究的时间序列，分析其自相关结构，而这在线性回归的标准应用中通常被忽略（或假设不存在）。

R语言实战：分析处方药价格数据

让我们应用所学知识分析一个真实数据集：prescript，它包含了美国1986年8月至1992年3月的月度平均处方药价格。

首先，我们绘制数据图。可以观察到明显的上升趋势，这意味着序列的均值是非平稳的。

第一步是尝试移除趋势。我们先用一个简单的线性回归模型（model1）来拟合时间趋势：

model1 <- lm(prescript ~ time(prescript))

模型摘要显示时间变量非常显著，斜率约为2.8，表明每月平均价格上涨约2.8美元。我们将拟合的直线添加到图中。

由于原始数据非平稳，直接计算其自相关函数意义不大。我们更应关注移除趋势后的残差序列。计算并绘制残差序列图，可以看到均值在0附近波动，但似乎存在周期性。

为了更正式地检验平稳性，我们可以使用增广迪基-富勒检验。该检验的原假设是序列非平稳（存在单位根），备择假设是平稳。我们对残差序列进行检验，得到了一个较小的p值，表明在移除线性趋势后，残差序列看起来是平稳的。

既然残差序列是平稳的，我们就可以计算其自相关函数。在ACF图中，我们观察到一个有趣的周期性模式，大约每12个月（一年）出现一个峰值。这强烈暗示数据中存在年度季节性。

我们可以尝试建立一个更复杂的模型（model2）来同时捕捉线性趋势和年度周期性，例如在回归中加入三角函数项：

model2 <- lm(prescript ~ time(prescript) + cos(2*pi*time(prescript)))

模型摘要显示时间和余弦项都是显著的。将拟合曲线添加到原始数据图上，可以看到它更好地捕捉了数据的波动。

然而，对model2的残差进行平稳性检验和观察，发现结果并不理想，残差序列的首尾部分似乎仍有非平稳的迹象。这个例子说明，时间序列建模有时需要反复尝试和判断。对于本数据，或许仅移除线性趋势（model1）后得到的残差序列就是一个合适的、平稳的“噪声”过程，我们可以在此基础上进一步分析其自相关结构（如季节性、ARMA成分等）。

本节课中我们一起学习了如何估计时间序列的自协方差和自相关函数，如何通过自相关图进行探索性分析并检测白噪声，并初步了解了如何将时间序列分解为确定性成分（如趋势、季节性）和平稳的随机成分进行建模分析。在接下来的课程中，我们将深入探讨平滑方法以及更强大的ARMA模型。

004：平滑方法

在本节课中，我们将学习时间序列分析中的平滑方法。这些方法用于减少数据中的噪声，帮助我们更清晰地观察数据的趋势和整体模式。我们将介绍几种不同的平滑技术，并通过R语言演示如何将它们应用于真实数据。

概述

时间序列数据通常包含大量噪声，这使得识别潜在趋势变得困难。平滑方法通过“平均”或“平滑”数据点来减少这种随机波动，从而揭示数据的基本结构。本节课将介绍几种常用的平滑方法，包括移动平均平滑、核平滑、局部加权散点图平滑（LOWESS）和三次平滑样条。这些方法主要用于探索性数据分析，帮助我们直观理解数据。

移动平均平滑器

上一节我们概述了平滑的目的，本节中我们来看看最基础的平滑方法：移动平均平滑器。

移动平均平滑器的核心思想是，通过计算一个滑动窗口内数据点的平均值来平滑序列。给定一个时间序列 X_t（t 从 1 到 T），我们定义平滑后的序列 M_t。

其公式如下：

M_t = Σ_{j=-r}^{r} θ_j * X_{t-j}

其中，r 是带宽参数，决定了窗口的大小（窗口总长度为 2r+1）。系数 θ_j 通常被选择为相等的权重，例如：

θ_j = 1 / (2r + 1)

这表示对窗口内的所有点进行简单算术平均。例如，在分析每日数据时，7日移动平均可以平滑掉单日的异常波动，更清晰地显示周趋势。

在R语言中，可以使用 filter 函数来实现移动平均平滑。

核平滑器

移动平均平滑器为所有窗口内的点分配了相同的权重。核平滑器在此基础上进行了改进，它根据数据点与中心点的距离来分配不同的权重，距离越近的点权重越高。

首先，我们需要定义一个核函数 K_h(x, x_0)。这是一个非负函数，以 x_0 为中心，h 为带宽。函数值随着 |x - x_0| 的增大而减小，并且其积分是有限的。

以下是几种常见的核函数：

高斯核：形状类似钟形曲线。
K_h(x, x_0) ∝ exp( - (x - x_0)^2 / (2h^2) )
三角核：形状为三角形。
K_h(x, x_0) = (1 - |x - x_0|/h)_+
（其中 (·)_+ 表示取正值部分，负值则设为0）
Epanechnikov核：形状为倒抛物线。
K_h(x, x_0) = (1 - (x - x_0)^2 / h^2)_+

那么，如何将核函数用于离散时间序列的平滑呢？核平滑器本质上是一种特殊权重的移动平均平滑器。

其系数 θ_j 由核函数在离散点上的取值归一化后得到：

θ_j = K_h(j, 0) / Σ_{i} K_h(i, 0)

这样，靠近当前时间点 t 的观测值会获得更大的权重，而远处的观测值权重较小。在R中，可以使用 ksmooth 函数进行核平滑。

局部加权散点图平滑（LOWESS）

LOWESS是另一种非常有效的非参数平滑方法。与核平滑类似，它也是局部加权的。

以下是LOWESS方法的核心思想：

它在一个滑动窗口内（例如，围绕点 X_t 的 X_{t-r} 到 X_{t+r}），拟合一个低阶（通常是线性或二次）多项式回归。拟合时，它会根据数据点与中心点的距离分配权重。然而，LOWESS算法的具体实现细节（如权重的精确计算和迭代过程）在R的文档中描述得比较模糊，被形容为一个“复杂算法”。

尽管如此，它仍然是一个强大且常用的探索性工具。在R中，可以直接使用 lowess 函数。

三次平滑样条

最后，我们介绍一种更结构化的平滑方法：三次平滑样条。它通过拟合一个全局的、分段的三次多项式来平滑数据，并加入一个惩罚项来控制平滑度。

首先，理解三次样条。我们将时间区间 [1, T] 划分为 K 个子区间，分界点称为“节点”。在每个子区间内，我们拟合一个三次多项式：

M_i(t) = β_{i0} + β_{i1}t + β_{i2}t^2 + β_{i3}t^3

为了确保曲线整体光滑，我们通常要求这些分段多项式在节点处连续且具有连续的一阶和二阶导数。

如果仅仅是这样，我们就是在做分段多项式回归。而三次平滑样条在此基础上，最小化一个包含惩罚项的目标函数：

最小化： Σ_{t=1}^T (X_t - M(t))^2 + λ ∫ [M''(s)]^2 ds

第一项是常见的残差平方和，衡量拟合程度。
第二项是惩罚项，它惩罚曲线的曲率（即二阶导数 M''）。λ 是平滑参数。
- 当 λ = 0 时，惩罚项无效，结果就是普通的最小二乘拟合，可能非常“崎岖”。
- 当 λ → ∞ 时，惩罚项占主导，迫使二阶导数处处趋近于0，结果就是一条直线（最平滑的状态）。

因此，通过调节 λ，我们可以在拟合优度和平滑度之间进行权衡。这类似于线性回归中的岭回归思想，都是通过添加惩罚项来约束模型复杂度。

在R中，可以使用 smooth.spline 函数来拟合三次平滑样条。

R语言应用演示

现在，我们将这些平滑方法应用于一个真实的时间序列数据集：美国处方药成本数据（1980年代中期至1990年代初期）。

以下是应用不同平滑器的步骤和效果：

移动平均平滑：使用 filter 函数。例如，7个月和23个月的移动平均。较长的窗口会产生更平滑但也更滞后的曲线，且会在序列两端产生缺失值。
核平滑：使用 ksmooth 函数，选择高斯核。通过调整带宽参数 bandwidth，可以控制平滑程度。带宽越小，曲线越贴近原始数据；带宽越大，曲线越平滑。
LOWESS平滑：使用 lowess 函数，调整参数 f（平滑跨度）。f 值越小，曲线越贴近数据细节；f 值越大，曲线越平滑。
三次平滑样条：使用 smooth.spline 函数，调整参数 spar（平滑参数）。spar 越大，曲线越平滑。

所有方法都揭示了数据中明显的上升趋势和周期性波动。尽管实现方式不同，但在选择合适的参数后，它们往往能给出视觉上相似的趋势估计。

重要注意事项：平滑的用途与局限

在结束之前，必须强调平滑方法的主要用途和局限。

平滑是一种强大的探索性数据分析工具。它帮助我们可视化趋势、发现周期性和异常点，为后续的正式建模提供直觉。例如，对处方药数据的一阶差分序列进行平滑，可以更清楚地看到其季节性成分。

然而，我们必须谨慎：

平滑后的数据不应直接用于严格的统计推断或模型拟合。 因为平滑过程移除了数据中的“噪声”，而许多时间序列模型（如ARIMA）正是利用这些噪声的结构（如自相关）进行建模和预测。如果对数据平滑后再进行自相关分析或回归，可能会得到误导性的结果，例如发现不存在的相关性或错误判断平稳性。

正确的做法是：使用平滑来探索数据、形成假设（例如，“数据可能存在12个月的季节性”），然后基于原始数据（或适当的差分后数据）来建立和检验正式的统计模型。

总结

本节课中我们一起学习了几种重要的时间序列平滑方法：

移动平均平滑器：通过窗口内简单平均来平滑。
核平滑器：根据距离加权平均，提供了更灵活的平滑方式。
LOWESS：一种局部加权多项式回归平滑方法。
三次平滑样条：通过带惩罚项的分段多项式拟合实现全局平滑。

我们了解了它们的原理、在R中的实现，并认识到它们主要是服务于探索性数据分析的视觉辅助工具。从下一讲开始，我们将进入课程的核心部分：学习用于预测和推断的正式时间序列模型，如ARMA模型。

005：AR与MA模型理论

在本节课中，我们将深入学习时间序列分析的核心模型：自回归模型和移动平均模型。我们将详细探讨它们的数学性质、因果性与可逆性，以及这些概念在统计建模中的重要性。

概述

自回归移动平均模型是时间序列建模的基石。为了理解如何用这些模型拟合数据，我们首先需要深入理解这些随机过程本身。本节将分别剖析AR和MA模型，为后续学习ARMA和ARIMA模型打下坚实基础。

自回归模型

自回归模型的核心思想是，当前时刻的观测值可以表示为过去若干个观测值的线性组合，再加上一个白噪声项。

一个 p阶自回归过程 可以定义为：
Xt = φ1 * Xt-1 + φ2 * Xt-2 + ... + φp * Xt-p + Wt
其中，Wt 是均值为0、方差为 σ² 的白噪声，φ1, ..., φp 是自回归系数。

为了书写方便，我们常使用后移算子B（B * Xt = Xt-1）和自回归算子来简化表达式。定义自回归算子 Φ(B) = 1 - φ1B - φ2B² - ... - φpB^p，则AR(p)模型可以简洁地写为：
Φ(B) * Xt = Wt

在分析中，我们通常假设过程均值为零。如果过程有非零均值μ，我们可以通过定义 X̃t = Xt - μ 将其转化为零均值过程，这不会影响模型的核心性质。

因果性

因果性是一个直观的概念：过去影响未来，但未来不影响过去。在时间序列中，我们希望模型是因果的。

考虑一个简单的AR(1)过程：Xt = φ1 * Xt-1 + Wt。

当 |φ1| < 1 时，该过程是平稳的，并且可以表示为仅依赖于过去白噪声的线性过程（即因果表示）：
Xt = Σ (i=0 to ∞) φ1^i * Wt-i
当 |φ1| > 1 时，过程在数学上仍然是平稳的，但其表示形式会依赖于未来的白噪声，成为非因果过程。
当 |φ1| = 1 时，过程变为随机游走，是非平稳的。

一个关键且有趣的现象是：对于一个非因果的AR(1)过程（|φ1| > 1），如果白噪声是高斯分布的，那么总存在一个因果的AR(1)过程（系数为 1/φ1），它与原过程具有完全相同的概率分布（即相同的均值和自协方差函数）。

这对统计学意味着可识别性问题。如果我们拿到一组数据，我们无法区分它来自系数为 φ 的非因果模型，还是来自系数为 1/φ 的因果模型。因此，在实践中，我们总是选择并估计因果的AR模型，这符合我们对“过去影响未来”的认知。

自回归算子的根

为了将AR(p)模型表示为因果的线性过程，我们需要研究自回归算子 Φ(B)。

我们将算子中的 B 替换为复数 z，得到多项式 Φ(z) = 1 - φ1*z - φ2*z² - ... - φp*z^p。根据代数基本定理，这个p次多项式可以分解为：
Φ(z) = (-1)^p * φp * (z - r1)(z - r2)...(z - rp)
其中 r1, r2, ..., rp 是该多项式的根（可能为复数）。

核心结论：当且仅当所有根的模长都大于1（即所有根都在复平面上的单位圆外）时，AR(p)过程是因果的，并且可以唯一地表示为一个仅依赖于过去白噪声的线性过程。这就是许多平稳性检验（如ADF检验）中“检验单位根”的数学来源——我们正是要避免根落在单位圆上或圆内。

移动平均模型

移动平均模型描述的是，当前观测值受到当前及过去若干期白噪声冲击的线性影响。

一个 q阶移动平均过程 定义为：
Xt = Wt + θ1 * Wt-1 + θ2 * Wt-2 + ... + θq * Wt-q
其中，Wt 是白噪声，θ1, ..., θq 是移动平均系数。

同样，我们可以定义移动平均算子 Θ(B) = 1 + θ1B + θ2B² + ... + θqB^q，从而将模型写为：
Xt = Θ(B) * Wt

与AR模型不同，MA模型对于任意系数都是平稳的。其自协方差函数在滞后 q 步之后截尾（变为0）。

可逆性

MA模型也存在类似AR模型因果性的可逆性问题。考虑一个MA(1)过程：Xt = Wt + θ1 * Wt-1。

我们可以证明，存在另一个MA(1)过程 Yt = Vt + (1/θ1) * Vt-1（其中 Vt 的方差经过调整），它与原过程 Xt 具有完全相同的自协方差结构。

这就产生了统计估计上的模糊性：我们估计的是 θ1 还是 1/θ1？为了解决这个问题，我们引入可逆性概念。如果一个MA过程可以“反转”，表示成当前观测值关于其过去观测值的无限加权和（即 Wt = Σ πi * Xt-i），则称它是可逆的。

类似于AR模型的因果性条件，MA(q)过程可逆的充要条件是：其移动平均算子 Θ(z) 的所有根的模长都大于1（位于单位圆外）。在实践中，我们总是选择并估计满足可逆性条件的MA模型，这通常也对应着 |θ1| < 1 的情形，意味着过去冲击的影响随时间衰减，这符合直观。

总结

本节课我们一起深入探讨了时间序列分析的两大核心组件：自回归模型和移动平均模型。

我们学习了AR和MA模型的定义与算子表示。
我们深入理解了AR模型的因果性概念及其与自回归算子根的位置（在单位圆外）的等价关系。统计建模中我们只考虑因果模型。
我们探讨了MA模型的可逆性概念，其同样与移动平均算子根的位置（在单位圆外）等价。统计建模中我们只考虑可逆模型。
我们看到了在特定条件下（如高斯白噪声），可能存在不同的模型参数生成完全相同的观测数据分布，这引出了统计模型的可识别性问题，而坚持因果性和可逆性是解决该问题的标准做法。

掌握这些理论基础至关重要。在下一讲中，我们将把AR和MA模型结合起来，构成更强大的ARMA模型，并继续探讨其因果性与可逆性条件。

006：ARMA与ARIMA模型 🧮

在本节课中，我们将继续学习时间序列的理论知识。上一讲我们分别介绍了自回归（AR）过程和移动平均（MA）过程。本节中，我们将把它们结合起来，探讨自回归移动平均（ARMA）模型，并进一步介绍自回归积分移动平均（ARIMA）模型。掌握这些理论工具，将为我们后续进行模型估计和拟合数据打下坚实的基础。

回顾：AR与MA过程 📚

在深入ARMA模型之前，我们先简要回顾一下上一讲的核心内容。

自回归（AR）过程

我们讨论了AR(p)过程，即p阶自回归过程。其定义如下：

公式：Φ(B) X_t = W_t

其中，Φ(B)是自回归算子，W_t是白噪声。我们得出的关键结论是：如果多项式Φ(z)的所有根r_i都位于单位圆外（即|r_i| > 1），则该AR(p)过程是因果的。这意味着X_t可以表示为一个因果线性过程，即仅依赖于当前和过去的白噪声。

移动平均（MA）过程

接着，我们讨论了MA(q)过程，即q阶移动平均过程。其定义如下：

公式：X_t = Θ(B) W_t

其中，Θ(B)是移动平均算子。类似地，如果多项式Θ(z)的所有根都位于单位圆外，则该MA(q)过程是可逆的。这意味着白噪声W_t可以表示为X_t过去值的线性组合。

结合：ARMA模型 ⚙️

上一节我们介绍了AR和MA过程，本节中我们来看看如何将它们融合成一个更强大的模型。

ARMA(p, q)模型，即自回归移动平均模型，结合了AR(p)和MA(q)的特点。其定义如下：

公式：Φ(B) X_t = Θ(B) W_t

或者，等价地写为：
X_t = φ_1 X_{t-1} + ... + φ_p X_{t-p} + W_t + θ_1 W_{t-1} + ... + θ_q W_{t-q}

其中，系数φ_1, ..., φ_p和θ_1, ..., θ_q是实数，且φ_p ≠ 0，θ_q ≠ 0，以确保模型的阶数准确。

模型的可辨识性与简化

在构建ARMA模型时，需要注意一个关键问题：非唯一性。有时，一个ARMA模型可能被过度参数化，可以通过约简得到更简单的形式。

示例：考虑过程 X_t = θ X_{t-1} - θ W_{t-1} + W_t。

初看之下，它像是一个ARMA(1,1)模型。
但我们可以将其重写为 (1 - θB) X_t = (1 - θB) W_t。
如果|θ| < 1，我们可以“约去”两边的(1 - θB)算子，得到 X_t = W_t。
这意味着，该过程本质上就是白噪声，而非真正的ARMA(1,1)模型。

这个例子说明，如果自回归多项式Φ(z)和移动平均多项式Θ(z)有公因子（即共同的根），那么模型就可以被简化。因此，在标准设定下，我们假设Φ(z)和Θ(z)是互质的（没有公因子），以确保模型是唯一且最简形式的。

因果性与可逆性

对于一个ARMA(p, q)过程，我们关心两个核心性质：

因果性：如果Φ(z)的根都在单位圆外，则该ARMA过程是因果的。这意味着X_t可以表示为当前和过去白噪声的线性组合（一个因果线性过程）：
公式：X_t = Σ_{j=0}^{∞} ψ_j W_{t-j}，其中系数ψ_j是绝对可和的。
可逆性：如果Θ(z)的根都在单位圆外，则该ARMA过程是可逆的。这意味着白噪声W_t可以表示为当前和过去X_t的线性组合：
公式：W_t = Σ_{j=0}^{∞} π_j X_{t-j}，其中系数π_j是绝对可和的。

在实际建模中，我们通常希望并且会寻找具有因果性和可逆性的ARMA模型表示，因为这样的表示在统计估计中更稳定、更易于处理。

扩展：ARIMA模型 📈

在现实数据中，我们经常遇到的情况是，序列Y_t本身并非平稳的ARMA过程，但经过差分处理后，可以转化为平稳的ARMA过程。这就引出了ARIMA模型。

差分算子

首先，定义d阶差分算子∇^d：
公式：∇^d = (1 - B)^d
其中，B是后移算子（B X_t = X_{t-1}）。一阶差分∇ X_t = X_t - X_{t-1}常用于消除线性趋势。

ARIMA模型定义

我们说一个序列X_t服从ARIMA(p, d, q)过程，如果它的d阶差分是一个ARMA(p, q)过程。

公式：Φ(B) (1 - B)^d X_t = Θ(B) W_t

其中：

p：差分后序列的自回归阶数。
d：使序列变得平稳所需的最小差分阶数。
q：差分后序列的移动平均阶数。

直观理解：许多时间序列数据（如股价、气温）可能包含确定性的趋势（如线性增长）。ARIMA模型通过差分运算(1-B)^d先移除这些趋势，然后对得到的平稳序列拟合ARMA模型。例如，如果数据有线性趋势，通常d=1就足够了；如果有二次趋势，可能需要d=2。

包含常数项的ARIMA

有时，即使经过d阶差分，序列的均值可能仍不为零。更一般的ARIMA模型可以包含一个常数项μ：

公式：Φ(B) (1 - B)^d X_t = δ + Θ(B) W_t
其中，δ = μ (1 - φ_1 - ... - φ_p)。

总结 🎯

本节课中我们一起学习了时间序列建模中两个核心的模型家族。

ARMA(p, q)模型：融合了自回归和移动平均成分，用于描述平稳时间序列。其关键性质是因果性和可逆性，这取决于相应多项式的根是否在单位圆外。
ARIMA(p, d, q)模型：ARMA模型的扩展，通过引入d阶差分来处理非平稳数据（特别是含有确定性趋势的数据）。它先对数据差分使其平稳，再应用ARMA模型。

这些模型为我们提供了描述各种时间序列动态结构的强大数学框架。在接下来的课程中，我们将学习如何利用自相关和偏自相关函数来识别模型的阶数(p, d, q)，并最终进入统计部分——如何从实际数据中估计这些模型的参数。

007：自相关与平稳性

在本节课中，我们将要学习如何检验时间序列数据中的自相关性以及判断序列是否平稳。这是进行有效统计推断的基础。

在之前的几讲中，我们深入探讨了AR、MA、ARMA和ARIMA模型的数学性质。现在，我们将暂时从繁重的理论中抽身，转而关注对这些模型的实际统计检验。具体来说，今天的课程有两个核心主题：自相关和平稳性。

检验自相关性

上一节我们介绍了ARMA模型的性质，本节中我们来看看如何检验数据中是否存在显著的自相关。除非数据是白噪声，否则我们通常会发现一些自相关。我们需要从统计上检验这些自相关是否显著不为零。

以下是两种常用的检验方法：

Box-Pierce 与 Ljung-Box 检验

这两种检验旨在回答同一个问题：在滞后阶数1到H之间，是否存在非零的自相关？其原假设是：所有滞后阶数（1到H）的自相关系数均为零。

Box-Pierce检验：其检验统计量 Q_BP = n * Σ_{j=1}^{H} (ρ̂(j))^2，在原假设下，它近似服从自由度为H的卡方分布。
Ljung-Box检验：这是对Box-Pierce检验的改进，尤其适用于有限样本。其检验统计量为 Q_LB = n(n+2) * Σ_{j=1}^{H} (ρ̂(j))^2 / (n - j)，同样近似服从自由度为H的卡方分布。

在实践中，我们通常不会直接将检验应用于原始时间序列，而是应用于拟合某个ARMA(p, q)模型后得到的残差序列。如果模型拟合良好，残差应近似为白噪声（即无自相关）。此时，检验统计量的自由度需要调整为 H - p - q。在R语言中，可以使用 Box.test() 函数进行这些检验。

Durbin-Watson 检验

Durbin-Watson检验专门用于检验一阶自相关。它通常应用于线性回归模型的残差。

该检验的统计量定义为：
Q_DW = Σ_{t=2}^{T} (r̂_t - r̂_{t-1})^2 / Σ_{t=1}^{T} (r̂_t)^2
其中，r̂_t 是回归残差。

Q_DW ≈ 2 表示无自相关。
Q_DW 显著接近 0 表示存在正自相关。
Q_DW 显著接近 4 表示存在负自相关。

其P值的计算涉及卡方随机变量的线性组合，计算较为复杂，但R语言中的 dwtest() 函数（来自 lmtest 包）可以方便地完成。

Breusch-Godfrey 检验

Breusch-Godfrey检验（BG检验）是Durbin-Watson检验的推广，它可以检验残差中是否存在高阶（P阶）的自回归过程。该检验通过比较 n * R²（其中R²是残差对其自身滞后项回归的判定系数）与自由度为P的卡方分布来进行。在R中，可以使用 bgtest() 函数。

检验平稳性

在学习了如何检验自相关之后，我们接下来探讨时间序列分析中另一个基石性的概念：平稳性。平稳性对于参数估计和统计推断至关重要。以下是两种常见的单位根检验方法，用于判断序列是否平稳。

(Augmented) Dickey-Fuller 检验

Dickey-Fuller检验及其增强版本（ADF检验）是典型的单位根检验。其核心思想是检验自回归多项式是否存在模为1的根（单位根），若存在，则序列非平稳。

基本DF检验：针对AR(1)过程 X_t = φ X_{t-1} + w_t。通过检验一阶差分方程 ∇X_t = φ' X_{t-1} + w_t 中的系数 φ' 是否为零（等价于检验原序列的φ是否为1）。原假设 H0: φ' = 0 (序列非平稳)，备择假设 H1: φ' < 0 (序列平稳)。注意，这里显著的P值意味着我们拒绝“非平稳”的原假设，支持序列是平稳的。
增强型ADF检验：将上述思想推广到AR(p)过程。检验方程扩展为：
∇X_t = φ'1 X + Σ_{i=1}^{p-1} φ'{i+1} ∇X + w_t
检验 φ'_1 是否为零，其假设与DF检验相同。R语言中的 adf.test() 函数（来自 tseries 包）在拟合时会先移除数据的线性趋势，然后再进行单位根检验。

Phillips-Perron 检验

Phillips-Perron检验是另一种单位根检验方法。它与ADF检验的目的相同，但试图对回归假设的偏离（如异方差和误差相关）更加稳健。它使用了Newey-West估计量来校正标准误。在R语言中，可以使用 pp.test() 函数（同样来自 tseries 包）。

总结

本节课中我们一起学习了时间序列分析中两个核心的检验主题：自相关与平稳性。

在自相关检验部分，我们介绍了 Box-Pierce、Ljung-Box、Durbin-Watson 和 Breusch-Godfrey 检验，它们主要用于检测序列（尤其是模型残差）中是否存在显著的自相关结构。
在平稳性检验部分，我们探讨了 （增强型）Dickey-Fuller 检验和 Phillips-Perron 检验这两种单位根检验方法，它们用于判断一个时间序列是否具有平稳性。

理解并正确应用这些检验，是构建可靠时间序列模型、进行有效预测和推断的关键步骤。在接下来的实践中，我们将有机会在R语言中应用这些检验函数来分析真实或模拟数据。

008：在R中检验平稳性与自相关函数

在本节课中，我们将学习如何在R软件中实际操作，检验时间序列数据的自相关性和平稳性。我们将使用模拟数据和真实数据，介绍几种关键的统计检验方法，并学习如何解读结果。

上一讲我们主要介绍了自相关检验和平稳性检验的理论知识。本节中，我们来看看如何在R中应用这些理论。

准备工作：生成模拟数据

首先，我们需要生成一些模拟的时间序列数据，以便在已知“正确答案”的情况下进行测试。

以下是生成三种不同类型时间序列的R代码：

# 生成500个时间点的白噪声
WN <- rnorm(500)

# 生成一个自回归系数为0.7的AR(1)过程
AR1 <- arima.sim(list(order = c(1,0,0), ar = 0.7), n = 500, method = "recursive")

# 生成一个更复杂的AR(4)过程
AR4 <- arima.sim(list(order = c(4,0,0), ar = c(0.7, -0.4, 0.2, -0.1)), n = 500)

模型拟合简介

在深入检验之前，我们先简要了解如何为时间序列数据拟合模型。这有助于我们后续分析模型的残差。

以下是使用R函数拟合AR模型的示例：

# 使用ar()函数自动拟合AR模型
ar(WN) # 拟合白噪声，应识别为近似阶数为0
ar(AR1) # 估计AR(1)系数，应接近0.7
ar(AR4) # 估计AR(4)系数，应接近我们设定的参数

# 使用arima()函数手动指定模型阶数进行拟合
# 为白噪声拟合一个不合适的ARMA(1,1)模型
arima(WN, order = c(1,0,1))

自相关检验

接下来，我们使用几种统计检验来探测时间序列中是否存在显著的自相关。

Box-Pierce 与 Ljung-Box 检验

Box.test函数可以执行Box-Pierce检验或Ljung-Box检验。其原假设是：在指定的滞后阶数内，所有自相关系数均为零。

以下是应用该检验的示例：

# 对白噪声序列进行检验（滞后阶数=1）
Box.test(WN, lag = 1, type = "Box-Pierce") # 预期p值不显著
Box.test(WN, lag = 1, type = "Ljung-Box")  # 预期p值不显著

# 对AR(1)序列进行检验
Box.test(AR1, lag = 1, type = "Box-Pierce") # 预期p值非常显著

# 对一个系数很小的AR(1)序列进行检验（难以检测）
AR1_small <- arima.sim(list(order = c(1,0,0), ar = 0.05), n = 500)
Box.test(AR1_small, lag = 1, type = "Ljung-Box") # 可能无法检测出显著性

# 对一个特殊的AR(2)序列进行检验（φ1=0， φ2≠0）
AR2_special <- arima.sim(list(order = c(2,0,0), ar = c(0, 0.7)), n = 500)
Box.test(AR2_special, lag = 1, type = "Ljung-Box")  # 滞后1阶不显著
Box.test(AR2_special, lag = 2, type = "Ljung-Box")  # 滞后2阶显著

重要提示：Box检验常应用于拟合模型后的残差。此时，需要从自由度中减去已估计的参数个数（使用fitdf参数）。

Durbin-Watson 检验

Durbin-Watson检验专门用于检测残差中是否存在一阶自相关。需要加载lmtest包。

以下是应用该检验的示例：

library(lmtest)
# 检验白噪声
dwtest(WN ~ 1) # 公式`~1`表示只拟合截距项，预期不显著

# 检验AR(1)过程
dwtest(AR1 ~ 1) # 预期非常显著

# 检验真实数据（处方药成本数据）及其线性趋势拟合后的残差
data(prescript) # 假设数据已加载
dwtest(prescript ~ time(prescript)) # 预期残差中仍存在显著自相关

Breusch-Godfrey 检验

Breusch-Godfrey检验是Durbin-Watson检验的推广，可以检验更高阶的序列相关。

以下是应用该检验的示例：

# 加载所需的包（如果未加载）
library(lmtest)
# 检验白噪声，滞后阶数设为6
bgtest(WN ~ 1, order = 6) # 预期不显著

# 检验AR(1)过程
bgtest(AR1 ~ 1, order = 1) # 预期显著

# 检验特殊的AR(2)过程，展示其可检测高阶相关的优势
bgtest(AR2_special ~ 1, order = 1) # 滞后1阶不显著
bgtest(AR2_special ~ 1, order = 2) # 滞后2阶显著

平稳性检验（单位根检验）

平稳性检验的核心是判断时间序列是否存在单位根。存在单位根意味着序列是非平稳的。

Augmented Dickey-Fuller (ADF) 检验

ADF检验的原假设是序列存在单位根（非平稳），备择假设是序列平稳。

以下是应用ADF检验的示例：

library(tseries)
# 检验白噪声（500个观测值）
adf.test(WN) # 预期p值很小，拒绝原假设，认为平稳

# 检验观测值较少时的情况
adf.test(WN[1:50])  # p值可能仍较小
adf.test(WN[1:20])  # p值可能变得不显著，数据不足难以判断

# 检验随机游走（非平稳过程的典型例子）
RW <- cumsum(rnorm(500)) # 累积求和生成随机游走
adf.test(RW) # 预期p值很大，无法拒绝非平稳的原假设

# 检验处方药成本数据
adf.test(prescript)
# 注意：adf.test默认会先尝试移除常数项和线性趋势，再对残差进行检验。
# 因此，即使原始序列有明显趋势，结果也可能显示残差序列是平稳的。

Phillips-Perron (PP) 检验

PP检验是另一种单位根检验方法，与ADF检验类似，但在处理序列相关时采用了不同的修正方法。

以下是应用PP检验的示例：

# 检验处方药成本数据
pp.test(prescript) # 结果应与ADF检验类似

# 检验白噪声
pp.test(WN) # 预期p值很大（不显著）

# 比较ADF和PP检验在小样本下的表现
pp.test(WN[1:20]) # 可能与adf.test(WN[1:20])的结果有所不同

自相关函数与偏自相关函数预览

在结束前，我们预览一下下一讲的重要内容：自相关函数和偏自相关函数图。它们能直观地提示时间序列模型的可能阶数。

以下是绘制ACF和PACF图的示例：

# 绘制AR(1)过程的ACF和PACF图
acf(AR1, main = “AR(1) ACF”)
pacf(AR1, main = “AR(1) PACF”) # PACF在滞后1阶后截尾，提示AR(1)

# 绘制AR(4)过程的ACF和PACF图
acf(AR4, main = “AR(4) ACF”)
pacf(AR4, main = “AR(4) PACF”) # PACF在滞后4阶后截尾，提示AR(4)

# 绘制特殊AR(2)过程的PACF图
pacf(AR2_special, main = “Special AR(2) PACF”) # PACF在滞后2阶后截尾

提示：偏自相关函数在滞后p阶后出现“截尾”现象，通常暗示这是一个p阶自回归过程。我们将在下一讲详细讨论。

总结

本节课中，我们一起学习了在R环境中进行时间序列分析的关键实操步骤：

我们学习了如何生成模拟的AR、MA等时间序列数据。
我们介绍了Box.test、dwtest、bgtest等多种自相关检验方法，并理解了它们在不同场景下的应用。
我们掌握了adf.test和pp.test这两种单位根检验方法，用于判断时间序列的平稳性。
我们了解到，许多检验（如ADF）在应用前会默认消除数据中的确定性趋势。
我们预览了自相关图和偏自相关图，认识到它们是识别模型类型的强大可视化工具。

通过结合模拟数据和真实数据练习这些检验，我们能够更好地理解统计检验的效力、局限性以及在实际分析中需要注意的细节（如自由度调整、样本量影响等）。这些工具为我们后续建立和诊断时间序列模型奠定了坚实的基础。

009：自相关与偏自相关函数（第一部分）

在本节课中，我们将学习时间序列分析中的两个核心工具：自相关函数（ACF）和偏自相关函数（PACF）。我们将从回顾自相关函数开始，然后深入探讨其“兄弟”——偏自相关函数的概念、理论性质及其在分析时间序列数据中的应用。

从相关到偏相关

上一节我们介绍了自相关函数，本节中我们来看看一个更复杂但密切相关的概念：偏相关。为了理解偏自相关，我们首先需要理解在一般随机变量背景下的偏相关概念。

什么是偏相关？

考虑两个实值随机变量 X 和 Y。它们之间的相关系数衡量的是它们之间的线性关系强度。然而，有时 X 和 Y 看起来相关，可能是因为它们都与第三个变量 Z 有关。

例如，假设：

X 代表房屋价格。
Y 代表汽车价格。
Z 代表收入。

房屋价格和汽车价格可能呈现正相关，但这很可能是因为它们都受到收入水平的影响。如果我们“控制”或“给定”收入 Z，那么房屋价格和汽车价格之间的相关性可能会消失。这种在给定 Z 的条件下，X 和 Y 之间的相关性，就是偏相关。

如何计算偏相关？

偏相关可以通过回归分析来计算。以下是计算步骤：

分别用 Z 对 X 和 Y 进行线性回归：
- X_hat = α_0 + α_1 * Z
- Y_hat = β_0 + β_1 * Z
计算两个回归的残差：
- Residual_X = X - X_hat
- Residual_Y = Y - Y_hat
偏相关系数就是这两个残差序列之间的普通相关系数。

核心公式：
偏相关系数 ρ_{XY·Z} = Corr( X - E[X|Z], Y - E[Y|Z] )

如果偏相关系数为零（或接近零），则表明在给定 Z 的条件下，X 和 Y 不相关。如果数据服从联合正态分布，零偏相关还意味着条件独立性，这是一个更强的结论。

回到时间序列：偏自相关函数（PACF）

现在，我们将偏相关的概念应用到时间序列中，定义偏自相关函数。

对于一个平稳时间序列 {X_t}，在滞后 h 处的偏自相关系数，记作 φ_{hh}，其定义如下：

当 h = 1 时，偏自相关系数就等于普通自相关系数：φ_{11} = ρ(1)。
当 h > 1 时，φ_{hh} 是 X_{t+h} 和 X_t 之间，在给定中间所有观测值 (X_{t+1}, X_{t+2}, ..., X_{t+h-1}) 的条件下的相关系数。

计算理解：
这相当于进行两次回归：

用 X_{t+1}, ..., X_{t+h-1} 预测 X_{t+h}，得到预测值 X_{t+h} 和残差 R_{t+h}。
用 X_{t+1}, ..., X_{t+h-1} 预测 X_t，得到预测值 X_t 和残差 R_t。
偏自相关系数 φ_{hh} 就是这两个残差序列 R_{t+h} 和 R_t 之间的相关系数。

直观意义：
PACF 试图揭示在排除中间所有滞后项的影响后，X_t 和 X_{t+h} 之间是否存在“直接”的关系。这对于识别自回归过程的阶数至关重要。

自回归过程的ACF特征

接下来，我们探讨不同类型时间序列过程中ACF的表现。首先从自回归过程开始。

AR(1) 过程的 ACF

考虑一个平稳、零均值的 AR(1) 过程：
X_t = φ X_{t-1} + W_t，其中 |φ| < 1。

通过计算，我们可以得到其自相关函数满足一个一阶差分方程，并具有简单的几何衰减形式：

核心公式：
ρ(h) = φ^h，对于 h ≥ 0。

这意味着 AR(1) 过程的 ACF 图从 ρ(0)=1 开始，在滞后 1 处为 φ，然后按 φ^2, φ^3, ... 的方式指数衰减到零。如果 φ 为正，则单调衰减；如果 φ 为负，则呈现正负交替的衰减模式。

AR(2) 及更高阶 AR(p) 过程的 ACF

对于 AR(2) 过程：X_t = φ_1 X_{t-1} + φ_2 X_{t-2} + W_t，其 ACF 满足二阶差分方程：ρ(h) = φ_1 ρ(h-1) + φ_2 ρ(h-2)。

解此方程得到的 ACF 形式取决于其特征多项式的根：

两个不等的实根：ACF 是两个几何衰减项的线性组合，整体仍呈现几何衰减趋势。
两个相等的实根（重根）：ACF 是 (C1 + C2 * h) * (几何衰减项)，衰减模式可能被一个线性项修饰。
一对共轭复根：ACF 的形式为 (几何衰减项) * cos(θh + 常数)。这意味着 ACF 在几何衰减的同时，还会产生阻尼振荡（正弦波式的波动），振荡频率由复根的幅角决定。

AR(p) 过程的一般结论：
无论特征根是实是复，是否有重根，AR(p) 过程的 ACF 通常表现为混合的几何衰减（可能伴有阻尼振荡），并且不会在有限滞后后突然截断，而是逐渐衰减至零。

滑动平均过程的ACF特征

现在，我们来看滑动平均过程 ACF 的特征，这与 AR 过程形成鲜明对比。

考虑一个 MA(q) 过程：
X_t = W_t + θ_1 W_{t-1} + ... + θ_q W_{t-q}，其中 θ_0 = 1。

其自协方差和自相关函数有一个非常简洁的性质：

核心公式：
对于 h > q，有 γ(h) = 0 且 ρ(h) = 0。

推导结果：
ρ(h) = (∑_{j=0}^{q-h} θ_j θ_{j+h}) / (∑_{j=0}^{q} θ_j^2)，当 0 ≤ h ≤ q。
ρ(h) = 0，当 h > q。

重要特征：
MA(q) 过程的 ACF 在滞后 q 处之后突然截断（严格为零）。这是一个非常关键的性质。

直观意义：
因为 MA 过程只依赖于最近 q 期的噪声，所以相隔超过 q 个时间点的观测值之间没有线性依赖关系。因此，我们可以通过观察 ACF 图在哪个滞后之后变得不显著（接近零），来初步判断 MA 过程的阶数 q。

例如，一个 MA(3) 过程的 ACF 可能在滞后 1、2、3 处显著非零，但从滞后 4 开始全部接近零。

总结与预告

本节课中我们一起学习了：

偏相关与偏自相关：理解了在控制其他变量影响后，两个变量之间“净”相关性的概念，并将其引入时间序列定义了偏自相关函数。
AR过程的ACF：AR过程的ACF通常表现为几何衰减（可能伴有振荡），且不会突然截断。衰减速度由自回归系数（或特征根）决定。
MA过程的ACF：MA(q)过程的ACF有一个标志性特征，即在滞后 q 之后突然截断为零。这可用于初步识别MA过程的阶数。

这些特征为我们识别时间序列模型提供了重要线索。在下节课中，我们将重点研究偏自相关函数的表现。你会发现一个有趣的对偶现象：对于AR过程，其PACF会在滞后超过阶数p后截断；而对于MA过程，其PACF则表现为几何衰减。结合ACF和PACF的图形，我们将能更有效地在AR、MA或ARMA模型之间做出初步判断。

010：ACF与PACF（第二部分）

在本节课中，我们将继续学习自相关函数（ACF）和偏自相关函数（PACF）。上一讲我们介绍了ACF和PACF的基本概念，本节我们将深入探讨PACF，特别是它在AR和MA过程中的表现，并学习如何利用这些工具来识别时间序列的模型类型和阶数。

回顾：偏自相关函数（PACF）的定义

上一节我们介绍了偏自相关的概念。为了确保概念清晰，我们在此简要回顾一下。

对于一个平稳时间序列 (X_t)，在滞后阶数 (h) 处的偏自相关函数 (\phi_{xx}(h)) 定义如下：

当 (h=1) 时，(\phi_{xx}(1)) 就是滞后1的自相关系数。
当 (h>1) 时，(\phi_{xx}(h)) 是 (X_{t+h}) 和 (X_t) 在剔除了中间观测值 (X_{t+1}, X_{t+2}, ..., X_{t+h-1}) 的线性影响后的相关系数。

更具体地说，我们通过将 (X_t) 和 (X_{t+h}) 分别对中间观测值进行回归，得到拟合值 (\hat{X}t) 和 (\hat{X})，然后计算残差 ((X_t - \hat{X}t)) 与 ((X - \hat{X}_{t+h})) 之间的相关系数。

如果偏自相关为零，意味着 (X_t) 和 (X_{t+h}) 之间的关系完全由中间的时间点所解释，它们之间没有直接的关联。

AR(1)过程的PACF

为了建立直观理解，我们先从最简单的AR(1)过程开始分析。

考虑一个平稳、因果的AR(1)过程：
[
X_t = \phi X_{t-1} + W_t, \quad |\phi| < 1
]
其中 (W_t) 是白噪声。

计算滞后2的偏自相关

根据定义，我们需要计算 (X_t) 和 (X_{t-2}) 在剔除 (X_{t-1}) 影响后的相关性。这等价于计算以下两个残差之间的相关性：

(X_t - \hat{X}_t)，其中 (\hat{X}t) 是基于 (X) 对 (X_t) 的最佳线性预测。
(X_{t-2} - \hat{X}{t-2})，其中 (\hat{X}) 是基于 (X_{t-1}) 对 (X_{t-2}) 的最佳线性预测。

对于AR(1)过程，可以证明，基于 (X_{t-1}) 预测 (X_t) 和预测 (X_{t-2}) 的最佳线性系数都是 (\phi)。因此：

(X_t - \hat{X}t = X_t - \phi X = W_t)
(X_{t-2} - \hat{X}{t-2} = X - \phi X_{t-1})

现在计算这两个残差的协方差。由于 (W_t) 是当前时刻的白噪声，而 (X_{t-2} - \phi X_{t-1}) 只依赖于 (t-1) 及更早的信息（因果性），因此它们不相关，协方差为零。

所以，对于AR(1)过程，有：
[
\phi_{xx}(1) = \phi
]
[
\phi_{xx}(h) = 0, \quad \text{对于所有 } h > 1
]

PACF的图像特征
AR(1)过程的PACF图会在滞后1处有一个显著的尖峰（值约为 (\phi)），而在所有更高的滞后阶数（(h \ge 2)）上，PACF值理论上为零（在实际样本中会接近零）。这是一个非常清晰的模式。

扩展到AR(p)过程

现在我们将结论推广到一般的AR(p)过程。

考虑一个平稳、因果的AR(p)过程：
[
X_t = \phi_1 X_{t-1} + \phi_2 X_{t-2} + ... + \phi_p X_{t-p} + W_t
]

PACF的行为
对于AR(p)过程，其偏自相关函数具有以下关键特征：

对于滞后阶数 (h \le p)，(\phi_{xx}(h)) 通常非零。
对于所有滞后阶数 (h > p)，有 (\phi_{xx}(h) = 0)。

直观解释
这是因为AR(p)过程具有“p阶马尔可夫性”。要预测 (X_t)，只需要最近p个历史值 (X_{t-1}, ..., X_{t-p})。当我们在计算滞后 (h > p) 的偏自相关时，例如 (X_t) 和 (X_{t-h})，在条件于中间所有的 (h-1) 个观测值后，(X_{t-h}) 中不包含任何预测 (X_t) 的额外信息（这些信息已完全包含在最近的p个值中）。因此，偏自相关为零。

实际意义
这一特性使得PACF成为识别AR过程及其阶数p的强大工具。在PACF图中，我们寻找一个“截尾”点：在该点之后，PACF值变得不显著（在置信区间内波动）。这个截尾点就提示了AR模型的阶数p。

MA(1)过程的PACF

接下来，我们看看移动平均过程的表现。以可逆的MA(1)过程为例：
[
X_t = W_t + \theta W_{t-1}, \quad |\theta| < 1
]

计算PACF
计算MA过程的PACF比AR过程更复杂，因为需要求解回归系数。我们省略繁琐的推导，直接给出MA(1)过程在滞后1处的偏自相关结果：
[
\phi_{xx}(1) = \frac{-\theta^2}{1 + \theta^2 + \theta^4}
]
对于更高的滞后阶数，PACF不会像AR过程那样在某个点后严格为零，而是会逐渐衰减（拖尾），通常其绝对值呈几何衰减趋势，并且可能正负交替。

PACF的行为
对于一般的MA(q)过程，其偏自相关函数 (\phi_{xx}(h)) 不会在有限滞后后截尾，而是逐渐衰减至零。这与AR过程的PACF形成鲜明对比。

总结：ACF与PACF的行为模式

根据以上分析，我们可以总结出AR和MA过程在ACF和PACF上的典型模式，这对于模型识别至关重要。

以下是ACF和PACF在不同模型中的行为对照表：

过程类型	自相关函数 (ACF)	偏自相关函数 (PACF)
AR(p)	逐渐衰减（拖尾）	在滞后p后截尾（即 (h>p) 时，(\phi_{xx}(h) \approx 0)）
MA(q)	在滞后q后截尾（即 (h>q) 时，(\rho(h) \approx 0)）	逐渐衰减（拖尾）
ARMA(p,q)	两者均表现为逐渐衰减（拖尾）	两者均表现为逐渐衰减（拖尾）

关键点

截尾：指函数值在某个滞后点之后突然变得非常接近零（落在置信区间内）。
拖尾：指函数值逐渐衰减至零，可能呈指数型或振荡型衰减，不会突然截断。

R语言实例演示

理论需要结合实际观察。让我们用R语言生成一些模拟数据，直观感受ACF和PACF图。

生成模拟数据

set.seed(479) # 设置随机种子保证结果可重复
n <- 500      # 生成500个观测值

# 1. 生成白噪声
wn <- rnorm(n)

# 2. 生成AR(1)过程，系数phi=0.7
ar1 <- filter(wn, filter=0.7, method="recursive")
ar1 <- ar1[!is.na(ar1)] # 移除初始的NA值

# 3. 生成MA(1)过程，系数theta=0.7
ma1 <- filter(wn, filter=c(1, 0.7), method="convolution", sides=1)
ma1 <- ma1[!is.na(ma1)]

绘制ACF和PACF图

# 设置图形布局，一行两列
par(mfrow=c(1,2))

# 绘制AR(1)的ACF和PACF
acf(ar1, main="ACF of AR(1)")
pacf(ar1, main="PACF of AR(1)")

# 绘制MA(1)的ACF和PACF
acf(ma1, main="ACF of MA(1)")
pacf(ma1, main="PACF of MA(1)")

par(mfrow=c(1,1)) # 恢复默认图形设置

观察结果

AR(1)过程：ACF图呈现缓慢的几何衰减。PACF图在滞后1处有一个明显的尖峰，之后的值基本在零附近波动（截尾）。
MA(1)过程：ACF图在滞后1处有一个尖峰，之后截尾。PACF图则呈现衰减模式，可能正负交替。

更高阶模型与注意事项
当我们模拟AR(3)或MA(3)过程时，模式依然成立，但有时不那么明显，尤其是当系数较小或样本量不足时。对于ARMA混合模型，ACF和PACF都会呈现拖尾现象，仅凭图形精确判断p和q的阶数会更具挑战性。

重要提示
这些图形模式是“经验法则”。在实际数据分析中，需要结合统计检验、信息准则（如AIC）以及模型诊断来综合判断最优模型。样本量越大，这些图形特征通常越清晰。

课程总结

本节课我们一起深入学习了偏自相关函数（PACF）的核心内容。

我们首先回顾了PACF的定义，它衡量了在剔除中间观测值影响后，两个时间点之间的纯相关性。接着，我们深入分析了PACF在AR和MA过程中的理论行为：

对于AR(p)过程，PACF在滞后p之后截尾。这是识别AR模型阶数p的关键特征。
对于MA(q)过程，PACF逐渐衰减（拖尾），而ACF在滞后q之后截尾。

通过R语言实例，我们直观验证了这些理论模式，并认识到在实际应用中，需要谨慎解释图形，并结合其他统计工具进行模型识别。

掌握ACF和PACF的这些特性，是进行时间序列模型识别和定阶的基础。下一讲，我们将开始一个激动人心的新篇章：学习如何从数据中估计这些模型的参数。

011：AR(p)模型的参数估计

在本节课中，我们将学习如何为自回归（AR）过程估计参数。我们将介绍两种主要方法：尤尔-沃克方程和最大似然估计，并比较它们的异同。

概述与目标

到目前为止，我们已经讨论了不同时间序列（如AR、MA和ARMA模型）的性质。现在，我们将进入统计学的核心环节：参数估计。具体来说，对于AR(p)模型，我们的目标是估计其自回归系数 φ₁, φ₂, ..., φₚ 以及白噪声过程的方差 σ²。

我们将从最简单的AR过程开始，逐步介绍两种最标准的估计方法。

第一步：估计均值

首先，我们需要估计时间序列的均值 μ。这很简单，使用样本均值即可：

公式：
x̄ = (1/T) * Σ_{t=1}^{T} x_t

其中，T 是观测数据的总数。计算出样本均值后，我们可以将原始序列减去均值，得到一个均值为0的新序列，从而简化后续计算。

第二步：估计方差与尤尔-沃克方程

上一节我们介绍了如何估计均值，本节中我们来看看如何估计白噪声的方差 σ² 以及自回归系数 φ。

我们的起点是AR(p)模型的定义：

公式：
x_t = w_t + φ₁ x_{t-1} + φ₂ x_{t-2} + ... + φ_p x_{t-p}

其中，w_t 是方差为 σ² 的白噪声。

为了估计参数，我们构建一个方程组。方法是将模型方程两边乘以 x_{t-k}（k=0, 1, ..., p）并取期望。这利用了自协方差函数 γ(k) = Cov(x_t, x_{t-k}) 的性质。

例如，当 k=0 时，我们得到：
γ(0) = σ² + φ₁ γ(1) + φ₂ γ(2) + ... + φ_p γ(p)

当 k=1, 2, ..., p 时，我们得到一系列方程：
γ(1) = φ₁ γ(0) + φ₂ γ(1) + ... + φ_p γ(p-1)
γ(2) = φ₁ γ(1) + φ₂ γ(0) + ... + φ_p γ(p-2)
...
γ(p) = φ₁ γ(p-1) + φ₂ γ(p-2) + ... + φ_p γ(0)

这个方程组被称为 尤尔-沃克方程。

我们可以用矩阵形式简洁地表示它：

公式：
Γ φ = γ

其中：

φ = [φ₁, φ₂, ..., φ_p]^T 是待估参数向量。
γ = [γ(1), γ(2), ..., γ(p)]^T 是自协方差向量（从滞后1到p）。
Γ 是一个 p×p 矩阵，其第 i 行第 j 列元素为 γ(|i-j|)。

在实际操作中，我们用样本自协方差 γ̂(k) 替换理论值 γ(k)，得到估计矩阵 Γ̂ 和向量 γ̂。然后，通过求解线性方程组来估计 φ：

公式：
φ̂_{YW} = Γ̂^{-1} γ̂

得到 φ 的估计值后，我们可以回代到第一个方程中，估计白噪声的方差：

公式：
σ̂²_{YW} = γ̂(0) - Σ_{i=1}^{p} φ̂_i γ̂(i)

最大似然估计方法

上一节我们介绍了基于矩估计的尤尔-沃克方法，本节中我们来看看基于概率模型的最大似然估计方法。

MLE要求我们指定误差项的分布。为了使推导可行，我们通常假设白噪声 w_t 服从高斯（正态）分布。我们以AR(1)模型为例进行说明，因为其形式更简单，但原理可以推广到AR(p)模型。

AR(1)模型为：
x_t - μ = φ (x_{t-1} - μ) + w_t, 其中 w_t ~ N(0, σ²)

我们的目标是基于观测数据 x₁, x₂, ..., x_T，找到参数 μ, φ, σ² 的值，使得数据的联合概率密度（即似然函数）最大。

由于AR(1)过程具有马尔可夫性（当前值只依赖于前一个值），我们可以将联合似然函数分解：

公式：
L(μ, φ, σ²) = f(x₁) * Π_{t=2}^{T} f(x_t | x_{t-1})

其中：

f(x₁) 是第一个观测值的边际分布，对于平稳AR(1)过程，x₁ ~ N(μ, σ²/(1-φ²))。
f(x_t | x_{t-1}) 是条件分布，x_t | x_{t-1} ~ N(μ + φ(x_{t-1} - μ), σ²)。

将具体的正态分布密度函数代入并取对数，我们得到对数似然函数。然而，由于 f(x₁) 项中包含 1/(1-φ²)，使得整个方程关于参数是非线性的，求解解析解非常困难。

条件最大似然估计

为了解决MLE中的非线性问题，一个常用的技巧是使用条件最大似然估计。其思想是：我们“忽略”第一个观测值 x₁ 的分布，只基于 x₁ 已知的条件下来估计后续观测值的似然。

条件似然函数为：
L_c(μ, φ, σ² | x₁) = Π_{t=2}^{T} f(x_t | x_{t-1})

代入条件正态分布，并取对数后，我们需要最大化的核心部分是一个条件平方和：

公式（对于AR(1)）：
S_c(μ, φ) = Σ_{t=2}^{T} [ (x_t - μ) - φ (x_{t-1} - μ) ]²

惊喜出现了！如果我们令 α = μ(1-φ)，那么上述条件平方和可以重写为：

公式：
S_c(α, φ) = Σ_{t=2}^{T} [ x_t - α - φ x_{t-1} ]²

这恰好是一个简单线性回归模型的残差平方和！其中，响应变量是 x_t，预测变量是 x_{t-1}，截距是 α，斜率是 φ。

因此，我们可以直接应用最小二乘法的公式来估计 α 和 φ：

公式：
[α̂, φ̂]^T = (X^T X)^{-1} X^T Y

其中，设计矩阵 X 的第一列全为1（对应截距），第二列为 [x₁, x₂, ..., x_{T-1}]^T；响应向量 Y = [x₂, x₃, ..., x_T]^T。

得到 φ̂ 和 α̂ 后，我们可以反解出均值估计 μ̂ = α̂ / (1 - φ̂)。最后，方差的CMLE估计为：

公式：
σ̂²_{CMLE} = S_c(μ̂, φ̂) / (T-1)

尤尔-沃克估计与最大似然估计的比较

现在，让我们对两种方法进行总结和比较。

对于AR(1)模型，尤尔-沃克估计量 φ̂_{YW} 本质上是样本自相关系数：

公式：
φ̂_{YW} = γ̂(1) / γ̂(0) = [ Σ (x_t - x̄)(x_{t-1} - x̄) ] / [ Σ (x_t - x̄)² ]
（求和范围通常为 t=2 到 T，但均值 x̄ 使用全部 T 个数据计算）

而条件最小二乘估计量 φ̂_{CLS} 形式上非常相似，但在计算均值时略有不同：

公式：
φ̂_{CLS} = [ Σ (x_t - x̄₂)(x_{t-1} - x̄₁) ] / [ Σ (x_{t-1} - x̄₁)² ]
其中 x̄₁ 是 x₁ 到 x_{T-1} 的均值，x̄₂ 是 x₂ 到 x_T 的均值。

核心区别在于如何处理数据的边界（首尾数据）。

尤尔-沃克方程：基于样本自协方差，隐含着对序列平稳性的假设，处理边界时有一种对称性。
条件MLE/最小二乘：明确地将第一个观测值作为条件，从第二个点开始进行回归，更直接地处理了有限数据集的起始问题。

结论：

当样本量 T 很大时，两种估计方法给出的结果非常接近，都是相合估计量。
对于小样本，它们可能会给出略有差异的估计。
尤尔-沃克方法总能保证估计出的模型是平稳的（即特征根在单位圆内），这是一个优点。
条件最小二乘/最大似然方法在计算上更直观，直接与回归框架挂钩，并且更容易推广到更复杂的模型（如ARMA）。

在实际应用中（例如R语言中的 ar() 函数），通常会提供多种方法供选择，理解其背后的原理有助于我们做出合适的选择。

总结

在本节课中，我们一起学习了为AR(p)过程估计参数的核心技术。

我们首先介绍了尤尔-沃克方程，它通过建立自协方差与自回归系数之间的线性方程组来求解参数，是一种基于矩估计的方法。

接着，我们探讨了最大似然估计。为了规避解析求解的复杂性，我们引入了条件最大似然估计，它将问题巧妙地转化为一个最小二乘回归问题，从而可以方便地求解。

最后，我们比较了这两种主流方法。它们本质上是相通的，主要区别在于对数据边界条件的处理方式。对于大样本，两者结果一致；对于小样本，则各有特点。

参数估计是时间序列建模的基础。在接下来的课程中，我们将利用估计好的模型，进行另一个至关重要的任务：预测未来值。

012：AR(p)模型的预测 🎯

在本节课中，我们将要学习如何对自回归过程进行预测。我们已经讨论了如何估计AR过程的参数，但通常我们不仅想知道参数是什么，还想知道过程未来的走向，这就是预测。

从估计到预测 🔄

上一节我们介绍了如何估计AR过程的系数，使用了尤尔-沃克方程和最大似然估计等方法。本节中，我们来看看如何利用这些估计结果来预测未来的值。

预测的核心问题是：给定从时间1到时间T的观测数据 (X_1, X_2, ..., X_T)，我们想要预测未来的值 (X_{T+1}, X_{T+2}, ..., X_{T+M})。

预测的基本思路 📈

我们假设时间序列是平稳的，但不一定是AR过程。目标是使用过去的数据来线性预测未来值。我们采用最小二乘法，这是统计学中常用的方法。

我们的预测模型形式如下：
[
\hat{X}{T+M} = a_0 + \sum^{T} a_t X_t
]
其中，(a_0, a_1, ..., a_T) 是需要估计的系数。

为了找到这些系数，我们求解以下最小二乘问题：
[
\min_{a_0, ..., a_T} E\left[ (X_{T+M} - \hat{X}_{T+M})^2 \right]
]

求解预测系数 🧮

通过对系数求导并令导数为零，我们得到一系列方程。首先，关于 (a_0) 求导得到：
[
E[X_{T+M} - \hat{X}{T+M}] = 0
]
这意味着预测值的期望必须等于过程本身的期望 (\mu)。由此可以推导出 (a_0 = \mu(1 - \sum^{T} a_t))。这表明，我们可以通过将过程中心化（减去均值 (\bar{X})）来简化问题，从而假设过程均值为零继续分析。

接下来，关于其他系数 (a_t) 求导，我们得到形如下式的方程：
[
E[(X_{T+M} - \hat{X}_{T+M}) X_t] = 0, \quad t = 1, 2, ..., T
]

将这些方程展开，会涉及到序列的自协方差。最终，我们得到一个庞大的线性方程组：
[
\mathbf{k} = \mathbf{\Gamma} \mathbf{a}
]
其中：

(\mathbf{k}) 是一个向量，其元素为 (Cov(X_{T+M}, X_t))，对于一步预测（M=1），这通常是自协方差函数在滞后 (T, T-1, ..., 1) 处的值。
(\mathbf{\Gamma}) 是一个 (T \times T) 的矩阵，其第 (i, j) 个元素为 (Cov(X_i, X_j) = \gamma_X(i-j))，即托普利茨矩阵。
(\mathbf{a}) 是系数向量 ([a_1, a_2, ..., a_T]^T)。

理论上，系数解为 (\mathbf{a} = \mathbf{\Gamma}^{-1} \mathbf{k})，一步预测值为 (\hat{X}_{T+1} = \mathbf{k}^T \mathbf{\Gamma}^{-1} \mathbf{x})，其中 (\mathbf{x}) 是观测数据向量。

预测的均方误差 📉

我们可以计算这种线性预测器的均方误差（MSE）：
[
MSE = E[(X_{T+1} - \hat{X}_{T+1})^2] = \gamma_X(0) - \mathbf{k}^T \mathbf{\Gamma}^{-1} \mathbf{k}
]
这个结果有很好的解释：

(\gamma_X(0)) 是过程的方差，代表没有任何信息时的预测误差。
(\mathbf{k}^T \mathbf{\Gamma}^{-1} \mathbf{k}) 是一个非负项，代表了利用历史数据 (X_1, ..., X_T) 所获得的信息，它减少了预测误差。
由于 (\mathbf{\Gamma}) 是半正定的，所以信息项非负，确保 MSE 不超过原始方差。

实践中的挑战与简化 💡

上述理论框架虽然优美，但在实践中面临巨大挑战：矩阵 (\mathbf{\Gamma}) 是 (T \times T) 的，对于长时间序列（如T=500或更大），求逆计算量巨大且数值不稳定。此外，我们还需要估计所有滞后阶数的自协方差，远距离滞后的估计本身就不准确。

然而，如果我们事先知道数据来自一个 (AR(p)) 过程，那么问题将大大简化。对于 (AR(p)) 过程，最优的一步预测器只依赖于最近的 (p) 个观测值：
[
\hat{X}{T+1} = \phi_1 X_T + \phi_2 X + ... + \phi_p X_{T-p+1}
]
其中系数 (\phi_i) 可以通过尤尔-沃克方程（一个仅 (p \times p) 的线性系统）估计得到。这意味着我们无需处理庞大的 (T \times T) 矩阵。

如果不知道过程是否是AR的，我们则需要更聪明的方法来求解大型托普利茨系统，这涉及到数值线性代数中的专门算法，我们将在后续课程中探讨。

R语言中的预测实践 🖥️

在R语言中，我们可以使用 stats 包中的 ar() 函数来拟合AR模型并进行预测。该函数默认使用尤尔-沃克方法，但也提供了最大似然估计（MLE）、普通最小二乘法（OLS）和伯格（Burg）算法等选项。通常，这些方法对于纯AR过程会给出相似的结果。

以下是使用 ar() 函数和 predict() 函数进行拟合和预测的基本代码示例：

# 设置随机种子以保证结果可重现
set.seed(479)

# 模拟一个AR(1)过程，系数为0.7
ar1_data <- arima.sim(model = list(ar = 0.7), n = 100)

# 使用ar函数拟合模型（默认YW方法）
fit_model <- ar(ar1_data)

# 查看拟合的系数和模型阶数
print(fit_model)

# 进行未来5步的预测
predictions <- predict(fit_model, n.ahead = 5)
print(predictions$pred)  # 预测值
print(predictions$se)    # 预测标准误

需要注意的是，ar() 函数会尝试用AR模型去拟合任何时间序列。如果数据本身是移动平均（MA）过程或ARMA过程，ar() 可能会拟合出一个高阶的AR模型，其系数会呈现几何衰减，这近似对应了可逆MA过程的无限阶AR表示。

总结 📚

本节课中我们一起学习了时间序列预测的核心概念。

我们从一般化的线性预测框架出发，推导了基于全部历史数据的最小二乘预测公式，它涉及求解一个以自协方差矩阵为核心的大型线性系统。
我们分析了预测的均方误差，并理解了历史信息如何降低预测的不确定性。
我们指出了直接求解该大型系统在计算上的不可行性。
作为一种重要的特例，我们了解到对于 (AR(p)) 过程，预测可以简化为仅使用最近 (p) 期数据，这通过尤尔-沃克方程等易于处理的方法实现。
最后，我们简要介绍了在R语言中执行AR模型拟合和预测的基本操作。

下一讲，我们将探讨如何高效地解决大型托普利茨系统问题，以实现更一般化场景下的预测。

013：Durbin-Levinson与创新算法

在本节课中，我们将进一步学习自回归过程的预测方法。上节课我们介绍了线性预测，但求解相关方程组较为复杂。本节课我们将介绍两种解决该问题的算法：Durbin-Levinson算法和创新算法。

上节课我们介绍了如何基于前T个数据点预测第T+1个点，这最终归结为求解一个大型线性方程组。直接求解这个系统（例如通过矩阵求逆）在计算上代价高昂。幸运的是，我们可以利用自协方差矩阵的特殊结构，发展出更高效的算法。

Durbin-Levinson算法

上一节我们介绍了线性预测的基本方程。本节中我们来看看Durbin-Levinson算法，它利用自协方差矩阵的托普利茨结构，以迭代方式高效求解预测系数。

我们想要预测X_{T+1}，其线性预测器形式为：

X̂_{T+1} = α_{1}X_T + α_{2}X_{T-1} + ... + α_{T}X_1

其中系数向量α需要求解方程 Γ_T α = k。这里Γ_T是T×T的自协方差矩阵，k是自协方差向量[γ_X(T), γ_X(T-1), ..., γ_X(1)]^T。直接求逆Γ_T需要O(T³)次运算。

Durbin-Levinson算法的核心思想是迭代求解。我们不是一次性求解所有系数，而是逐步构建从1步到T步的预测器。

以下是该算法的迭代步骤：

初始化：
- 设 α_{1,1} = 0
- 设 P_{1,0} = γ_X(0) （即用零信息预测X_1的均方预测误差，就是方差）

迭代更新 (对于 n = 2, 3, ..., T)：

计算第n个预测器的第一个系数：

α_{n, n} = (γ_X(n) - Σ_{i=1}^{n-1} α_{n-1, i} γ_X(n-i)) / P_{n-1, n}

更新其他系数 (对于 i = 1, 2, ..., n-1)：

α_{n, i} = α_{n-1, i} - α_{n, n} α_{n-1, n-i}

更新均方预测误差：

P_{n, n+1} = P_{n-1, n} (1 - α_{n, n}²)

该算法将计算复杂度从O(T³)降低到了O(T²)。对于自回归过程，由于只需要回顾有限阶数（p阶），计算可以进一步简化。

创新算法

上一节我们介绍了基于原始观测值的Durbin-Levinson算法。本节中我们来看看另一种思路：基于预测残差（称为“创新”）进行预测的创新算法。

一个“创新”是指一步超前预测的残差：

ν_t = X_t - X̂_{t | t-1}

其中X̂_{t | t-1}是基于直到时间t-1的所有观测值对X_t的预测。关键性质是，不同时间的创新ν_t和ν_s（当t ≠ s时）是不相关的。因此，创新向量的协方差矩阵是对角阵，这大大简化了计算。

创新算法的目标是使用过去的创新来预测未来值，其预测形式为：

X̂_{T+1 | T} = Σ_{j=1}^{T} θ_{T, j} ν_{T+1-j}

以下是该算法的迭代步骤：

初始化：
- 设 X̂_{1 | 0} = 0 （假设过程均值为零）
- 设 P_{1,0} = γ_X(0)

迭代更新 (对于 T = 1, 2, ...)：

计算系数θ_{T, T-j} (对于 j = 0, 1, ..., T-1)：

θ_{T, T-j} = (γ_X(T-j) - Σ_{k=0}^{j-1} θ_{T, T-k} θ_{j, j-k} P_{k+1, k}) / P_{j+1, j}

更新预测：

X̂_{T+1 | T} = Σ_{j=0}^{T-1} θ_{T, T-j} ν_{T+1-j}

更新均方预测误差：

P_{T+1, T} = γ_X(0) - Σ_{j=0}^{T-1} θ_{T, T-j}² P_{j+1, j}

虽然公式看起来复杂，但其核心思想是用不相关的创新序列替代相关的原始观测序列来构建预测器。

算法比较与应用场景

我们已经介绍了两种预测算法。本节中我们来看看它们各自适用的模型类型。

Durbin-Levinson算法更适用于AR过程：对于AR(p)过程，预测X_{T+1}实际上只需要最近的p个观测值。Durbin-Levinson算法能有效利用这一特性。
创新算法更适用于MA和ARMA过程：对于可逆的MA(q)过程，它可以表示为过去所有观测值的无限加权和。创新算法通过使用不相关的残差，避免了处理长程相关的困难，因此对MA和ARMA模型特别有效。

为了更直观地理解，我们考虑一个MA(1)过程的简单例子：X_t = θ W_{t-1} + W_t，其中W_t是白噪声。

应用创新算法后，我们会发现大多数系数θ_{T, j}（当j > 1时）为零。预测方程简化为：

X̂_{T+1 | T} ≈ (某个标度因子) * θ * ν_T

这直观地说明，对于MA(1)过程，下一步的预测主要取决于当前的一步预测误差（创新）。我们试图用当前的残差来估计“隐藏”的白噪声W_t，进而进行预测。相比之下，AR(1)过程X_t = φ X_{t-1} + W_t的预测则直接基于上一个观测值：X̂_{T+1 | T} ≈ φ X_T。

本节课中我们一起学习了两种重要的时间序列预测算法：Durbin-Levinson算法和创新算法。Durbin-Levinson算法通过迭代方式高效求解托普利茨方程组，适用于AR类过程。创新算法则通过将预测建立在互不相关的创新序列之上，特别适用于MA和ARMA类过程。理解这些算法的基本思想，有助于我们在使用统计软件进行预测时，了解其背后的计算逻辑。下节课我们将进入R语言实践，学习如何应用这些方法进行实际的时间序列预测。

014：ARMA模型的估计与预测

在本节课中，我们将学习如何对ARMA模型进行参数估计和预测。我们将沿用处理AR模型时使用的工具，特别是最大似然估计法，但会引入移动平均（MA）部分带来的新挑战。课程将涵盖估计方法、估计量的渐近性质、模型阶数选择以及预测的基本原理。

估计ARMA(p, q)模型

上一节我们介绍了AR模型的估计方法，本节中我们来看看如何将其扩展到ARMA模型。ARMA模型包含自回归（AR）和移动平均（MA）两部分，因此我们需要估计的参数更多，包括均值μ、白噪声方差σ²、自回归参数φ₁到φ_p以及移动平均参数θ₁到θ_q。

最大似然估计法

我们将继续使用最大似然估计法。这需要一个关键假设：驱动时间序列的白噪声过程是高斯（正态）的。如果白噪声不是高斯的，那么ARMA过程本身也不是高斯的，我们就不能直接使用最大似然估计。

与AR模型类似，我们不直接写出整个时间序列的联合密度函数，而是将其分解为一系列条件密度的乘积：
L(μ, σ², φ, θ) = ∏_{t=1}^{T} f(x_t | x_{t-1}, ..., x_1)
其中，每个条件密度f都是高斯分布。

对于ARMA模型，给定过去所有观测值x_{t-1}, ..., x_1，x_t的条件分布是一个正态分布：

均值：一步超前预测值 \hat{x}_{t-1}^{t}
方差：一步超前预测的均方误差 P_{t}^{t-1}

因此，似然函数可以写为：
L ∝ ∏_{t=1}^{T} (P_{t}^{t-1})^{-1/2} exp{ - (x_t - \hat{x}_{t-1}^{t})² / (2 P_{t}^{t-1}) }

处理预测误差方差

我们需要计算每个时间点的预测误差方差P_{t}^{t-1}。对于t=1，P_{1}^{0}就是过程的平稳方差，可以表示为σ²乘以一个由ψ系数（来自模型的因果线性表示）决定的常数。

对于t>1，我们可以使用Durbin-Levinson算法进行递归更新。关键点在于，所有的P_{t}^{t-1}都可以写成σ² * R_t的形式，其中R_t是一个不依赖于σ²的量。

将这个关系代入似然函数，我们可以将其重写，从而将σ²与其他参数分离开。通过标准的微积分运算（对σ²求导并令导数为零），我们可以得到方差σ²的最大似然估计量：
\hat{σ}² = S(μ, φ, θ) / T
其中S(μ, φ, θ)是一个与数据和参数有关的平方和项。

集中似然函数

得到σ²的估计量后，我们可以将其代回似然函数，从而得到一个集中似然函数，它现在只依赖于参数μ, φ, θ：
L_c(μ, φ, θ) ∝ [S(μ, φ, θ)/T]^{-T/2} * (∏_{t=1}^{T} R_t)^{-1/2}
然后，我们通过最大化这个集中似然函数（或其对数形式）来求解μ, φ, θ的估计值。这个过程通常需要数值优化方法来完成。

以下是估计ARMA模型参数的核心步骤总结：

设定模型阶数：指定自回归阶数p和移动平均阶数q。
构建似然函数：基于高斯白噪声假设，写出条件似然函数。
分离方差参数：将预测误差方差表示为σ² * R_t，简化似然函数。
求解方差估计：首先得到白噪声方差σ²的估计量\hat{σ}²。
优化其他参数：将\hat{σ}²代回，得到集中似然函数，并优化求解μ, φ, θ的估计值。

在实际应用中，我们通常使用R语言中的arima()等函数来执行这些复杂的计算。

估计量的渐近性质

在统计学中，我们估计的参数通常具有渐近正态性，ARMA模型的估计量也不例外。这意味着当样本量T足够大时，估计量的分布近似于正态分布。这为我们进行假设检验（例如，检验某个参数是否显著不为零）和构建置信区间提供了理论基础。

具体来说，令β表示所有参数（φ和θ）组成的向量，\hat{β}是其最大似然估计量。那么，当T → ∞时，有：
√T (\hat{β} - β) → N(0, σ² Γ_{p,q}^{-1})
其中，Γ_{p,q}是一个分块矩阵，其元素由两个辅助过程y_t和y'_t的自协方差和互协方差函数决定：

y_t：原ARMA过程x_t中仅保留AR部分所构成的AR过程，即φ(B)y_t = w_t。
y'_t：原ARMA过程x_t中仅保留MA部分所构成的MA过程，即y'_t = θ(B)w_t。

示例：AR(1)和MA(1)

对于AR(1)模型x_t = φ x_{t-1} + w_t，其参数φ的估计量满足：
√T (\hat{φ} - φ) → N(0, 1 - φ²)
有趣的是，当φ接近±1时（但仍保持平稳），估计量的渐近方差会变小。当φ=0时，过程退化为白噪声，估计量的方差就是白噪声的方差。

对于可逆的MA(1)模型x_t = w_t + θ w_{t-1}，通过将其视为一个AR过程，我们可以得到完全对称的结果：
√T (\hat{θ} - θ) → N(0, 1 - θ²)

如何选择模型阶数p和q？

在拟合ARMA模型之前，我们需要确定自回归阶数p和移动平均阶数q。常用的方法是使用模型选择准则，例如AIC（赤池信息准则）或BIC（贝叶斯信息准则）。

以下是选择模型阶数的常用方法：

自动拟合函数：R语言forecast包中的auto.arima()函数可以自动搜索一系列(p, q)组合，并选择AIC或BIC值最小的模型。
手动比较：使用stats包中的arima()函数手动拟合多个不同(p, q)组合的模型，然后比较它们的AIC或BIC值，选择最小的一个。

通常AIC和BIC会给出相同的结果，但并非总是如此。

ARMA模型的预测

现在，我们转向ARMA模型的预测问题。我们的目标是基于到时间T为止的观测数据x_1, ..., x_T，预测未来h步的值x_{T+h}。最优预测（在均方误差最小意义下）是条件期望：
\hat{x}_{T+h} = E[ x_{T+h} | x_T, x_{T-1}, ..., x_1 ]

基于无限过去数据的预测（理论构建）

为了数学上的便利，我们首先考虑一个理论情形：假设我们可以观测到无限过去的数据..., x_0, x_1, ..., x_T。基于此的预测我们记为\tilde{x}_{T+h}。当实际数据量T很大时，\tilde{x}_{T+h}和\hat{x}_{T+h}会非常接近。

我们利用ARMA模型的两种表示：

因果表示（无限MA）：x_{T+h} = ∑_{j=0}^{∞} ψ_j w_{T+h-j}
可逆表示（无限AR）：w_{T+h} = ∑_{j=0}^{∞} π_j x_{T+h-j} （其中π_0 = 1）

对因果表示两边取条件期望，由于未来白噪声w_{T+1}, ..., w_{T+h}与过去数据独立，其条件期望为0。因此：
\tilde{x}_{T+h} = ∑_{j=h}^{∞} ψ_j w_{T+h-j}
预测误差x_{T+h} - \tilde{x}_{T+h} = ∑_{j=0}^{h-1} ψ_j w_{T+h-j}，对应的预测均方误差为：
P_{T+h}^{T} = σ² ∑_{j=0}^{h-1} ψ_j²

对可逆表示两边取条件期望，我们可以得到\tilde{x}_{T+h}的一个递归表达式，它依赖于已知数据x_T, x_{T-1}, ...以及之前步的预测值\tilde{x}_{T+1}, ..., \tilde{x}_{T+h-1}。这说明多步预测可以递归进行：先预测一步，然后用这个预测值去预测两步，依此类推。

长期预测行为

当预测步数h非常大时，会发生什么？由于因果表示中的系数ψ_j是绝对可和的，其尾部∑_{j=h}^{∞} |ψ_j|会随着h增大而趋于0。这意味着：

预测值趋近于均值：\tilde{x}_{T+h} → μ （当h → ∞）
预测误差趋近于过程方差：P_{T+h}^{T} → σ² ∑_{j=0}^{∞} ψ_j² = Var(x_t) （当h → ∞）

直观解释是：对于一个平稳过程，随着我们预测得越来越远，过去数据的影响逐渐消失，我们能做的最好预测就是过程的长期均值μ，而预测的不确定性则趋近于过程本身的波动性（方差）。

在实际中，我们只有有限数据。下一讲我们将探讨如何将理论预测\tilde{x}_{T+h}转化为基于有限数据x_1, ..., x_T的实用预测\hat{x}_{T+h}。

总结

本节课中我们一起学习了ARMA模型的估计与预测基础：

参数估计：我们使用最大似然估计法，需要假设高斯白噪声。通过将联合似然分解为条件似然、分离方差参数并构建集中似然函数，可以估计所有参数（μ, σ², φ, θ）。
渐近性质：在大样本下，参数估计量服从渐近正态分布，这为统计推断（如置信区间、假设检验）提供了依据。
模型选择：可以使用AIC或BIC等准则来选择ARMA模型的合适阶数(p, q)。
预测原理：最优预测是条件期望。基于无限过去数据的理论预测\tilde{x}_{T+h}可以通过模型的因果或可逆表示来计算，其预测误差方差可以明确写出。长期预测会收敛到过程均值。
实际应用：尽管推导复杂，但R语言中的现有函数（如arima(), auto.arima()）可以方便地完成估计、模型选择和预测工作。理解其背后的原理有助于我们正确使用和解释这些工具的结果。

015：ARMA模型预测 📈

在本节课中，我们将要学习如何对ARMA模型进行预测。上一讲我们介绍了ARMA过程，讨论了如何估计系数，并开始探讨预测的数学原理。本节中，我们将深入探讨如何利用真实数据进行预测。

回顾与数学基础

首先，我们回顾上一讲引入的符号和概念。对于一个ARMA(p, q)过程，其模型为：
Φ(B)X_t = Θ(B)W_t
我们假设该过程是因果且可逆的。

在预测中，我们有两个核心预测量：

X̂_{T+h}：基于观测数据 X_T, X_{T-1}, ..., X_1 对 X_{T+h} 的预测。
X̃_{T+h}：基于观测数据及假设延伸到无限过去的“完整”数据对 X_{T+h} 的预测。

上一讲我们主要讨论了 X̃_{T+h}。我们将其写成因果形式：
X̃_{T+h} = μ + Σ_{j=h}^{∞} ψ_j W_{T+h-j}

由此，我们得到了两个重要性质：

当预测步长 h → ∞ 时，预测值 X̃_{T+h} 会收敛于过程的均值 μ。
预测的均方误差 P_{T+h|T} 为：
P_{T+h|T} = σ_w² Σ_{j=0}^{h-1} ψ_j²
随着 h 增大，它将收敛于过程的方差。

为了更严谨地证明第一个性质（收敛于均值），我们可以使用切比雪夫不等式。我们想证明对于任意 ε > 0，当 h → ∞ 时，概率 P(|X̃_{T+h} - μ| > ε) 趋于 0。利用切比雪夫不等式，该概率的上界为 Var(X̃_{T+h} - μ) / ε²，即 (σ_w² / ε²) Σ_{j=h}^{∞} ψ_j²。由于因果表示要求 Σ_{j=0}^{∞} ψ_j² < ∞，因此其尾部求和 Σ_{j=h}^{∞} ψ_j² 必然随着 h 增大而趋于 0，从而证明了依概率收敛。

处理有限数据：截断无限过去

然而，现实问题是我们只有有限的数据（X_1 到 X_T），而非无限的历史数据。因此，我们需要处理如何用有限数据计算 X̂_{T+h}。

主要有两种方法：

如果数据量 T 较小，我们可以直接求解一个线性方程组。
如果 T 较大，我们需要迭代地进行计算。

我们将利用ARMA模型的可逆表示进行迭代预测。可逆表示将白噪声 W_t 表示为 X_t 的无限线性组合。对模型两边取条件期望，我们可以得到预测的递推公式。

具体而言，对于已知系数 φ₁...φ_p 和 θ₁...θ_q 的ARMA模型，h步预测可以写成：
X̂_{T+h} = Σ_{j=1}^{p} φ_j X̂_{T+h-j} + Σ_{j=1}^{q} θ_j Ŵ_{T+h-j}
其中，对于 j < h，Ŵ_{T+h-j} = 0（未来白噪声的期望为0）；对于 j ≥ h，Ŵ_{T+h-j} = W_{T+h-j}（过去白噪声已知或可估计）。

而白噪声项 Ŵ_t 本身也可以通过一个类似的递推公式，利用观测数据 X_t 和之前的白噪声估计值来计算。

因此，预测是一个逐步进行的过程：先利用数据估计出必要的白噪声序列，然后从 h=1 开始，依次计算 h=2, 3, ... 的预测值，每一步都依赖于前一步的预测结果。

ARMA(1,1) 预测示例

让我们通过一个ARMA(1,1)模型的例子来具体说明这个过程。模型为：
X_{t+1} = φ X_t + W_{t+1} + θ W_t

一步预测 (h=1)：
X̂_{T+1} = φ X_T + θ Ŵ_T
我们需要知道 Ŵ_T。
多步预测 (h ≥ 2)：
X̂_{T+h} = φ X̂_{T+h-1}
因为对于 h≥2，未来白噪声项的期望为0。

那么如何得到 Ŵ_T 呢？我们可以从模型方程迭代求解。从 t=0 开始，假设 Ŵ_0 = 0，然后利用公式 Ŵ_t = X_t - φ X_{t-1} - θ Ŵ_{t-1} 依次计算出 Ŵ_1, Ŵ_2, ..., Ŵ_T。

对于这个ARMA(1,1)模型，我们也可以计算其预测均方误差。首先写出其因果表示：X_t = W_t + Σ_{i=1}^{∞} ψ_i W_{t-i}，其中 ψ_i = (φ + θ) φ^{i-1}。那么 h 步预测的均方误差为：
P_{T+h|T} = σ_w² [1 + (φ+θ)² Σ_{j=1}^{h-1} φ^{2(j-1)}]
这是一个有限几何级数，可以求和。当 h → ∞ 时，它收敛于过程的方差 σ_w² [1 + (φ+θ)²/(1-φ²)]，这与我们之前对ARMA(1,1)方差的分析一致。

反向预测

除了预测未来，有时我们也可能希望预测过去，这被称为“反向预测”。其目标是根据数据 X_1, ..., X_T 来预测 X_{1-h}。

其数学原理与正向预测完全对称。我们寻找一个线性估计量 X̂_{1-h} = Σ_{i=1}^{T} α_i X_i。通过最小化均方误差，我们同样会得到一个由自协方差函数构成的线性方程组 k = Γα，其中 k 是包含 γ(h), γ(h+1), ... 的向量，Γ 是自协方差矩阵。求解这个方程组即可得到系数 α，进而得到预测值。

一个有趣的事实是：对于一个平稳高斯过程，向量 (X_{T+1}, ..., X_1) 的联合分布与向量 (X_0, ..., X_{T}) 的联合分布相同。这意味着，对于平稳高斯过程，时间正向运行与反向运行在概率意义上是不可区分的，因此正向预测和反向预测在本质上是相同的问题。

R语言实战：估计与预测

理论之后，让我们通过R语言来看看如何对真实数据和模拟数据进行ARMA模型的估计与预测。

首先，我们模拟一些数据进行分析。

set.seed(2479) # 设置随机种子保证结果可重复
# 生成AR(2)过程
ar2_data <- arima.sim(model = list(ar = c(0.7, 0.2)), n = 500)
# 生成MA(2)过程
ma2_data <- arima.sim(model = list(ma = c(0.7, 0.2)), n = 500)
# 生成ARMA(2,2)过程
arma22_data <- arima.sim(model = list(ar = c(0.7, 0.2), ma = c(0.7, 0.2)), n = 500)

模型估计

我们可以使用Arima函数（来自forecast包）或arma函数（来自TSA包）来拟合模型。选择模型时，可以依赖AIC（赤池信息准则）或BIC（贝叶斯信息准则）等准则，数值较小的模型通常更优。

library(forecast)
# 为AR(2)数据拟合正确的AR(2)模型
m_correct <- Arima(ar2_data, order = c(2, 0, 0))
summary(m_correct)
# 尝试拟合错误的ARMA(2,2)模型
m_wrong <- Arima(ar2_data, order = c(2, 0, 2))
summary(m_wrong)

比较两个模型的AIC，通常正确模型的AIC会更小。我们也可以使用auto.arima函数让R自动寻找最优的ARIMA模型。

auto_model <- auto.arima(ar2_data)
summary(auto_model)

auto.arima会尝试多种模型组合（包括差分阶数），并返回AIC最小的模型。这对于初步探索非常有用。

进行预测

使用forecast函数可以对拟合好的模型进行预测。

# 对ARMA(2,2)模拟数据进行未来20步的预测
fc <- forecast(arma22_data, h = 20)
plot(fc, main = "ARMA(2,2)过程预测",
     xlim = c(400, 520), ylim = c(-5, 5))

预测图会显示历史数据、点预测值以及80%和95%的预测区间。可以看到，随着预测步长增加，预测区间会迅速变宽，点预测最终会趋向于序列的均值。

评估预测精度

我们可以使用accuracy函数来评估模型的预测精度。它返回一系列误差指标，如平均绝对误差（MAE）、均方根误差（RMSE）等。这有助于在不同模型之间进行比较。

# 假设我们有两个候选模型：m1和m2
accuracy(m1) # 训练集精度
# 比较两个模型在训练集上的精度
accuracy(m1, arma22_data)
accuracy(m2, arma22_data)

真实数据示例：油价数据

最后，我们使用forecast包内置的oil.price数据集进行实战。

data(oil.price)
oil <- window(oil.price, end = 1990) # 截取1990年之前的数据以避开后期巨大波动
plot(oil, main = "油价数据（截取）")
# 自动寻找合适模型
oil_model <- auto.arima(oil)
summary(oil_model)
# 进行预测并绘图
fc_oil <- forecast(oil_model, h = 24) # 预测未来两年
plot(fc_oil, main = "油价预测")

通过这个过程，我们可以看到时间序列分析是一个探索性很强的领域。我们需要尝试不同的模型，结合AIC/BIC、系数显著性、预测精度等多种信息，来选择一个相对合理且有用的模型。没有绝对正确的答案，关键在于理解数据并找到能合理解释和预测其行为的模型。

总结

本节课中我们一起学习了ARMA模型的预测。

我们首先回顾了基于无限历史数据的理想预测公式及其性质（收敛于均值）。
接着，我们重点讨论了如何利用有限数据进行实际预测，核心方法是利用模型的可逆表示进行迭代计算。
我们还简要介绍了反向预测的概念及其在平稳高斯过程中的对称性。
最后，通过R语言实战，我们演示了如何对模拟数据和真实油价数据进行模型估计、预测以及精度评估的完整流程。

下一讲，我们将探讨带有季节性成分的时间序列模型。

016：季节性ARIMA模型

在本节课中，我们将学习季节性ARIMA模型。这些模型用于分析具有周期性或季节性趋势的时间序列数据，例如月度销售额或年度温度变化。我们将探讨如何将之前学到的ARIMA模型扩展到包含季节性成分，并理解其数学表示和性质。

引言

上一节我们介绍了标准的ARIMA模型。本节中，我们将探讨如何将这些模型扩展以捕捉数据中的季节性模式。许多现实世界的数据，如零售销售额或气温，会表现出以固定周期（如月、季、年）重复的模式。季节性ARIMA模型正是为此类数据设计的。

季节性ARIMA模型的基本概念

季节性ARIMA模型的核心思想是，在标准ARIMA模型的基础上，额外考虑时间点之间相隔一个完整季节周期（记为 S）的依赖关系。

纯季节性AR(1)模型

一个最简单的季节性模型是纯季节性AR(1)模型。其公式如下：
X_t = φ X_{t-S} + W_t
其中，W_t 是白噪声序列，S 是季节周期长度（例如，月度数据中 S=12）。

这个模型意味着当前值 X_t 仅依赖于 S 个时间单位之前的值 X_{t-S}，而与中间的所有值无关。从某种意义上说，这相当于同时运行 S 个相互独立的AR(1)过程。

一般季节性ARMA模型

更一般地，我们可以定义季节性ARMA(P’, Q’)_S模型。其公式使用后移算子 B 表示如下：
Φ_S(B^S) X_t = Θ_S(B^S) W_t
其中：

Φ_S(B^S) = 1 - Φ_1 B^S - Φ_2 B^{2S} - ... - Φ_{P'} B^{P'S} 是季节性自回归多项式。
Θ_S(B^S) = 1 + Θ_1 B^S + Θ_2 B^{2S} + ... + Θ_{Q'} B^{Q'S} 是季节性移动平均多项式。
P’ 和 Q’ 分别是季节性自回归和移动平均的阶数。

结合非季节性成分：SARIMA模型

通常，时间序列同时包含季节性成分和非季节性（短期）成分。我们可以将两者结合，得到乘积季节性ARIMA模型，记作 ARIMA(p,d,q)×(P’,D’,Q’)_S。其完整的模型方程为：
Φ_S(B^S) φ(B) (∇_S^{D’} ∇^d X_t) = Θ_S(B^S) θ(B) W_t
其中：

φ(B) 和 θ(B) 是非季节性ARIMA(p,d,q)模型中的多项式。
∇^d = (1-B)^d 是非季节性差分算子。
∇_S^{D’} = (1-B^S) 是季节性差分算子。
d 和 D’ 分别是非季节性和季节性差分的阶数。

这个模型结合了短期波动和长期季节性模式，参数较多，但能更全面地描述复杂的时间序列。

季节性模型的特性分析

为了理解季节性模型的行为，我们来分析两个具体例子。

示例1：纯季节性AR(1)_12模型

考虑模型 X_t = φ X_{t-12} + W_t，其中 |φ| < 1。

自协方差函数的计算：
通过递归代入，可以将 X_t 表示为白噪声的线性过程。计算其自协方差函数 γ(h) 可得：

γ(0) = σ_w² / (1 - φ²)
对于 h 不是12的整数倍的情况，γ(h) = 0。
对于 h = 12m (m=1,2,3...)，γ(12m) = σ_w² φ^{|m|} / (1 - φ²)。

关键结论：
该模型的自协方差仅在季节周期（12的倍数）的滞后处非零。这印证了其本质是多个独立AR(1)过程的叠加。在实践中的一个重要启示是：对于此类纯季节性数据，应直接拟合季节性模型，而非一个高阶的AR(12)模型，后者会引入大量不必要的参数。

示例2：混合季节性模型 ARIMA(0,1,1)×(1,0,0)_12

考虑一个结合了非季节性MA(1)和季节性AR(1)12的模型：
**X_t = φ X + θ W_{t-1} + W_t**

方差与自协方差的计算：
利用该模型的方程以及各项之间的不相关性，我们可以推导出：

方差 γ(0) = σ_w² (1 + θ²) / (1 - φ²)
在季节周期滞后处（h = 12m）：γ(12m) = σ_w² φ^{|m|} / (1 - φ²)，其中 m ≠ 0。
在紧邻季节周期的滞后处（h = 12m ± 1）：γ(12m ± 1) = σ_w² θ φ^{|m|} / (1 - φ²)。
对于其他滞后 h，γ(h) = 0。

自相关函数（ACF）的形态：
根据上述自协方差和方差，可以画出ACF图。它将呈现以下特征：

在滞后0处，自相关为1。
在滞后1处，有一个由 θ 决定的小峰。
在滞后2至11处，自相关接近于0。
在滞后12处，出现一个由 φ 决定的主峰。
在滞后11和13处，会出现较小的侧峰，其高度是滞后12处主峰高度的 θ 倍。
这种“主峰-侧峰”的模式会在24、36等季节周期倍数处重复出现，且峰值按 φ^m 几何衰减。

实践意义：
当我们对具有季节性的真实数据绘制样本ACF图时，如果观察到这种在固定周期间隔处出现峰值并几何衰减的模式，就是存在季节性成分的强烈信号。这为我们识别和建模季节性提供了直观的工具。

建模实践与软件实现

从数学上讲，季节性ARIMA模型是标准ARIMA模型的直接扩展。参数估计、模型选择（如使用AIC/BIC准则）、诊断检验和预测的所有原理仍然适用。

主要的挑战在于确定合适的季节周期 S。通常可以基于数据背景知识（如月度、季度、周度数据）进行选择。对于更复杂的情况，可以使用统计软件中的自动模型选择功能。

例如，在R语言的 forecast 包中，auto.arima() 函数可以自动尝试包含季节性成分的模型，并帮助确定一个合适的 S 值以及模型的阶数 (p,d,q)×(P’,D’,Q’)。

总结

本节课中我们一起学习了季节性ARIMA模型。我们从最简单的纯季节性AR模型入手，逐步介绍了通用季节性ARMA模型的表示方法，以及如何将其与非季节性ARIMA模型结合为SARIMA模型。通过分析两个具体例子的自协方差函数，我们理解了季节性模型在ACF图上表现出的独特模式——在季节周期及其附近出现峰值。最后，我们指出了建模实践中的关键点，即季节周期 S 的识别，并提及了可利用现有软件工具（如R中的auto.arima）来辅助完成这一过程。在接下来的课程中，我们将使用R软件对包含季节性的实际数据进行完整的分析与建模。

017：在R Studio中使用季节性ARIMA模型

在本节课中，我们将继续学习季节性ARIMA模型，但这次我们将使用R Studio进行实际操作。我们将利用之前课程中提到的处方药数据集，应用我们学到的所有工具，包括拟合ARMA过程、ARIMA过程、季节性模型以及线性模型等多种方法来分析时间序列数据。时间序列分析本质上是一个探索性的过程，我们需要花费大量时间拟合模型，以找出最适合我们数据集的方法。本节课我们将通过R Studio代码来实践这一过程。

数据加载与初步观察

首先，我们需要加载一些标准的时间序列分析包，并导入数据集。

# 加载必要的R包
library(TSA)
library(tseries)
library(forecast)

# 加载处方药数据集
data(prescription)

处方药数据集记录了从1986年8月到1992年3月期间，美国每月的平均处方药成本。通过绘制时间序列图，我们可以观察到数据有明显的上升趋势，并且似乎存在周期性的波动，这暗示数据中可能存在季节性成分。

平稳性检验

在拟合模型之前，我们需要检查时间序列的平稳性。我们使用增广迪基-富勒检验和菲利普斯-佩龙检验。

# 执行增广迪基-富勒检验
adf_test <- adf.test(prescription)
print(adf_test)

# 执行菲利普斯-佩龙检验
pp_test <- pp.test(prescription)
print(pp_test)

两个检验都给出了较小的p值，这意味着在去除线性趋势后，数据看起来是平稳的。因此，我们的第一步是使用线性模型去除数据中的线性趋势。

去除线性趋势

我们使用lm函数拟合一个线性模型，以去除数据中的线性趋势。

# 拟合线性模型去除趋势
linear_model <- lm(prescription ~ time(prescription))
prescription_detrended <- residuals(linear_model)

去除趋势后，我们得到了残差序列。接下来，我们需要检查这些残差是否还存在任何可识别的模式或自相关性。

残差分析

我们绘制去除趋势后残差的时序图、自相关函数图和偏自相关函数图。

# 绘制残差时序图
plot(prescription_detrended, main="去除线性趋势后的残差", ylab="残差")

# 绘制自相关函数图
acf(prescription_detrended, main="残差的ACF图")

# 绘制偏自相关函数图
pacf(prescription_detrended, main="残差的PACF图")

ACF图显示出明显的周期性波动，而PACF图在滞后1阶处有一个显著的尖峰，之后迅速衰减。这提示我们，残差中可能包含一个非季节性的AR(1)成分，同时ACF的周期性暗示可能存在季节性成分。

拟合非季节性ARIMA模型

基于PACF图的观察，我们首先尝试拟合一个非季节性的AR(1)模型。同时，为了比较，我们也拟合MA(1)和ARMA(1,1)模型。

# 拟合AR(1)模型
ar1_model <- arima(prescription_detrended, order=c(1,0,0))
print(summary(ar1_model))

# 拟合MA(1)模型
ma1_model <- arima(prescription_detrended, order=c(0,0,1))
print(summary(ma1_model))

# 拟合ARMA(1,1)模型
arma11_model <- arima(prescription_detrended, order=c(1,0,1))
print(summary(arma11_model))

通过比较AIC和BIC值，我们发现AR(1)模型在信息准则上表现最好，且更为简洁。因此，我们暂时选择AR(1)模型。

检查模型残差

拟合模型后，我们必须检查其残差是否像白噪声。我们提取AR(1)模型的残差并进行检验。

# 提取AR(1)模型残差
ar1_residuals <- residuals(ar1_model)

# 绘制残差ACF和PACF图
acf(ar1_residuals, main="AR(1)模型残差的ACF图")
pacf(ar1_residuals, main="AR(1)模型残差的PACF图")

尽管ACF和PACF图中的大部分尖峰都落在了置信区间内，但在滞后5、10等处仍能看到一些周期性模式的迹象。为了更正式地检验残差的自相关性，我们使用Ljung-Box检验。

# 执行Ljung-Box检验，考虑已拟合的参数
box_test_result <- Box.test(ar1_residuals, lag=5, type="Ljung-Box", fitdf=1)
print(box_test_result)

检验结果显示p值较小，表明残差中仍然存在显著的自相关性。这意味着我们的模型尚未完全捕捉数据中的所有结构，需要考虑加入季节性成分。

引入季节性成分

根据ACF图中在滞后5、10等处出现的尖峰，我们尝试为模型添加一个周期为5的季节性成分。我们拟合包含季节性AR(1)、季节性MA(1)以及季节性ARMA(1,1)的模型。

# 拟合非季节性AR(1) + 季节性AR(1)，周期为5
sar1_model <- arima(prescription_detrended,
                    order=c(1,0,0),
                    seasonal=list(order=c(1,0,0), period=5))
print(summary(sar1_model))

# 拟合非季节性AR(1) + 季节性MA(1)，周期为5
sma1_model <- arima(prescription_detrended,
                    order=c(1,0,0),
                    seasonal=list(order=c(0,0,1), period=5))
print(summary(sma1_model))

# 拟合非季节性AR(1) + 季节性ARMA(1,1)，周期为5
sarma_model <- arima(prescription_detrended,
                     order=c(1,0,0),
                     seasonal=list(order=c(1,0,1), period=5))
print(summary(sarma_model))

比较这些模型的AIC和BIC，我们发现包含季节性ARMA(1,1)成分的模型表现最佳。其残差的ACF和PACF图显示所有尖峰均落在置信区间内，Ljung-Box检验的p值也不再显著。这表明该模型已能很好地捕捉数据中的相关结构。

使用差分法替代线性回归

上一节我们使用线性回归去除趋势。另一种处理趋势的方法是直接对原始序列进行差分。我们计算处方药数据的一阶差分，并对其进行分析。

# 计算一阶差分
prescription_diff <- diff(prescription)

# 绘制差分后序列的ACF和PACF图
acf(prescription_diff, main="一阶差分序列的ACF图")
pacf(prescription_diff, main="一阶差分序列的PACF图")

差分后序列的ACF图仍然显示出在滞后5等处的周期性尖峰。我们尝试为差分后的序列拟合包含季节性成分的ARIMA模型。

自动模型选择

为了更高效地找到最优模型，我们可以使用forecast包中的auto.arima函数进行自动模型选择。

# 使用auto.arima自动选择最佳ARIMA模型
auto_model <- auto.arima(prescription, seasonal=TRUE, stepwise=FALSE, approximation=FALSE)
print(summary(auto_model))

auto.arima函数选择了一个模型：对原始数据做一阶差分，并包含一个周期为12的季节性AR(1)项。这个模型的AIC值很低，且其预测图能合理地延续数据的上升趋势和季节性波动，残差检验也表现良好。

包含漂移项的模型

最后，我们尝试拟合一个包含漂移项的ARIMA模型，而不是先进行差分。漂移项可以捕捉数据中的线性趋势。

# 拟合包含漂移项、非季节性AR(1)和季节性AR(1)的模型
drift_model <- arima(prescription,
                     order=c(1,0,0),
                     seasonal=list(order=c(1,0,0), period=5),
                     include.mean=TRUE,
                     xreg=time(prescription))
print(summary(drift_model))

该模型也给出了合理的拟合结果和预测。其估计出的月度增长率与最初线性模型得到的结果非常接近。

总结

在本节课中，我们一起探索了如何使用R Studio对时间序列数据拟合季节性ARIMA模型。我们从一个处方药成本数据集出发，完整地实践了数据分析流程：

数据观察与平稳性检验：识别趋势和季节性。
趋势处理：通过线性回归或差分方法去除趋势。
模型识别与拟合：依据ACF/PACF图选择模型阶数，并拟合非季节性和季节性ARIMA模型。
模型诊断：检查残差是否为白噪声，使用Ljung-Box等统计检验。
模型比较与选择：利用AIC、BIC等信息准则，并借助auto.arima进行自动化探索。
预测：使用选定模型进行未来值预测。

时间序列建模是一个迭代和探索的过程，通常没有唯一正确的答案。关键在于结合统计工具、图形化分析和领域知识，不断尝试和比较不同模型，直到找到一个既能充分捕捉数据特征、残差检验良好，又能做出合理预测的模型。希望本节课的实践能为你分析自己的时间序列数据提供一个清晰的思路框架。

018：周期过程与离散傅里叶变换

在本节课中，我们将开启课程最后一个主要章节的学习，即频域分析。我们将学习如何从频率的视角来审视时间序列数据，这对于识别数据中的周期性波动（如季节性模式）非常有用。我们将从周期过程的概念入手，并介绍一个强大的工具——离散傅里叶变换。

从时域到频域

到目前为止，我们所有的分析都集中在时域，即直接观察数据随时间的变化。然而，任何随时间变化的信号都可以从频域的角度进行分析。频域分析能帮助我们识别时间序列中哪些频率成分最为显著或“强大”。

周期过程

我们从一个基础的周期过程模型开始。它本质上是一种平稳过程，其形式如下：

X_t = A * cos(2πωt + φ)

其中：

A 是振幅，决定了信号的强度。
ω 是频率，例如，对于月度数据中的年度周期，ω = 1/12。
φ 是相位，决定了信号在时间轴上的偏移。

为了便于分析，我们利用三角恒等式将其重写为：

X_t = U₁ * cos(2πωt) + U₂ * sin(2πωt)

这里，我们将 U₁ 和 U₂ 视为均值为零、互不相关的随机变量。这使得 X_t 成为一个随机过程。如果我们进一步假设 U₁ 和 U₂ 的方差相同，都为 σ²，那么可以计算出该过程的自协方差函数：

γ(h) = Cov(X_{t+h}, X_t) = σ² * cos(2πωh)

这个结果非常有趣：自协方差只依赖于时间间隔 h，而与具体的时间点 t 无关，这证明了该过程是平稳的。同时，自协方差本身也以余弦波形式随 h 变化。

更一般的周期过程

为了使模型更贴近现实，我们可以考虑多个频率的叠加。一个更一般的周期过程可以表示为多个不同频率的余弦和正弦波之和：

X_t = Σ_{j=1}^{Q} [ U_{j1} * cos(2πω_j t) + U_{j2} * sin(2πω_j t) ]

其中，对于每个频率 ω_j，对应的随机变量 U_{j1} 和 U_{j2} 都是均值为零、互不相关且方差为 σ_j² 的。此时，总的自协方差函数是各个频率成分的自协方差之和：

γ(h) = Σ_{j=1}^{Q} σ_j² * cos(2πω_j h)

这样，我们的信号就可以包含多个不同强度的周期性成分。

⚠️ 警告：混叠现象
在离散采样时需要注意混叠现象。当采样频率不够高时，一个高频信号可能会在采样数据中表现为一个低频信号，导致频率识别错误。

回归与傅里叶变换

上一节我们介绍了周期过程模型，本节中我们来看看如何将其应用于实际数据。给定一个长度为 T 的时间序列数据 X₁, ..., X_T，我们可以尝试用一组余弦和正弦函数的线性组合来“拟合”或“表示”它。这本质上是一个线性回归问题。

对于 T 为奇数的情况，模型可写为：
X_t = β₀ + Σ_{j=1}^{(T-1)/2} [ β_{j1} * cos(2πt * j/T) + β_{j2} * sin(2πt * j/T) ]

对于 T 为偶数的情况，形式略有不同，会多出一项。

为什么这样做？因为余弦和正弦函数在一定条件下构成一组正交基。这意味着我们可以完美地（或近似地）用它们表示数据。更重要的是，系数 β_{j1} 和 β_{j2} 的大小直接反映了频率 j/T 在数据中的显著程度。

以下是求解这些回归系数的关键点：

设计矩阵（由所有余弦和正弦项构成）具有正交列。
这使得系数估计变得非常简单，无需进行复杂的矩阵求逆，可以直接计算：

β̂_{j1} = (2/T) * Σ_{t=1}^{T} X_t * cos(2πt * j/T)
β̂_{j2} = (2/T) * Σ_{t=1}^{T} X_t * sin(2πt * j/T) (对于 j 不为 0 或 T/2)
β̂₀ = X̄ (样本均值)

周期图

为了综合衡量某个频率 j/T 的强度，我们定义缩放周期图：

P(j/T) = β̂_{j1}² + β̂_{j2}²

P(j/T) 的值越大，表明频率 j/T 在时间序列中越显著。 直观上，它可以被理解为归属于该频率成分的方差。

离散傅里叶变换

直接按上述公式计算所有频率的系数和周期图在数据量大时计算量巨大。幸运的是，这可以通过离散傅里叶变换 高效完成。

DFT 的定义如下（其中 i 是虚数单位）：

D(j/T) = (1/√T) * Σ_{t=1}^{T} X_t * exp(-2πi * t * j / T)

利用欧拉公式 exp(iθ) = cos(θ) + i sin(θ)，可以将 DFT 分解为余弦和正弦部分。关键的联系在于，缩放周期图可以通过 DFT 的模平方来计算：

P(j/T) = (4/T) * |D(j/T)|²

因此，计算 DFT 就等价于计算我们之前回归分析所关心的所有频率成分的强度。

快速傅里叶变换

虽然 DFT 给出了计算方法，但直接计算一个长度为 T 的序列的 DFT 需要大约 O(T²) 次运算。对于大规模时间序列，这依然很慢。解决方案是使用快速傅里叶变换。

FFT 是一种高效算法，它利用了 DFT 计算中的对称性和冗余性。其核心思想是分而治之：

将长度为 T 的时间序列按奇偶索引分成两个子序列。
原始序列的 DFT 可以通过这两个子序列的 DFT 组合得到。
递归地对子序列进行同样的操作。

当 T 是 2 的整数次幂时，这种算法最为高效，能将计算复杂度从 O(T²) 降低到 O(T log₂ T)。现代计算库中的 FFT 算法还处理了 T 不是 2 的幂次的情况。

总结

本节课中我们一起学习了时间序列频域分析的基础：

我们引入了周期过程模型，它由不同频率的余弦和正弦波叠加而成，并且是平稳的。
我们展示了如何通过线性回归到一组正交的余弦/正弦基上来拟合数据，其系数揭示了数据中的显著频率。
我们定义了缩放周期图 P(j/T) 作为衡量频率显著性的工具。
我们介绍了离散傅里叶变换，它是计算周期图的等价且更概念化的方法。
最后，我们了解了快速傅里叶变换算法如何通过分而治之的策略，极大地提升了 DFT 的计算效率，使得对大规模时间序列进行频域分析成为可能。

在下一讲中，我们将深入探讨谱分布和谱密度的概念，并建立它们与我们熟悉的自协方差函数之间的重要联系。

019：谱密度与谱分布

在本节课中，我们将学习时间序列分析中谱方法的核心概念：谱密度与谱分布。我们将探讨它们如何与自协方差函数相关联，以及如何帮助我们识别时间序列数据中最显著的频率成分。

概述

上一节我们介绍了周期过程，其核心思想是将过程表示为正弦和余弦函数的组合，并通过系数大小来判断哪些频率在时间序列中占主导地位。本节中，我们将进一步深入，学习谱分布和谱密度这两个概念。它们的命名与概率论中的累积分布函数和概率密度函数类似。我们将发现，这两个概念与自协方差函数紧密相关，正如本课程一再强调的，自协方差是时间序列分析的基础。

谱分布与维纳-辛钦定理

首先，我们介绍维纳-辛钦定理。该定理建立了平稳时间序列的自协方差函数与谱分布之间的联系。

定理 1 (维纳-辛钦定理第一部分)：
令 {X_t} 为一个平稳时间序列，其自协方差函数为 K_X(h) = Cov(X_{t+h}, X_t)。那么，存在一个唯一的单调递增函数 F_X(ω)，称为谱分布函数，满足：

F_X(-1/2) = 0
F_X(1/2) = K_X(0) （即过程的方差）
对于任意滞后 h，自协方差函数可以表示为：
K_X(h) = ∫_{-1/2}^{1/2} e^{2πiωh} dF_X(ω)
这里的积分是黎曼-斯蒂尔杰斯积分。

这个定理表明，自协方差函数可以通过对复指数函数关于谱分布函数积分得到。如果谱分布函数 F_X(ω) 是可微的，那么我们可以定义其导数。

谱密度

上一节我们介绍了谱分布的存在性，本节我们来看看在什么条件下可以定义其导数，即谱密度。

定理 2 (维纳-辛钦定理第二部分)：
令 {X_t} 为一个平稳时间序列，其自协方差函数 K_X(h) 是绝对可和的，即：
∑_{h=-∞}^{∞} |K_X(h)| < ∞
那么，谱分布函数 F_X(ω) 是绝对连续的，其导数 f_X(ω) = dF_X(ω)/dω 几乎处处存在，并称为谱密度函数。此时，自协方差函数可表示为：
K_X(h) = ∫_{-1/2}^{1/2} e^{2πiωh} f_X(ω) dω
更重要的是，我们有一个逆变换，可以直接从自协方差得到谱密度：
f_X(ω) = ∑_{h=-∞}^{∞} K_X(h) e^{-2πiωh}，其中 ω ∈ [-1/2, 1/2]。

谱密度具有以下重要性质：

定义域：ω 在 [-1/2, 1/2] 区间内。
偶函数：f_X(-ω) = f_X(ω)。这是因为自协方差函数是偶函数，且复指数项在求和时对称。
方差分解：当 h=0 时，K_X(0) = Var(X_t) = ∫_{-1/2}^{1/2} f_X(ω) dω。这意味着过程的方差可以分解为不同频率成分的贡献之和，类似于方差分析中将总平方和分解为各因素平方和的思想。

示例分析

为了更直观地理解，让我们看两个简单的例子。

示例 1：周期过程
考虑过程 X_t = U_1 cos(2πω_0 t) + U_2 sin(2πω_0 t)，其中 U_1, U_2 是不相关的零均值随机变量。其自协方差为 K_X(h) = σ² cos(2πω_0 h)。这个自协方差序列不是绝对可和的，因此不存在谱密度函数，但存在谱分布函数 F_X(ω)。它是一个阶梯函数，在频率 ±ω_0 处有跳跃。这直观地反映了该过程的所有“能量”都集中在频率 ω_0 上。

示例 2：白噪声
白噪声过程 W_t 的自协方差为：K_W(0) = σ²，K_W(h) = 0 (当 h ≠ 0)。这显然是绝对可和的。其谱密度为：
f_W(ω) = ∑_{h=-∞}^{∞} K_W(h) e^{-2πiωh} = σ²
这是一个常数函数。表明白噪声在所有频率上具有均匀的“能量”密度，这解释了其名称中“白”的由来（类比白光包含所有颜色的光）。

滤波与ARMA过程的谱密度

我们已经知道如何从自协方差计算谱密度。但对于像ARMA这样的常见模型，是否有更简便的方法呢？答案是肯定的，这涉及到滤波的概念。

考虑一个“过滤”过程：Y_t = ∑_{j=-∞}^{∞} a_j X_{t-j}，其中系数 {a_j} 绝对可和。定义脉冲响应函数 A(j) = a_j 及其频率响应函数（傅里叶变换）：
A(ω) = ∑_{j=-∞}^{∞} a_j e^{-2πiωj}

定理 3 (滤波的谱密度)：
如果 {X_t} 具有谱密度 f_X(ω)，那么过滤后过程 {Y_t} 的谱密度为：
f_Y(ω) = |A(ω)|² f_X(ω)
这意味着，滤波操作对谱密度的影响，仅仅是乘以频率响应函数模的平方。

这个定理可以优雅地应用于ARMA过程。任何一个因果可逆的ARMA(p,q)过程 φ(B)X_t = θ(B)W_t 都可以表示为白噪声的无限滑动平均形式：X_t = ∑_{j=0}^{∞} ψ_j W_{t-j}。这里，系数 {ψ_j} 对应滤波系数 {a_j}，而 W_t 是白噪声。频率响应函数 Ψ(ω) 恰好是移动平均算子与自回归算子在 z = e^{-2πiω} 处取值的比：
Ψ(ω) = θ(e^{-2πiω}) / φ(e^{-2πiω})
由于白噪声的谱密度是常数 σ²，根据定理3，ARMA过程的谱密度为：
f_X(ω) = |Ψ(ω)|² σ² = |θ(e^{-2πiω}) / φ(e^{-2πiω})|² σ²
这个公式让我们无需计算复杂的自协方差，直接通过模型参数就能得到谱密度。

实际应用与代码演示

在理论介绍之后，让我们通过R代码来看一个实际例子。我们生成一个包含多个特定频率成分的周期信号，并加入一些白噪声。

set.seed(479)
n <- 512
t <- 0:(n-1)
signal <- sin(2*pi*32/n * t) + 0.5*cos(2*pi*16/n * t) + 0.8*sin(2*pi*16/n * t + pi/6) + 1.2*cos(2*pi*96/n * t)
x <- signal + rnorm(n, sd=0.5)

# 计算周期图 (谱密度的估计)
xx_fft <- fft(x)
xx_periodogram <- (2/n) * Mod(xx_fft)^2
freq <- (0:(n-1))/n

# 绘制周期图
plot(freq, xx_periodogram, type='l', main='Periodogram', xlab='Frequency', ylab='Intensity')
abline(v=0.5, lty=2) # 标记 Nyquist 频率 (1/2)

# 找出最强的频率成分
sorted_periodogram <- sort(xx_periodogram, decreasing=TRUE, index.return=TRUE)
top_freq_indices <- sorted_periodogram$ix[1:16]
top_cycle_lengths <- 1/freq[top_freq_indices]
print(top_cycle_lengths[1:8])

运行代码后，我们会在周期图中观察到几个突出的尖峰。通过识别这些尖峰对应的频率（或周期长度），我们可以发现数据中隐藏的周期性成分（例如，对应周期长度约为16、32和5.33）。即使数据被噪声污染，主要的频率成分依然可以通过谱分析被检测出来。

总结

本节课我们一起学习了时间序列谱分析的核心：谱分布与谱密度。

我们首先通过维纳-辛钦定理理解了谱分布与自协方差函数之间的等价关系。
在自协方差绝对可和的条件下，我们定义了谱密度 f_X(ω)，它代表了方差在不同频率上的分布，并且是偶函数。
通过滤波定理，我们掌握了如何计算线性过滤后过程的谱密度，并推导出ARMA过程谱密度的简洁公式：f_X(ω) = σ² |θ(e^{-2πiω}) / φ(e^{-2πiω})|²。
最后，通过代码演示，我们看到了如何从实际数据中估计谱密度（周期图），并识别出主导的频率成分。

谱密度为我们分析时间序列提供了全新的视角，即从“时域”转换到“频域”。它特别擅长揭示数据中的周期性、季节性以及各种频率的振荡模式。在接下来的课程中，我们将以此为基础，学习如何进行谱估计和统计推断。

020：谱密度估计 🎯

在本节课中，我们将学习如何利用谱方法进行实际的统计分析，核心目标是估计谱密度函数。我们将发现，我们已经掌握的工具——周期图——正是这个估计器。本节课将深入探讨其细节，并介绍谱方差分析以及谱密度的置信区间构建方法。

离散傅里叶变换回顾 📊

上一节我们介绍了谱密度的数学概念，但在实际中，我们面对的是数据而非数学模型。因此，我们需要从数据出发来估计谱密度。这一切的起点是离散傅里叶变换。

给定一个时间序列数据 ( x_1, x_2, \ldots, x_T )，其离散傅里叶变换定义为：
[
d(\omega_j) = \frac{1}{\sqrt{T}} \sum_{t=1}^{T} x_t e^{-2\pi i \omega_j t}
]
其中，频率 ( \omega_j = j / T )，( j = 0, 1, \ldots, T-1 )。

这个变换可以分解为实部和虚部，分别对应余弦变换和正弦变换：

余弦变换：( d_c(\omega_j) = \frac{1}{\sqrt{T}} \sum_{t=1}^{T} x_t \cos(2\pi \omega_j t) )
正弦变换：( d_s(\omega_j) = \frac{1}{\sqrt{T}} \sum_{t=1}^{T} x_t \sin(2\pi \omega_j t) )

它们的关系是：
[
d(\omega_j) = d_c(\omega_j) - i \cdot d_s(\omega_j)
]

离散傅里叶变换是可逆的，我们可以通过逆变换从频域数据恢复原始时域数据：
[
x_t = \sum_{j=0}^{T-1} d(\omega_j) e^{2\pi i \omega_j t}, \quad t = 1, \ldots, T
]

周期图作为谱密度估计器 🔍

现在，我们来看看如何从数据估计谱密度。核心工具是周期图。

对于频率 ( \omega_j )，周期图定义为离散傅里叶变换模的平方：
[
I(\omega_j) = |d(\omega_j)|^2 = [d_c(\omega_j)]^2 + [d_s(\omega_j)]^2
]

通过一系列代数推导（包括对数据进行中心化处理，即减去均值 ( \bar{x} )），我们可以将周期图表示为样本自协方差函数的余弦变换：
[
I(\omega_j) = \sum_{h=-(T-1)}^{T-1} \hat{\gamma}(h) \cos(2\pi \omega_j h)
]
其中，( \hat{\gamma}(h) ) 是在滞后 ( h ) 处的样本自协方差函数估计。

这个形式让我们联想到维纳-辛钦定理。该定理指出，对于一个平稳过程，其谱密度 ( f(\omega) ) 是其真实自协方差函数 ( \gamma(h) ) 的傅里叶变换：
[
f(\omega) = \sum_{h=-\infty}^{\infty} \gamma(h) e^{-2\pi i \omega h}
]

结论：周期图 ( I(\omega_j) ) 是谱密度 ( f(\omega) ) 的一个估计器。它用有限的样本自协方差 ( \hat{\gamma}(h) ) 代替了无限的真实自协方差 ( \gamma(h) )。

然而，对于较大的滞后 ( h )，样本自协方差 ( \hat{\gamma}(h) ) 的估计质量会变差（因为用于估计的数据点变少）。因此，在实践中，我们常使用一个截断版本，只对 ( |h| \leq M ) 的滞后进行求和，其中 ( M \ll T )：
[
I(\omega_j) \approx \sum_{h=-M}^{M} \hat{\gamma}(h) \cos(2\pi \omega_j h)
]

谱方差分析 📈

从方差分析的角度，我们可以将时间序列的总变异分解到不同频率成分上。这提供了理解周期图的另一种视角。

我们可以将中心化后的数据 ( x_t - \bar{x} ) 表示为一组余弦和正弦函数的线性组合（一个饱和模型）。通过最小二乘法，可以精确拟合数据，其系数正是我们之前定义的余弦变换和正弦变换。

当我们计算数据的总平方和（总变异）时，一个关键的性质是：不同频率对应的余弦和正弦基函数是正交的。这使得交叉项在求和时为零。最终，总平方和可以分解为各个频率成分的贡献之和：
[
\sum_{t=1}^{T} (x_t - \bar{x})^2 = 2 \sum_{j=1}^{(T-1)/2} I(\omega_j)
]

这意味着：我们将数据的总变异（总平方和）分解为一系列频率 ( \omega_j ) 对应的平方和。每个频率成分 ( I(\omega_j) ) 可以视为拥有2个自由度的平方和。这与经典方差分析的思想一致，只不过这里的“处理”是不同频率的周期成分。

注意：虽然数学上等价，但在实际计算中，我们不会真的去解这个巨大的线性回归系统，而是使用高效的快速傅里叶变换算法。

大样本性质与置信区间 📉

理解了周期图是谱密度的估计器后，一个自然的问题是：这个估计器有多好？我们在大样本下（即数据量 ( T \to \infty )）研究其性质。

我们假设时间序列 ( X_t ) 是一个平稳过程，具有绝对可和的自协方差函数，因此谱密度存在。

期望：首先，在一定的正则条件下，可以证明周期图的期望值渐近地趋近于真实的谱密度：
[
\mathbb{E}[I(\omega_j)] \to f(\omega) \quad \text{当 } T \to \infty，且 \omega_j \to \omega
]

分布：为了构建置信区间，我们需要知道周期图的抽样分布。对于一个线性过程（例如ARMA过程），在更强的自协方差衰减条件下，我们可以得到一个极限定理。

具体来说，对于一系列趋向于固定频率 ( \omega ) 的离散频率 ( \omega_{j,T} )，有以下联合收敛分布结果：
[
\frac{2 I(\omega_{j,T})}{f(\omega)} \xrightarrow{d} \chi^2_2
]
其中，( \chi^2_2 ) 表示自由度为2的卡方分布。

这个定理非常强大，因为它允许我们为谱密度构造置信区间。

置信区间：基于上述极限分布，我们可以推导出谱密度 ( f(\omega) ) 的近似 ( 1-\alpha ) 置信区间。对于某个我们感兴趣的频率 ( \omega )，选择数据中与之最接近的离散频率 ( \omega_j )，则置信区间为：
[
\left[ \frac{2 I(\omega_j)}{\chi^2_{2, 1-\alpha/2}}, \ \frac{2 I(\omega_j)}{\chi^2_{2, \alpha/2}} \right]
]
这里，( \chi^2_{2, \alpha/2} ) 和 ( \chi^2_{2, 1-\alpha/2} ) 是自由度为2的卡方分布的上侧分位数。

总结 🎓

本节课中，我们一起学习了谱密度估计的核心统计方法。

估计器：我们明确了周期图 ( I(\omega_j) ) 是谱密度 ( f(\omega) ) 的自然估计器，它本质上是样本自协方差函数的傅里叶变换。
方差分析视角：我们从方差分析的角度理解了周期图，它将数据的总变异分解到不同频率上，每个频率成分贡献一个具有2个自由度的平方和。
统计推断：我们探讨了周期图的大样本性质。最重要的是，对于线性过程（如ARMA模型），标准化后的周期图渐近服从卡方分布。这为我们构建谱密度的置信区间提供了理论基础。

通过本节课的学习，我们不仅知道了如何用周期图估计谱密度，还掌握了评估这种估计不确定性的统计工具。在接下来的课程中，我们将继续探讨更复杂的谱估计方法，并学习如何在R语言中实现这些分析。

021：谱密度估计 II

在本节课中，我们将继续讨论谱密度的估计方法。上一讲我们介绍了周期图，它是一种有效的估计器，但仍有改进空间。本节我们将探讨多种平滑和优化估计的技术，包括Bartlett方法、Welsh方法、频带平均、加权平均、锥化以及相关核函数的概念。这些方法旨在通过减少方差或聚焦于特定频率来获得更好的谱密度估计。

Bartlett方法

上一节我们介绍了周期图，本节我们来看看如何通过Bartlett方法来改进估计。Bartlett方法的核心思想是将时间序列分割成多个不重叠的等长子序列，分别计算每个子序列的周期图，然后对这些周期图进行平均。

对于一个长度为 ( T ) 的时间序列 ( X_1, \ldots, X_T )，我们将其分割成 ( K ) 个长度为 ( M ) 的不重叠子序列，其中 ( M = T/K )。每个子序列可以表示为：

第1段：( X_1, \ldots, X_M )
第2段：( X_{M+1}, \ldots, X_{2M} )
...
第K段：( X_{T-M+1}, \ldots, X_T )

对于每个子序列，我们计算其离散傅里叶变换（DFT）并得到周期图 ( I^{(i)}(\omega_j) )，其中 ( i = 1, \ldots, K )。Bartlett估计器定义为这些周期图的平均值：

[
\bar{I}(\omega_j) = \frac{1}{K} \sum_{i=1}^{K} I^{(i)}(\omega_j)
]

使用Bartlett方法的效果是，我们考虑的频率点变少了（因为 ( \omega_j = j/M )），但估计的方差降低了。这是因为我们对 ( K ) 个独立估计进行了平均。此外，计算 ( K ) 个较小的DFT通常比计算一个大型DFT更快，且易于并行化。

Welsh方法

Welsh方法与Bartlett方法思路相似，但它允许子序列之间存在重叠。这意味着在分割时间序列时，各段可以部分重叠，从而可能利用更多数据窗口来获得更平滑的估计。虽然具体实现略有不同，但其核心目标与Bartlett方法一致：通过平均多个估计来降低方差。

频带平均（Banding）

接下来我们探讨在频域进行平均的方法，即频带平均。与在时域分割数据不同，频带平均是将频率 ( \omega ) 划分成多个频带，然后在每个频带内对周期图值进行平均。

我们定义一个以 ( \omega_j ) 为中心、带宽为 ( (2M+1)/T ) 的频带 ( B )：

[
B = \left{ \omega : \omega_j - \frac{M}{T} \le \omega \le \omega_j + \frac{M}{T} \right}
]

其中，( 2M+1 ) 是频带内包含的离散频率点数。带宽 ( (2M+1)/T ) 表示该频带覆盖的单位区间比例。

对于给定的频带 ( B \，我们计算平均周期图：

[
\bar{I}(\omega_j) = \frac{1}{2M+1} \sum_{i=-M}^{M} I\left(\omega_j + \frac{i}{T}\right)
]

这种方法的优点是，如果谱密度在频带 ( B ) 内近似为常数，那么通过平均可以减少估计的方差。从理论上讲，经过适当标准化后，( 2(2M+1) \bar{I}(\omega_j) / f_X(\omega_j) ) 会依分布收敛于自由度为 ( 4M+2 ) 的卡方分布，据此可以构建置信区间。

加权频带平均

我们还可以对频带内的周期图值进行加权平均，而不仅仅是简单平均。这引入了更大的灵活性。

定义加权平均周期图：

[
\tilde{I}(\omega_j) = \sum_{i=-M}^{M} c_i I\left(\omega_j + \frac{i}{T}\right)
]

其中，权重系数 ( c_i ) 满足 ( \sum_{i=-M}^{M} c_i = 1 )。当所有权重相等，即 ( c_i = 1/(2M+1) ) 时，就退化为简单的频带平均。

在满足一定条件（如 ( M/T \to 0 ) 且 ( \sum c_i^2 \to 0 )）时，该估计量是渐近无偏的，并且不同频率处的估计值渐近不相关。其极限分布是一个加权卡方和，但可以近似看作自由度为 ( 2L ) 的卡方分布，其中 ( L = 1 / \sum_{i=-M}^{M} c_i^2 )。当权重相等时，( L = 2M+1 )，与之前的结果一致。

核函数选择：Daniel核

如何选择权重系数 ( c_i ) 呢？一个常见的选择是使用Daniel核。其思想是反复对序列应用均匀平滑操作。

一次应用：产生权重为 ( [1/3, 1/3, 1/3] ) 的均匀核。
二次应用：相当于对均匀核进行卷积，产生形如 ( [1/9, 2/9, 3/9, 2/9, 1/9] ) 的三角权重。
多次应用：会产生更高阶（如抛物线形）的平滑权重。

这类似于均匀分布的反复卷积。在R中，可以通过指定平滑次数来控制核的形状。此外，还有修正Daniel核，它初始就给中心点更高的权重（例如 ( [1/4, 2/4, 1/4] )），从而更快地形成平滑效果。

锥化（Tapering）

前面介绍的方法主要在频域进行操作。现在我们将视角转回时域，介绍锥化技术。锥化是指对原始时间序列数据本身施加一个权重函数（锥化窗口），通常给予序列中间部分较高的权重，两端部分较低的权重。

设 ( X_t ) 是零均值平稳时间序列，其谱密度为 ( f_X(\omega) )。定义锥化后的序列 ( Y_t = a_t X_t )，其中 ( a_t ) 是实数权重系数。( Y_t ) 的离散傅里叶变换为：

[
D_Y(\omega_j) = \frac{1}{\sqrt{T}} \sum_{t=1}^{T} a_t X_t e^{-2\pi i \omega_j t}
]

可以证明，锥化后序列 ( Y_t ) 的周期图的期望值不再直接等于 ( f_X(\omega_j) )，而是原谱密度与一个谱窗函数 ( W_T(\omega) ) 的卷积：

[
\mathbb{E}[I_Y(\omega_j)] = \int_{-1/2}^{1/2} W_T(\omega_j - \lambda) f_X(\lambda) d\lambda
]

其中，谱窗 ( W_T(\omega) = |A(\omega)|^2 )，而 ( A(\omega) ) 是权重序列 ( {a_t} ) 的离散傅里叶变换。这意味着，对时域数据锥化的效果，在频域上相当于用一个特定的窗函数对谱密度进行局部平滑。

一个经典的例子是Fejér或Bartlett核（注意与之前的Bartlett方法区分），它对应于不对数据做任何锥化（即 ( a_t = 1 )）时产生的谱窗，其形式为：

[
W_T(\omega) = \frac{\sin^2(T \pi \omega)}{T \sin^2(\pi \omega)}
]

有趣的是，之前讨论的频带平均方法，也可以等价地表示为使用某种特定谱窗的锥化操作，这体现了时域和频域方法之间的深刻联系。

锥化通常能改善估计性能，因为它减少了因时间序列截断（两端突然变为0）而引起的频谱泄漏。一个常用的锥化函数是余弦锥，在时域定义为：

[
a_t = \frac{1}{2} \left[ 1 + \cos\left( \frac{2\pi (t - (T+1)/2)}{T} \right) \right]
]

这个函数在序列中心达到最大值1，向两端逐渐平滑地衰减到0，从而强调中间数据，削弱边界数据的影响。

总结

本节课我们一起学习了多种谱密度估计的改进方法。

Bartlett和Welsh方法通过在时域分割并平均多个周期图来降低方差。
频带平均及其加权版本则在频域对相邻频率的周期图进行平均或加权平均，同样旨在减少方差，并可通过选择核函数（如Daniel核）来调整平滑方式。
锥化技术通过给原始时间序列施加一个两端衰减的权重窗口，在时域进行预处理，其频域效应等价于用谱窗对谱密度进行平滑，有助于减少频谱泄漏。

所有这些方法都没有绝对的“最佳”选择，它们提供了从不同角度审视和估计谱密度的工具。在实际分析中，我们常常需要尝试多种方法并调整参数（如带宽 ( M )、权重核、锥化函数等），以探索数据中隐藏的周期特性。下一讲我们将通过R语言代码来实践这些方法，直观地观察它们的效果。

022：谱密度估计 III 🎯

在本节课中，我们将学习如何使用参数化方法来估计时间序列的谱密度。到目前为止，我们主要使用了周期图及其平滑版本，这些都属于非参数统计的范畴。然而，如果我们假设数据来自一个自回归（AR）或自回归滑动平均（ARMA）过程，我们就可以通过拟合模型来直接计算其谱密度。本节课我们将深入探讨这种方法，并通过R代码进行实践。

参数化谱密度估计 📈

上一节我们介绍了非参数化的谱密度估计方法。本节中，我们将探讨一种不同的思路：参数化估计。

其核心思想是：与其直接从数据计算周期图，不如先为数据拟合一个参数化时间序列模型（如ARMA模型），然后基于拟合好的模型来计算其理论谱密度。

具体步骤如下：

为观测数据拟合一个ARMA(p, q)模型，得到参数估计值 (\hat{\phi}) 和 (\hat{\theta})。
使用这些估计参数，计算该ARMA过程的谱密度。

一个更常见且在实践中易于实现的方法是：仅拟合一个自回归（AR）模型。R语言中的许多函数默认采用此方法。

AR模型谱密度公式

假设我们为数据拟合了一个AR(p)模型，并得到了参数估计 (\hat{\phi}_1, ..., \hat{\phi}_p) 以及白噪声方差估计 (\hat{\sigma}^2)。那么，在频率 (\omega) 处的谱密度估计为：

[
\hat{f}x(\omega) = \frac{\hat{\sigma}^2}{|1 - \sum^{p} \hat{\phi}_k e^{-2\pi i k \omega}|^2}
]

这个公式与我们之前推导的AR过程理论谱密度公式一致，只是将真实参数替换为了其估计值。

置信区间

与周期图方法类似，参数化估计也能给出谱密度的渐近置信区间。对于较大的样本量 (T)，有：

[
\frac{\hat{f}x(\omega)}{1 + \sqrt{2p/T} , z{1-\alpha/2}} \leq f_x(\omega) \leq \frac{\hat{f}x(\omega)}{1 - \sqrt{2p/T} , z{\alpha/2}}
]

其中 (z_{\alpha}) 是标准正态分布的分位数。这与非参数方法中使用的卡方分布置信区间形成了对比。

理论基础与注意事项

参数化方法有一个强大的理论支撑：对于任何平稳过程，都可以找到一个有限阶的因果AR过程，使其谱密度任意接近原过程的谱密度。这意味着，即使真实过程不是AR过程，用AR模型来近似其谱密度在理论上是可行的。

然而，需要注意以下几点：

上述结果是渐近的，要求样本量 (T \to \infty)。
同时，自回归阶数 (p) 的增长速度必须满足 (p^3 / T \to 0)，即不能增长过快。
在实践中，即使真实过程是ARMA，其可逆表示可以写成一个无限阶AR过程，但系数通常会快速衰减。因此，用一个中等阶数（如5或10）的AR模型来近似通常是足够的。

R语言实践 💻

理论部分之后，让我们进入RStudio，看看如何实际进行谱密度估计，并比较参数化与非参数化方法。

首先，我们需要生成一些示例数据来进行分析。

# 设置随机种子以保证结果可重复
set.seed(2479)

# 加载必要的包
library(TSA)
library(astsa)
library(forecast)

# 模拟数据
# 1. AR(1) 过程，系数为 0.8
ar1 <- arima.sim(model = list(ar = 0.8), n = 500)
# 2. 季节性 AR(4) 过程 (系数仅在滞后4处非零)
ar4_seasonal <- arima.sim(model = list(ar = c(0, 0, 0, 0.8)), n = 500)
# 3. MA(1) 过程，系数为 0.8
ma1 <- arima.sim(model = list(ma = 0.8), n = 500)
# 4. 季节性 MA(4) 过程
ma4_seasonal <- arima.sim(model = list(ma = c(0, 0, 0, 0.8)), n = 500)

非参数化估计：`spectrum()` 与 `spec.pgram()`

R基础包中的 spectrum() 函数（或直接使用 spec.pgram()）提供了非参数化的谱密度估计。默认使用修正的Daniel平滑器对周期图进行平滑。

spans 参数用于控制平滑窗的宽度。数值越大，平滑程度越高。

# 对AR1数据使用非参数化估计
# 默认设置，平滑较少
spec_raw <- spectrum(ar1, plot = FALSE)
# 使用 spans 参数增加平滑 (例如，使用宽度为31的窗)
spec_smooth <- spectrum(ar1, spans = 31, plot = FALSE)

# 绘制对比
plot(spec_raw$freq, spec_raw$spec, type = "l", col = "blue",
     main = "非参数化谱估计对比 (AR1数据)",
     xlab = "频率", ylab = "谱密度 (对数尺度)")
lines(spec_smooth$freq, spec_smooth$spec, col = "red", lwd = 2)
legend("topright", legend = c("原始周期图", "平滑后 (spans=31)"),
       col = c("blue", "red"), lty = 1, lwd = c(1,2))

参数化估计：`spec.ar()`

spec.ar() 函数通过拟合AR模型来进行参数化谱估计。默认使用AIC准则自动选择AR模型的阶数 (p)。

# 对AR1数据进行参数化谱估计
spec_ar <- spec.ar(ar1, plot = FALSE)
# 绘制结果
plot(spec_ar$freq, spec_ar$spec, type = "l", lwd = 2, col = "darkgreen",
     main = "参数化(AR)谱估计 (AR1数据)",
     xlab = "频率", ylab = "谱密度")

方法比较与叠加

将参数化估计的结果叠加到非参数化估计的图上，可以直观比较。

# 在同一图中比较非参数平滑估计与参数化估计
spectrum(ar1, spans = 31, main = "AR1数据谱估计对比",
         xlab = "频率", ylab = "谱密度 (对数尺度)")
# 添加参数化估计曲线
lines(spec.ar(ar1, plot = FALSE)$freq,
      spec.ar(ar1, plot = FALSE)$spec,
      col = "red", lwd = 2)
legend("topright", legend = c("非参数 (spans=31)", "参数化 (AR)"),
       col = c("black", "red"), lty = 1, lwd = 2)

分析不同数据生成过程

让我们将这两种方法应用于我们生成的所有数据集，观察其表现。

以下是针对不同数据的观察：

AR(1) 数据：参数化方法 (spec.ar) 能给出非常平滑且准确的估计，因为它正确指定了模型形式。非参数方法需要合适的平滑才能接近真实形状。
季节性 AR(4) 数据：参数化方法成功捕捉到了季节频率（0.25和0.5 cycles/time unit）处的峰值。未经平滑的非参数周期图非常嘈杂，而过度平滑可能会抹平这些重要峰值。
MA(1) 数据：参数化方法需要拟合一个相对高阶的AR模型来近似MA过程，其估计曲线可能不够平滑或包含虚假波动。非参数方法通过适度平滑，有时能更好地展现真实的谱形状（MA(1)的谱通常较为平坦或单调）。
季节性 MA(4) 数据：情况类似，参数化估计可能产生复杂的波动模式。需要谨慎解释这些波动是真实特征还是模型近似带来的伪影。

查看真实谱密度

对于ARMA模型，我们可以使用 armaspec 函数（来自 TSA 包）计算其真实的理论谱密度，以便与估计值比较。

# 计算 AR(1) 系数为 0.8 的真实谱密度
true_spec_ar1 <- armaspec(model = list(ar = 0.8), plot = FALSE)
# 计算并绘制估计谱与真实谱
spec_est <- spec.ar(ar1, plot = FALSE)
plot(true_spec_ar1$freq, true_spec_ar1$spec, type = "l", lwd = 3,
     main = "AR1: 真实谱 vs 参数化估计谱",
     xlab = "频率", ylab = "谱密度")
lines(spec_est$freq, spec_est$spec, col = "red", lwd = 2, lty = 2)
legend("topright", legend = c("真实谱密度", "估计谱密度 (AR)"),
       col = c("black", "red"), lty = c(1, 2), lwd = c(3,2))

总结与要点 🎓

本节课我们一起学习了时间序列谱密度估计的参数化方法。

核心思路：通过为数据拟合一个AR（或ARMA）模型，然后利用模型的理论公式计算谱密度。
主要函数：在R中，spec.ar() 用于参数化估计，spectrum() 或 spec.pgram() 用于非参数化估计。
方法比较：
- 参数化方法 在模型指定正确时效率高、结果平滑，并能提供理论上的置信区间。即使模型不完全正确，AR模型也能很好地近似许多平稳过程的谱密度。
- 非参数化方法 更灵活，不依赖于特定的模型假设。但其结果严重依赖于平滑程度（spans参数）的选择，面临着偏差-方差权衡：平滑不足则估计方差大、曲线嘈杂；过度平滑则会引入偏差，可能掩盖真实的谱特征（如尖锐峰值）。
实践建议：没有一种方法总是最优。建议在实际分析中同时尝试多种方法（参数化、不同平滑程度的非参数化），并将结果进行对比。通过观察谱密度图中的峰值、谷值和平坦区域，可以推断时间序列中存在的周期性、季节性或白噪声成分。

参数化谱估计为我们提供了另一种强大的工具来理解和揭示时间序列中的频率域特征。结合非参数方法，我们可以更全面、更稳健地进行分析。

023：指数平滑法 📈

在本节课中，我们将学习一种名为指数平滑的时间序列建模与预测方法。我们将探讨其基本概念、与ARIMA模型的联系，并通过R语言中的实际案例来应用这些方法。

概述 📋

指数平滑法是一种广泛使用的预测技术，尤其适用于没有明显季节性但可能包含趋势的数据。与之前课程中深入讨论的ARIMA模型不同，指数平滑法采用了一种加权平均的思维方式，对近期观测值赋予更高的权重。本节课我们将介绍简单指数平滑、双指数平滑（霍尔特-温特斯法）的基本原理，并了解它们如何与特定的ARIMA模型相关联。

简单指数平滑 📉

上一节我们回顾了课程初期提到的移动平均平滑器。本节中我们来看看一种不同的平滑方法：简单指数平滑。

简单指数平滑通过一个递归公式来创建平滑序列。其核心思想是，新的平滑值是当前观测值与上一期平滑值的加权平均。

公式如下：
S_t = α * X_t + (1 - α) * S_{t-1}
其中：

S_t 是时间 t 的平滑值。
X_t 是时间 t 的原始观测值。
α 是平滑参数，取值范围在0到1之间。
初始值通常设为 S_1 = X_1。

参数 α 控制平滑程度：

α 接近1时，模型几乎完全依赖最新观测值，平滑效果弱。
α 接近0时，模型赋予历史平滑值很大权重，结果更为平滑。

这种方法之所以被称为“指数”平滑，是因为如果将公式迭代展开，会发现每个历史观测值 X_{t-j} 的权重按 (1-α)^j 呈指数衰减。因此，它是一个使用全部历史数据的加权移动平均平滑器。

从信号处理的角度看，这种平滑相当于一个低通滤波器，旨在滤除时间序列中的高频随机波动。

双指数平滑（霍尔特线性趋势法）📈

简单指数平滑适用于没有趋势的平稳序列。如果数据包含趋势（持续上升或下降），人们通常使用双指数平滑，也称为霍尔特线性趋势法。

这种方法在简单指数平滑的基础上，引入了一个额外的方程来捕捉和更新趋势成分。

以下是模型的初始化与更新方程：

初始化：
S_1 = X_1
B_1 = X_2 - X_1 （初始趋势估计，可理解为斜率）

平滑方程：
S_t = α * X_t + (1 - α) * (S_{t-1} + B_{t-1})
B_t = β * (S_t - S_{t-1}) + (1 - β) * B_{t-1}

其中：

S_t 是截距项（序列水平）。
B_t 是斜率项（趋势）。
α 和 β 都是介于0和1之间的平滑参数。

B_t 的方程本身也是一个指数平滑形式，它平滑的是连续两期水平值 (S_t - S_{t-1}) 的变化量。通过这种方式，模型能够动态地估计并更新数据中的趋势。

与ARIMA模型的联系 🔗

一个有趣的事实是，简单指数平滑在数学上等价于一个特定的ARIMA模型。

结论：简单指数平滑对应于 ARIMA(0,1,1) 模型。

让我们来验证一下。ARIMA(0,1,1)模型的形式为：
(1 - B) X_t = (1 + θ B) W_t
其中 B 是后移算子，W_t 是白噪声。

将其展开并重新排列，可以得到预测公式：
X_t = X_{t-1} + W_t + θ W_{t-1}

在进行一步向前预测时，我们利用直到 t-1 时刻的信息。未知的 W_t 期望为0，而 W_{t-1} 可以用上一期的预测误差 (X_{t-1} - X_{t-1的预测值}) 来估计。经过代数变换，最终的预测公式会呈现为：
X_t的预测值 = α * X_{t-1} + (1 - α) * X_{t-1的预测值
这正是简单指数平滑的递归形式，其中 α = 1 + θ（注意 θ 通常为负值）。

类似地，双指数平滑（霍尔特线性法）与 ARIMA(0,2,2) 模型存在联系。虽然公式推导更复杂，但核心思想一致：这些指数平滑方法都可以在ARIMA模型的框架下找到对应的形式。

R语言实践：霍尔特-温特斯法 🛠️

理论之后，我们来看看如何在R中应用这些方法。我们将使用阿尔伯特省2020年2月至8月的新冠病例数数据。

R的基础stats包中提供了HoltWinters()函数。该函数可以实现加法或乘法形式的霍尔特-温特斯法，后者适用于具有乘法性季节性的数据。

对于我们的数据（无明显季节性），我们使用加法模型并禁用季节性成分。

以下是关键步骤：

拟合模型：如果不指定参数alpha（水平平滑）和beta（趋势平滑），函数会自动优化它们以最小化一步预测误差。
```
# 拟合霍尔特-温特斯模型，无季节性
hw_model <- HoltWinters(ab_cases_ts, gamma = FALSE)
```
查看参数：运行后可以查看优化得到的alpha和beta值。
进行预测与绘图：使用plot()函数可以绘制原始序列及拟合的平滑序列。
```
plot(hw_model)
```
手动调整参数：你可以尝试不同的alpha和beta值，观察平滑效果的变化。例如，alpha=1会使预测几乎等于上一期观测值（类似随机游走），而alpha=0.1会产生一个非常平滑但可能滞后于实际变化的序列。

扩展：ETS模型框架 📦

除了HoltWinters()，forecast包中的ets()函数提供了更现代、更灵活的指数平滑模型框架。

ets()函数允许通过模型代码字符串来指定误差类型（加法A/乘法M）、趋势类型（无N/加法A/乘法M）和季节性类型。例如，“AAN”表示具有加法误差、加法趋势且无季节性的模型。

其优势包括：

自动模型选择功能。
提供AIC等模型选择准则。
输出分解后的水平、趋势成分图，便于理解。

使用示例：

library(forecast)
ets_model <- ets(ab_cases_ts)
plot(ets_model) # 绘制序列分解图

方法对比：指数平滑 vs. ARIMA ⚖️

最后，一个自然的问题是：哪种方法更好？我们可以通过预测精度指标来比较。

我们同时对阿尔伯特省新冠病例数据拟合自动选择的ARIMA模型（使用auto.arima()）和ETS模型，并计算它们的预测误差。

代码示例：

# 拟合ARIMA模型
arima_model <- auto.arima(ab_cases_ts)
# 拟合ETS模型
ets_model <- ets(ab_cases_ts)

# 比较预测精度
accuracy(arima_model)
accuracy(ets_model)

在本案例中，auto.arima()选择的ARIMA(1,1,1)模型在均方根误差（RMSE）和平均绝对误差（MAE）上略优于ETS模型。对于测试数量数据，ARIMA(1,1,1)也表现稍好。

这提醒我们，没有放之四海而皆准的“最佳”模型。最佳实践是：

根据数据特性（趋势、季节性等）尝试多种模型。
不要完全依赖算法的自动选择，应结合统计检验和业务理解进行判断。
使用样本外预测精度作为重要的评估标准。

正如统计学家乔治·博克斯所言：“所有的模型都是错的，但有些是有用的。” 关键在于为你的具体问题找到那个“有用”的模型。

总结 🎯

本节课我们一起学习了时间序列分析中的指数平滑方法。

我们从简单指数平滑开始，理解了其作为加权移动平均平滑器的基本原理。
接着，我们介绍了双指数平滑（霍尔特线性法），它通过引入趋势方程来处理具有趋势的数据。
我们探讨了这些平滑方法与ARIMA模型（如ARIMA(0,1,1)）之间的深刻联系。
最后，我们在R中实践了HoltWinters()函数和更先进的ets()函数，并使用阿尔伯特省新冠疫情数据比较了指数平滑与ARIMA模型的预测性能。

指数平滑法提供了不同于经典ARIMA模型的建模视角，是时间序列预测工具箱中的一个重要组成部分。理解其原理和适用场景，将帮助你在面对实际数据时拥有更丰富的分析手段。

024：GARCH过程

在本节课中，我们将学习一个特殊的主题：GARCH模型。GARCH是“广义自回归条件异方差”的缩写。这个模型扩展了我们之前学习的自回归和ARMA模型，但这次我们将重点关注方差的变化。在本课程中，我们一直假设白噪声过程的方差是恒定的，但在现实中，时间序列的方差（或波动性）常常会随时间变化，这在金融时间序列中尤为常见。本节课将介绍如何对这种时变方差进行建模。

模型背景与动机

上一节我们回顾了ARMA等模型，它们主要关注序列的均值结构。本节中，我们来看看如何将注意力转向方差。

在标准的时间序列模型中，我们通常假设误差项（白噪声）的方差是恒定的，即同方差性。然而，在许多实际应用中，特别是金融领域，序列的波动性会随时间发生显著变化。例如，股票价格可能在某个时期平稳波动，而在另一个时期剧烈震荡。GARCH模型就是为了捕捉和建模这种“条件异方差性”而设计的。

从ARCH模型开始

为了理解GARCH，我们首先从更简单的ARCH模型入手。ARCH代表“自回归条件异方差”。

ARCH(1)模型定义

ARCH(1)模型由两个方程定义。第一个方程描述观测到的“收益率”或“增长率”序列 ( R_t )：

[
R_t = \sigma_t W_t
]

其中，( W_t ) 是独立同分布的标准正态白噪声，即 ( W_t \sim N(0, 1) )。关键点在于，标准差 ( \sigma_t ) 是随时间变化的。

第二个方程定义了时变方差 ( \sigma_t^2 ) 如何依赖于过去的收益率：

[
\sigma_t^2 = \alpha_0 + \alpha_1 R_{t-1}^2
]

这里，( \alpha_0 > 0 ) 且 ( \alpha_1 \geq 0 ) 是模型参数。( \alpha_0 ) 可以看作是方差的“基线”或下界。( \alpha_1 ) 衡量了前一期收益率的平方对当前方差的影响程度。

ARCH(1)模型的性质

以下是ARCH(1)模型的一些关键性质：

条件分布：给定过去的信息（如 ( R_{t-1} )），( R_t ) 的条件分布是均值为0、方差为 ( \sigma_t^2 ) 的正态分布。这正是“条件异方差”的含义——方差依赖于过去的值。
无条件性质：虽然条件方差在变化，但收益率序列 ( R_t ) 本身的无条件均值仍为0，并且序列是不相关的。然而，这并不意味着 ( R_t ) 是独立的，因为方差之间存在依赖关系。
平方序列：如果我们考虑平方收益率 ( R_t^2 )，可以证明它服从一个AR(1)过程：
[
R_t^2 = \alpha_0 + \alpha_1 R_{t-1}^2 + u_t
]
其中 ( u_t ) 是一个均值为0的噪声项（但不再是正态分布）。这为我们分析方差动态提供了另一个视角。
参数约束：为了保证过程平稳且方差有限，需要 ( 0 \leq \alpha_1 < 1 )。如果进一步要求 ( R_t^2 ) 的四阶矩有限（即方差有限），则需要更强的约束 ( 0 \leq \alpha_1 < 1/\sqrt{3} )。违反此约束会导致分布的尾部比正态分布更“厚”。

扩展到GARCH模型

ARCH模型的一个自然扩展是GARCH模型，它在方差方程中加入了自回归项，使其更灵活、通常也更简洁。

GARCH(1,1)模型定义

GARCH(1,1)模型是实践中最常用的形式。其方程如下：

收益率方程：
[
R_t = \sigma_t W_t
]

方差方程：
[
\sigma_t^2 = \alpha_0 + \alpha_1 R_{t-1}^2 + \beta_1 \sigma_{t-1}^2
]

与ARCH(1)相比，方差方程中多了一项 ( \beta_1 \sigma_{t-1}^2 )。这意味着当前方差不仅受上一期收益率冲击（( R_{t-1}^2 )）的影响，还受上一期方差（( \sigma_{t-1}^2 )）本身的影响。

为何称为“广义”自回归？

类似于ARCH模型，GARCH模型的平方序列 ( R_t^2 ) 可以表示为一个ARMA过程。具体来说，GARCH(1,1)对应的 ( R_t^2 ) 是一个ARMA(1,1)过程。参数 ( \alpha_1 ) 扮演了“自回归”部分的角色，而参数 ( \beta_1 ) 扮演了“移动平均”部分的角色。这就是“广义自回归”中“自回归”一词的由来。

模型估计与应用

理论介绍完毕，接下来我们看看如何在实践中应用GARCH模型。

模型拟合与R语言实现

在R语言中，我们可以使用 fGarch 包中的 garchFit() 函数来拟合GARCH模型。该函数默认使用最大似然估计法。以下是基本步骤：

数据准备：对于价格等原始序列 ( X_t )，我们通常先计算其对数收益率作为建模对象 ( R_t )：
[
R_t = \log(X_t) - \log(X_{t-1})
]
这近似于百分比收益率，并且通常能获得更稳定的序列。
模型设定与拟合：我们可以指定均值方程和方差方程。例如，拟合一个包含AR(1)均值项和GARCH(1,1)方差项的模型：
```
# 假设 r 是对数收益率序列
library(fGarch)
model <- garchFit(formula = ~ arma(1,0) + garch(1,1), data = r)
summary(model)
```
结果解读：summary() 输出会提供参数估计值、标准误、显著性检验，以及一系列关于标准化残差的诊断检验（如Ljung-Box检验、ARCH-LM检验等），用于评估模型拟合效果。信息准则（AIC, BIC）可用于模型比较。

实例：COVID-19病例数据的波动性分析

我们以阿尔伯塔省的COVID-19每日新增病例数据为例。原始序列显示，病例数在不同波次期间波动性明显不同。我们对其取对数并计算一阶差分（近似增长率），然后尝试用GARCH类模型进行拟合。

通过尝试不同的模型（如纯GARCH(1,1)、AR(1)-GARCH(1,1)、MA(1)-GARCH(1,1)），并比较信息准则和残差诊断图，我们可以选择一个相对合适的模型来描述病例增长率的波动性模式。拟合好的模型可以给出条件标准差的估计序列，直观展示波动性如何随时间演变。

模型扩展

GARCH模型家族非常庞大，有许多变体以适应不同的数据特征：

高阶模型：可以定义ARCH(p)或GARCH(p, q)模型，让方差依赖于更多期的过去信息。
集成GARCH (IGARCH)：用于建模方差中的持久性冲击。
指数GARCH (EGARCH)：允许正面和负面冲击对波动性有不对称的影响。
GARCH-M模型：在均值方程中引入条件方差或条件标准差项，用于研究风险与收益的关系。
与其他模型结合：GARCH结构可以作为误差项，与任何描述序列均值结构的模型（如ARIMA、回归模型）结合使用，形成完整的均值-方差模型。

总结

本节课中我们一起学习了GARCH过程。我们从恒定方差的假设局限出发，引入了能够刻画时变波动性的ARCH模型。通过分析ARCH(1)模型的条件与无条件性质，我们理解了其核心思想。接着，我们将其推广到更强大、更简洁的GARCH(1,1)模型，并解释了其名称的由来。最后，我们介绍了如何在R语言中实现GARCH模型的拟合与诊断，并以COVID-19数据为例进行了演示。GARCH模型是金融计量经济学和时间序列分析中不可或缺的工具，它为理解和预测波动性提供了坚实的框架。

025：STAT479 课程总览

在本节课中，我们将对 STAT479 时间序列分析课程的全部内容进行一次全面的回顾。我们将从课程起点开始，梳理我们学习过的所有核心概念、模型和方法，并总结它们之间的联系。

课程概述与起点

欢迎来到 STAT479 时间序列分析的最后一讲。祝贺大家完成了这门课程。今天我们将进行课程回顾，总结我们从哪里开始，以及经过超过30小时的讲座视频后，我们最终到达了何处。

时间序列是统计学中一个非常精妙的领域。它融合了许多不同的内容。一方面，它有点像回归分析，但回归分析假设所有观测值相互独立，而时间序列则明确知道观测值之间存在时间相关性，我们必须处理这些相关性。另一方面，我们试图拟合模型来观察数据，但很多时候，我们的模型只是在尝试对过程的“噪声”进行建模。对于一个平稳过程，平稳的噪声过程很像纯粹的噪声，它是一系列随时间的随机波动，遵循某种平稳分布。然而，根据平稳过程的行为，我们可以了解很多关于时间序列的信息，包括如何理解、估计、建模和预测它，这些都是关键部分。

我们是从观察不同类型的“噪声”开始这门课程的。那时，我们介绍了白噪声、自回归过程和移动平均过程。我们最终得到了白噪声、AR过程和MA过程，这些在课程一开始就介绍了，后来你会发现，这基本上就是我们整个课程要讨论的核心，因为它们是三个主要的构建模块。

我们在课程开始时引入的另一个主要概念是自协方差函数。自协方差是本课程的关键部分，因为课程中几乎所有内容都源于自协方差函数。无论是绘制ACF图、估计模型参数，还是估计谱密度，都离不开自协方差。

我们还讨论了互协方差和平稳性的概念。在实践中，我们经常有两个时间序列，因此可能对互协方差感兴趣，即两个序列如何相互关联或相关。我们花了大量时间讨论平稳性。

平稳性是本课程的另一个关键点，因为我们需要它来进行统计分析。只有当我们要估计的对象不随时间变化时，我们才能进行统计。如果均值随时间变化，我们就无法计算样本均值。自协方差也是如此，我们只能从数据中估计自协方差或查看ACF图，前提是我们的过程是平稳的。因为最终为了估计这些量，我们需要进行某种平均，而只有当我们试图估计的量在不同时间点上保持不变时，跨时间点平均才有意义。

第一章：基础估计与检验

在第一章中，我们讨论了均值估计、方差估计和自协方差估计。我们还研究了ACF图和白噪声检测。在实践中，我们无法先验地知道噪声是白噪声，还是来自某个移动平均过程或AR模型。从本课程中可能学到的一点是，如果你绘制出一个平稳过程，它可能看起来就像一堆静态噪声，但其中可能隐藏着有趣的特征，也可能没有。如果是白噪声，从某种意义上说，它就不那么有趣，因为它已经尽可能地随机了。

第二章：统计模型与平滑

在第二章中，我们开始研究时间序列的统计模型。我们首先看的是线性回归。正如之前提到的，尝试用线性回归来拟合数据中的线性趋势是可以的，但必须意识到，最小二乘估计器并非为处理相关观测值而设计。这里的关键点是，最小二乘法假设观测值不相关。因此，在处理时间序列数据应用线性回归时，必须记住这一点。你可能会得到与尝试用带漂移项的时间序列模型拟合相同的结果，但在线性模型中，我们假设误差或观测值之间不相关，而这在有自回归或移动平均过程的情况下显然不成立。

第二章中还有一个额外的部分，即平滑的概念。平滑是一个很好的工具，但它更像是一个定性主题，因为没有绝对“正确”的平滑方法或平滑量。它更像是一种可以用来发现时间序列数据中趋势的技术。使用平滑时需要小心，因为你通常不希望平滑数据后再进行统计分析，因为这会改变数据的概率特性。平滑更像是在最后阶段，你可以观察平滑后的结果，看看趋势是上升还是下降。它是一个非常有用的工具，我们讨论了许多不同的方法，如局部加权散点图平滑法、三次样条、核平滑器，以及简单的移动平均平滑器，后者使用移动平均过程来过滤时间序列。

平滑以这种方式起作用，正如我们在谱分析部分提到的，它就像一个低通滤波器，去除了数据中的高频或快速波动，只留下平滑的低频趋势。

第三章：ARIMA模型

之后，我们开始讨论ARIMA模型，这是我们课程大部分时间讨论的内容。ARIMA模型，即自回归积分移动平均模型。在这一部分，我们讨论了自回归部分、移动平均部分，并将它们组合起来引入了ARIMA模型，讨论了其各个方面。一个需要记住的关键点是，可能有不止一个模型可以拟合你的数据。因为最终你需要担心参数冗余、因果性或非因果性、可逆性或不可逆性等问题。我们发现，在随机意义上，可能存在不同类型但行为完全相同、无法区分的时间序列过程。因此，在使用这些模型时，我们必须非常小心，这就是为什么我们通常强制实施因果性和可逆性等概念。我们还希望确保AR部分和MA部分没有公因子，因为如果有公因子，在查看特征多项式时会导致问题。

在这一部分，我们还讨论了平稳性检验和自相关检验。本课程中我们感兴趣的主要问题包括：我的数据是平稳的吗？或者我的数据是不相关的吗？对于数据不相关的问题，我们回顾了Box-Pierce和Ljung-Box检验，它们都用于检验数据中是否存在显著的自相关。我们还研究了Durbin-Watson检验和Breusch-Godfrey检验，这些检验用于查看线性模型的残差。许多这些工具不仅适用于时间序列，也适用于一般的线性模型，这就是为什么你可以在R的lmtest包中找到它们，因为你是在检验线性模型的残差，想知道残差是否看起来不相关。如果残差看起来相关，那么可能需要做更多工作来找出这些相关性的来源或如何建模这些相关性。

我们还讨论了平稳性检验，本课程主要讨论的是增广迪基-富勒检验。它是众多所谓“单位根检验”中的一种。我们之前讨论过，如果我们有非平稳数据或非平稳ARMA过程，那么自回归多项式的根会落在复平面的单位圆上，这将给我们一个非平稳过程。因此，当你进行单位根检验时，你的零假设是认为根在单位圆上，备择假设是根不在单位圆上。所以零假设是数据非平稳，备择假设是数据平稳，即所有根都离单位圆足够远，我们可以安全地声称拥有一个平稳时间序列过程。这有点令人困惑，因为人们总倾向于认为零假设是平稳性，但实际上是非平稳性，这是在自己数据上运行此类检验时需要记住的一点。

还有Phillips-Perron检验，这是另一种基本做相同事情、检验平稳性的方法，类似于增广迪基-富勒检验。Phillips-Perron检验使用略有不同的方式计算检验统计量，但基本思想相同，都是问：我的特征多项式是否至少有一个根在单位圆上？如果没有，我们就拒绝零假设，并声称我们有一个平稳过程。

另一个需要记住的是，在进行这些检验时，必须注意涉及的参数，如自由度、滞后阶数等。能够检验这些不同方面总是好的，因为可能在一个设定下看不到显著结果，但调整检验中的一些参数后，就可能发现显著结果。

我们还讨论了ACF和PACF，因为它们是时间序列中非常有用的工具，可以直观地理解数据发生了什么。ACF就是自相关函数，它是对称的，具有许多良好的数学性质。绘制数据的ACF图也很有趣，它可以告诉我们很多关于显著自相关的信息，这些显著自相关至少能提示我们数据中可能存在自回归或移动平均过程，或两者兼有。我们还可以使用互补工具——偏自相关函数来分析数据。偏自相关考察的是偏相关的概念，并将其扩展到自相关。偏相关基本上是说，在给定另一个变量的条件下，两个变量不相关。我的标准例子是房屋价格和汽车价格，它们可能相关，但如果在考虑个人收入后，它们可能就不相关了。如果有人更有钱，他可能在房子和汽车上都花更多钱，反之亦然。因此，在这种情况下，以个人收入为条件，汽车和房屋价格可能看起来彼此不相关。偏自相关函数基本上是说，如果我有一个时间序列数据，我想知道某些时间点在给定所有中间点的情况下是否相关。这为我们提供了一种查看时间序列数据中相关性的新方法。

我们发现了一个很好的小图表，基本上说明了类似这样的内容：对于AR过程，ACF呈几何衰减，而PACF在滞后p阶后消失。当然，对于数据来说，它永远不会精确为零，我们仍然会在零附近有一些噪声或小的波动，但从数学上讲，在滞后p阶后它等于零。对于移动平均过程，我们看到相反的情况：ACF在滞后q阶后消失，而PACF呈几何衰减。因此，这两个工具对于真正理解给定时间序列数据集的情况非常有用。再次强调，需要是平稳的才能查看这些图。如果你的时间序列不平稳，你就无法正确进行这种估计，你可能会在ACF和PACF中看到到处都是极大的峰值。从某种意义上说，你可以用它来确定是否拥有非平稳时间序列，因为对于非平稳序列，你通常会看到到处都是巨大的相关性，通常是巨大的正相关性。

第四章：估计与预测

接下来我们在本课程中讨论的大事是估计和预测，这是处理时间序列时想要做的两件大事。我们有时序数据，当然可以观察并试图理清发生了什么，但最终我们想要估计模型的参数，并且我们想要进行预测，因为几乎可以肯定，如果你在观察某物随时间演变，你不仅关心它从哪里开始、如何到达这里，还关心它未来会走向何方。我几乎无法想象任何时间序列场景下你会对预测不感兴趣，无论是金融（我想知道我的投资是否会赚钱）、环境（我们几年后是否会生活在沙漠中，而不是这个宜人的埃德蒙顿？），还是COVID病例（我是否需要继续待在后院的房子里躲避所有人，或者一切都会好转？）。每当我们观察随时间变化的事物时，我们都想进行预测，而在尝试进行预测时会出现一些有趣的复杂性。

我们在这一部分研究了对AR过程的预测，然后更一般地对ARMA过程以及季节性模型进行了预测。回想一下，对于AR过程，我们有我们的朋友——尤尔-沃克方程。一般来说，AR过程的估计和预测是时间序列中最简单的设置，因为尤尔-沃克方程非常好用，易于计算，你可以手动完成，就像你在一次作业中必须做的那样。你可以估计参数，可以使用AIC等信息准则来确定模型中应包含多少项，最终你可以进行预测。关于AR过程预测的一个好处是，对于预测，我们只需要前p个观测值，假设我们正在处理一个AR(p)过程。当然，我们仍然需要整个数据集，因为我们将使用整个数据集来估计模型参数φ1到φp，所以更多的数据仍然有好处。但好处是，我们不必使用整个数据集来尝试估计或预测未来的时间点，这使得运行起来非常容易。

这与对ARMA过程做同样的事情形成对比。因为如果我们有一个具有非平凡移动平均分量的ARMA过程，我们就无法使用尤尔-沃克方程这种好用的设置，而是最终使用最大似然估计。同样，我们有两种不同的方法：条件最大似然和非条件最大似然。区别在于，当你使用条件方法时，你会得到一个漂亮的线性方程组；如果不使用条件方法，你会得到一个非线性系统。如今计算机可以处理非线性系统，但在过去，你真的希望一切都是线性方程组，因为如果是线性的，你可以快速求解；如果是非线性的，你就会陷入困境。如今，非线性问题对计算机来说仍然很难，但可以做到，而且经常被做到，比如在R中拟合ARMA或ARIMA模型的包中。

另一个问题是预测。因为预测现在需要所有数据。当我们试图处理这些移动平均部分时，我们需要数据在某种意义上追溯到尽可能远。我们曾有一个有点滑稽的设定，即想象我们拥有无限的数据，这在数学上非常好，但在实践中行不通。不过至少在数学上，我们可以考虑拥有无限的数据，数据追溯到无限的过去，然后我们可以考虑基于此进行建模。在实践中我们没有无限数据，这导致了几种我们可以用于预测的方法。

我们有Durbin-Levinson算法。这基本上是如何使用我们的自协方差矩阵来求解方程组。它在数值线性代数以及数值ODE等领域也有应用，时间序列中也是如此。我们还研究了创新算法。关于创新算法的关键一点是，这些基本上是一步超前预测的残差，它们的行为很像白噪声。所以它们有点像白噪声的估计。解释它的最佳方式是，如果我们观察我们认为的第t个观测值与实际看到的第t个观测值之间的差异，剩下的就是残差，它很像白噪声过程。其思想是，如果我们试图估计这些对象，相对于这些创新（即残差）进行预测，在某些情况下通常比尝试Durbin-Levinson等方法更容易。同样，这两种方法在方程表达上都非常复杂和繁琐。

第五章：季节性模型与谱方法

我们还有季节性ARIMA或SARIMA模型。从某种意义上说，季节性模型并没有真正太新的东西，它实际上只是说我们有另一个AR和/或MA部分，可能还有一个季节性差分，我们只每隔一定时间点应用它。这在数据中非常常见，因为通常你观察的数据具有某种周期性，无论是温度的年度周期、金融的季度周期，还是五月份这里日照的每日周期。你的时间序列中经常有这些周期性行为，寻找季节性模型很重要，但最终季节性模型与非季节性模型并没有本质上的深刻不同。不过拟合它们确实需要一些努力，你必须考虑什么是合理的周期。你可以像auto.arima函数那样从数据中尝试估计它，或者你可以根据问题的背景尝试做出有根据的猜测。如果你有每日观测值，每周周期很常见；如果你有月度数据，每年12个月可能有一个周期。在尝试查看真实数据、真实时间序列数据时，总是需要注意这些事情。

之后，我们进入了频域，事情变得更有趣了。在频域中，我们可以做很多以前尝试做的事情，但我们可以从频率的角度而不是时间的角度来做。我们主要研究的是周期图。记住，周期图是离散傅里叶变换应用于我们时间序列后的平方幅度。离散傅里叶变换就像本课程中的许多东西一样，只是一个巨大的线性方程组，你通常不想直接求解。所以你做什么？你使用快速傅里叶变换。在频率设置中，我们试图弄清楚的关键是周期图。原因之一是它会告诉我们时间序列数据中存在的不同频率的强度，我们可以用它来寻找周期性行为。我们非常简要地讨论了周期过程，这是一种与我们在本课程中花费大量时间研究的ARMA过程略有不同的平稳随机过程。这些过程具有非常强的频率，我们可能试图通过周期图来检测。

但周期图也是估计谱密度的一种方式。谱密度是另一个关键点，以至于我们花了三讲不同的课来讨论估计谱密度的不同方法。有一种非参数方法，它基于周期图。还有一种参数方法，它基于AR过程。同样，有两种不同的方法来估计谱密度。我们可以直接计算周期图并使用它，但这通常不是好的方法。我们可以做的是尝试不同的技术来改进我们的估计，比如分组、Bartlett方法（将数据分成几块，对每块拟合周期图，然后平均周期图）。记住，我们假设有一个平稳过程，所以数据前半部分的周期图应该与后半部分的周期图相同。这是平稳性成为贯穿整个课程关键主题的另一个原因，如果我们要使用任何这些技术，我们确实需要一个平稳过程。

你可以尝试将多个周期图平均在一起，可以使用分组（即将多个频率平均在一起），可以使用平滑（即使用Daniel核方法），这只是应用不同类型权重的另一种方式，是一种更广义的分组方式，以尽可能平滑周期图作为谱密度的非参数估计量。

你也可以采用参数方法，这是人们也喜欢使用的一种。我不是一个应用型的、处理大量时间序列数据的工程师类人物，所以我不知道不同领域（金融数据、信号处理、环境数据）的人们在实践中喜欢使用什么，他们可能有自己的偏好。但我能说的是，有不同的方法，有些人喜欢基于AR过程的非参数方法，有些人喜欢参数方法。在那种情况下，你只需将AR过程拟合到你的数据，R中的函数可能使用尤尔-沃克方程，我认为它只是使用AR函数拟合那些参数，使用AIC，然后直接从该过程计算谱密度，因为正如我们之前讨论的，我们可以过滤时间序列，并确定这将如何影响谱密度。

过滤的概念。如果我们有一些Y_t，它等于从负无穷到正无穷对a_i * X_{t-i}的求和，现在我们有了基于某些系数a的时间序列X的过滤版本。我们发现，对于这种过滤将如何影响谱密度，有一个简单的公式。如果我们有白噪声，谱密度是常数，每个频率都与其他频率同等表示。但我们可以使用白噪声转变为自回归和移动平均以及更一般的ARMA过程的方式来过滤它。如果我们这样做，那么我们可以快速找出一般ARMA过程的谱密度，正如我们在谱方法的一讲中所做的那样。

最终，所有这些谱方法实际上只是为我们提供了一种不同的方式来查看时间序列数据，通过查看频谱（即周期图或谱密度的某种估计），我们可以了解时间序列数据中哪些频率最突出，这又会给我们一些关于如何思考这些数据的思路。它可以为拟合季节性模型找到有用的周期，也可以发现与白噪声的偏差，因为如果我们有一个白噪声过程，我们应该看到一个平坦的周期图或谱密度；如果我们没有看到，我们就有一个相当强的迹象表明我们没有白噪声。

第六章：扩展主题与总结

课程最后还有一些扩展主题，它们本身也相当精妙。我们有指数平滑的概念，它为我们提供了另一种拟合时间序列模型和进行预测等工作的工具。这是一种与我们之前所做的略有不同的方法，但我们确实看到它与ARMA过程密切相关。我们还研究了GARCH过程。到目前为止，在本课程的所有其他部分，我们假设背后有白噪声驱动一切，总是有白噪声，并且白噪声具有恒定方差。这在实践中不一定发生，事实上，方差通常会随时间变化。GARCH过程试图考虑到这一点，这是金融人士喜欢它的原因之一，因为波动性（无论你持有股票还是债券）就像随机游走一样，没有太多可做的，但你仍然可以对波动性建模：我们是在小幅上下波动，还是在大幅上下波动？这在尝试基于历史时间序列数据做出决策时会产生很大差异。

以上就是对本课程所做内容的回顾。这是一个非常大的主题，时间序列是这样一个领域，有大量的研究投入其中，因为每当你在统计学中提出一个有趣的东西时，你几乎可以通过说“如果我们加上时间维度，然后在时间上再做一遍会怎样？”来把它变成另一个完整的研究计划或博士论文。这基本上就是时间序列，但事情在涉及时间相关性时变得相当复杂。经典统计学中有太多内容基于不相关的独立数据，而时间序列为我们提供了一些工具，用于思考当数据可能存在相关性时的情况。是的，整个课程基本上就是ARIMA模型和自相关，我们所做的一切中大约90%就是这两件事，但其中涉及如此多的细微之处：如何与之合作、如何估计、如何建模、如何思考、如何处理真实数据。这是一个非常有趣的探索性主题，我希望你在这门课程中学到了很多，希望你能走出去，尝试自己处理更多时间序列。也许有人会付钱让你去为某个公司做时间序列分析，因为时间序列数据无处不在，知道如何处理它非常有用。

再次感谢大家参加所有这些讲座。如果你没有参加所有讲座，请回去看讲座，然后再回到这里，我会感谢你参加了所有讲座。对着我的相机讲话，想象着一大群明亮、快乐的面孔一直对我微笑，这绝对是一种体验。希望你在电脑前盯着屏幕时就是这样做的。如果没有，嗯，我实际上没有任何数据表明情况相反，所以我对自己的感觉相当好。无论如何，再次感谢，祝你在未来查看所有这些数据以及在实际世界中使用所学知识时好运。

posted @ 2026-03-26 12:21 布客飞龙III 阅读(55) 评论(0) 收藏举报

刷新页面返回顶部

龙哥盟

阿尔伯塔时间序列分析笔记-全-

阿尔伯塔时间序列分析笔记（全）

001：噪声过程

概述：什么是时间序列分析？

白噪声：时间序列的基石

自回归过程：用过去预测现在

移动平均过程：平滑噪声

重要概念与定义

线性过程：一个统一的框架

示例：将AR(1)过程表示为线性过程

总结

002：自协方差与平稳性

自协方差与自相关

示例：带漂移项的 AR(1) 过程

平稳性

联合平稳性与样本均值的方差

总结

003：估计与线性模型

回顾平稳过程

估计自协方差与自相关函数

在R中分析自相关函数

检测白噪声

时间序列的统计模型

时间序列模型 vs. 线性回归模型

R语言实战：分析处方药价格数据

004：平滑方法

概述

移动平均平滑器

核平滑器

局部加权散点图平滑（LOWESS）

三次平滑样条

R语言应用演示

重要注意事项：平滑的用途与局限

总结

005：AR与MA模型理论

概述

自回归模型

因果性

自回归算子的根

移动平均模型

可逆性

总结

006：ARMA与ARIMA模型 🧮

回顾：AR与MA过程 📚

自回归（AR）过程

移动平均（MA）过程

结合：ARMA模型 ⚙️

模型的可辨识性与简化

因果性与可逆性

扩展：ARIMA模型 📈

差分算子

ARIMA模型定义

包含常数项的ARIMA

总结 🎯

007：自相关与平稳性

检验自相关性

Box-Pierce 与 Ljung-Box 检验

Durbin-Watson 检验

Breusch-Godfrey 检验

检验平稳性

(Augmented) Dickey-Fuller 检验

Phillips-Perron 检验

总结

008：在R中检验平稳性与自相关函数

准备工作：生成模拟数据

模型拟合简介

自相关检验

Box-Pierce 与 Ljung-Box 检验

Durbin-Watson 检验

Breusch-Godfrey 检验

平稳性检验（单位根检验）

Augmented Dickey-Fuller (ADF) 检验

Phillips-Perron (PP) 检验

自相关函数与偏自相关函数预览

总结

009：自相关与偏自相关函数（第一部分）

从相关到偏相关

什么是偏相关？

如何计算偏相关？