数学建模实战：基于ARMA与ARIMA模型的时间序列预测

柴木头 B2B电商

1370人浏览 · 2025-09-26 13:52:37

柴木头 B2B电商 · 2025-09-26 13:52:37 发布

简介：在数学建模中，ARMA（自回归移动平均）和ARIMA（自回归整合移动平均）模型是分析和预测时间序列数据的核心工具，尤其适用于具有周期性与趋势性特征的数据，如生物种群数量变化。ARMA模型结合自回归与移动平均机制，有效捕捉数据的动态依赖关系；ARIMA模型通过差分处理非平稳序列，提升建模精度。本项目以MATLAB为实现平台，利用 arima 、 estimate 、 autocorr 等函数完成数据预处理、阶数识别、模型构建与性能评估，通过AIC/BIC准则和残差分析优化模型效果，助力生态预测与科学决策。

1. ARMA模型基本原理与适用场景

1.1 ARMA模型的基本构成与数学表达

自回归移动平均（ARMA）模型是时间序列分析中的核心工具，适用于平稳序列的建模。其结构由两部分组成：自回归（AR）部分描述当前值与历史值的线性关系，数学形式为 $ X_t = c + \sum_{i=1}^p \phi_i X_{t-i} + \varepsilon_t $；移动平均（MA）部分刻画当前值对过去误差的响应，表示为 $ X_t = \mu + \varepsilon_t + \sum_{j=1}^q \theta_j \varepsilon_{t-j} $。二者结合形成ARMA(p, q)模型：

X_t = c + \sum_{i=1}^p \phi_i X_{t-i} + \varepsilon_t + \sum_{j=1}^q \theta_j \varepsilon_{t-j}

其中 $\varepsilon_t$ 为白噪声过程。

1.2 模型适用条件与现实应用场景

ARMA模型要求时间序列具有 弱平稳性 ——均值、方差恒定且自协方差仅依赖于滞后阶数。典型应用场景包括金融收益率预测、气候变量短期建模及工业传感器数据滤波等。对于非平稳序列，需通过差分预处理转化为ARIMA框架。该模型在MATLAB中可通过 arima(p,0,q) 直接定义，便于后续参数估计与诊断分析。

2. ARIMA模型结构与非平稳序列处理方法

时间序列分析在现代数据科学中占据着至关重要的地位，尤其在金融、生态学、气象学和工业控制等领域，大量观测数据呈现出明显的动态演化特征。其中，自回归积分移动平均模型（Autoregressive Integrated Moving Average, 简称 ARIMA）作为经典的时间序列建模工具，因其能够有效处理非平稳序列而被广泛采用。与传统的ARMA模型仅适用于弱平稳序列不同，ARIMA通过引入差分机制，实现了对趋势性、随机漂移等非平稳成分的系统化解析。该模型的核心思想在于“先使序列平稳，再进行建模”，从而将复杂的现实问题转化为可操作的统计推断任务。

ARIMA模型的形式化表达为ARIMA(p, d, q)，其中三个参数分别代表自回归阶数 $ p $、差分阶数 $ d $ 和移动平均阶数 $ q $。这一结构不仅扩展了ARMA模型的应用边界，也带来了更为精细的数据预处理要求和建模逻辑设计。理解ARIMA的理论构成，是掌握其应用能力的前提；识别非平稳性的来源并合理选择差分阶数，则直接决定模型的有效性；最终，整个建模流程必须遵循从数据诊断到模型评估的闭环逻辑，确保预测结果具备统计稳健性和实际解释力。

本章将深入剖析ARIMA模型的内在结构，系统阐述其三大组成部分——自回归（AR）、差分整合（I）与移动平均（MA）的数学基础与作用机制，并结合典型非平稳序列的识别方法，探讨如何通过单位根检验与差分技术实现序列平稳化。在此基础上，构建完整的ARIMA建模流程框架，涵盖模型识别、参数估计与诊断检验的关键环节，并通过跨领域案例对比揭示其普适性与局限性。

2.1 ARIMA模型的理论构成

ARIMA模型的本质是对原始时间序列进行差分变换后，在新生成的平稳序列上建立ARMA模型。因此，其理论构成可以分解为三个相互关联的部分：自回归（AR）部分描述当前值与历史值之间的线性依赖关系；移动平均（MA）部分刻画当前观测值与过去扰动项之间的动态反馈；而差分整合（I）则负责消除原始序列中的趋势性或随机游走成分，使其满足平稳性假设。这三个模块共同构成了ARIMA模型的完整表达体系。

为了更清晰地展现这一结构，考虑一个一般形式的ARIMA(p, d, q)过程：

\phi(B)(1 - B)^d X_t = \theta(B)\varepsilon_t

其中：
- $ X_t $ 是时间序列在时刻 $ t $ 的观测值；
- $ B $ 是滞后算子（Backshift Operator），满足 $ BX_t = X_{t-1} $；
- $ (1 - B)^d $ 表示 $ d $ 阶差分操作；
- $ \phi(B) = 1 - \phi_1 B - \phi_2 B^2 - \cdots - \phi_p B^p $ 是自回归多项式；
- $ \theta(B) = 1 + \theta_1 B + \theta_2 B^2 + \cdots + \theta_q B^q $ 是移动平均多项式；
- $ \varepsilon_t \sim WN(0, \sigma^2) $ 是白噪声误差项。

该公式表明，ARIMA模型首先对原序列进行 $ d $ 次差分得到平稳序列 $ Y_t = (1 - B)^d X_t $，然后在 $ Y_t $ 上拟合一个ARMA(p, q)模型。这种“差分+ARMA”的两步策略，使得ARIMA能够灵活应对多种非平稳模式。

下面分别从三个子模块出发，逐层解析其数学表达与统计意义。

### 2.1.1 自回归（AR）部分的数学表达与意义

自回归模型的基本思想是利用变量自身的过去值来预测当前值。对于一个AR(p)过程，其定义如下：

Y_t = \phi_1 Y_{t-1} + \phi_2 Y_{t-2} + \cdots + \phi_p Y_{t-p} + \varepsilon_t

该模型假设当前时刻的值 $ Y_t $ 是前 $ p $ 个时刻值的线性组合加上一个独立同分布的误差项 $ \varepsilon_t $。系数 $ \phi_i $ 反映了第 $ i $ 阶滞后项对当前值的影响强度，其符号和大小决定了序列的记忆特性和动态响应速度。

例如，当 $ \phi_1 > 0 $ 且接近1时，序列表现出较强的正相关性，即高值之后倾向于继续出现高值，形成持续的趋势惯性；而若 $ |\phi_1| > 1 $，则会导致系统不稳定，序列发散，不符合平稳性要求。因此，AR模型的平稳性条件要求其特征方程 $ 1 - \phi_1 z - \phi_2 z^2 - \cdots - \phi_p z^p = 0 $ 的所有根都在单位圆外。

在ARIMA框架中，AR部分的作用在于捕捉序列内部的长期记忆结构。它能有效建模如经济增长惯性、种群增长延迟效应等具有持续影响的现象。值得注意的是，AR项并不直接作用于原始序列 $ X_t $，而是施加在差分后的平稳序列 $ Y_t $ 上，这意味着所建模的是“变化量”之间的依赖关系，而非原始水平值。

以下是一个简单的AR(2)模型实现代码示例（使用Python）：

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.arima_process import ArmaProcess

# 定义AR系数：Y_t = 0.6*Y_{t-1} - 0.3*Y_{t-2} + ε_t
ar_coefs = np.array([1, -0.6, 0.3])  # 注意：ArmaProcess要求首项为1，其余为负号
ma_coefs = np.array([1])             # MA部分为空

# 生成AR(2)过程
process = ArmaProcess(ar_coefs, ma_coefs)
y_ar2 = process.generate_sample(nsample=500)

# 绘图
plt.figure(figsize=(10, 4))
plt.plot(y_ar2)
plt.title("Simulated AR(2) Process")
plt.xlabel("Time")
plt.ylabel("Y_t")
plt.grid(True)
plt.show()

代码逻辑逐行解读：
- 第4行：定义AR系数数组， [1, -0.6, 0.3] 对应多项式 $ 1 - 0.6B + 0.3B^2 $，注意 statsmodels 库要求输入格式为带首1的滞后多项式。
- 第5行：MA系数设为 [1] ，表示无MA成分。
- 第8行：调用 ArmaProcess 创建过程对象， generate_sample 方法生成500个样本点。
- 后续绘图展示模拟序列的时间路径，可见其具有一定的波动周期性和衰减振荡特征。

此代码可用于教学演示或模型仿真测试，帮助理解AR过程的动力学行为。

### 2.1.2 移动平均（MA）部分的随机误差建模机制

移动平均模型关注的是当前观测值如何受到过去随机冲击的影响。一个MA(q)过程定义为：

Y_t = \varepsilon_t + \theta_1 \varepsilon_{t-1} + \theta_2 \varepsilon_{t-2} + \cdots + \theta_q \varepsilon_{t-q}

这里，$ \varepsilon_t $ 是白噪声序列，代表不可预测的外部扰动。MA模型的核心在于：当前值是由最近 $ q $ 个扰动项的加权和构成的。与AR模型不同，MA过程总是平稳的（只要参数有限），因为它本质上是一个滑动窗口对噪声的滤波。

MA项的意义在于解释“意外事件”的持续影响。例如，在生态系统中，某年极端气候导致生物死亡率上升，这一冲击可能不会立即消失，而是通过繁殖延迟等方式在未来几期继续影响种群数量，这就可以用MA项来建模。

此外，MA模型还常用于修正AR模型残差中的短期相关性。由于MA(q)的自相关函数（ACF）在滞后 $ q $ 步后截尾，这一特性成为判断MA阶数的重要依据。

下表总结了AR与MA模型的主要区别：

特征	AR(p) 模型	MA(q) 模型
平稳性条件	存在（特征根在单位圆外）	总是平稳
ACF 行为	拖尾（指数衰减或振荡）	在 $ q $ 步后截尾
PACF 行为	在 $ p $ 步后截尾	拖尾
参数解释	历史值的影响	过去扰动的影响
记忆长度	无限（理论上）	有限（最多 $ q $ 期）

该对比有助于在实际建模中根据ACF/PACF图形特征初步判断模型类型。

以下是一个MA(3)过程的Python模拟代码：

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import lfilter

np.random.seed(42)
eps = np.random.normal(0, 1, 500)  # 白噪声序列

# MA(3): Y_t = ε_t + 0.5ε_{t-1} - 0.3ε_{t-2} + 0.2ε_{t-3}
theta = [1, 0.5, -0.3, 0.2]
y_ma3 = lfilter(theta, 1, eps)

plt.figure(figsize=(10, 4))
plt.plot(y_ma3)
plt.title("Simulated MA(3) Process")
plt.xlabel("Time")
plt.ylabel("Y_t")
plt.grid(True)
plt.show()

参数说明与逻辑分析：
- lfilter(b, a, x) 实现数字滤波，此处 $ b = \theta $ 为MA系数向量，$ a = 1 $ 表示无自回归分母。
- 系数 [1, 0.5, -0.3, 0.2] 对应 $ \theta_0=1, \theta_1=0.5, \theta_2=-0.3, \theta_3=0.2 $。
- 输出序列显示出短记忆特征：单个冲击只影响未来3期。

该代码可用于验证MA模型的截尾性质，辅助学生理解其与AR的根本差异。

### 2.1.3 差分整合（I）在时间序列平稳化中的作用

差分操作是ARIMA模型区别于ARMA的关键所在。设原始序列为 $ X_t $，其一阶差分为：

Y_t = (1 - B)X_t = X_t - X_{t-1}

若一阶差分仍不平稳，可进行二阶差分：

Z_t = (1 - B)^2 X_t = (X_t - X_{t-1}) - (X_{t-1} - X_{t-2}) = X_t - 2X_{t-1} + X_{t-2}

依此类推，直到差分后的序列满足平稳性为止。这里的整数 $ d $ 即为差分阶数，表示需要多少次差分才能使序列平稳。

差分的本质是一种去趋势手段。对于含有线性趋势的序列，一阶差分即可去除；对于二次趋势，则需二阶差分。然而，过度差分会带来信息损失和方差放大风险，因此 $ d $ 的选择需谨慎。

下图展示了差分操作的效果（使用Mermaid流程图表示）：

graph TD
    A[原始非平稳序列 X_t] --> B{是否存在趋势?}
    B -- 是 --> C[进行一阶差分 Y_t = X_t - X_{t-1}]
    C --> D{Y_t 是否平稳?}
    D -- 否 --> E[进行二阶差分 Z_t = Y_t - Y_{t-1}]
    E --> F{Z_t 是否平稳?}
    F -- 是 --> G[在Z_t上拟合ARMA(p,q)]
    F -- 否 --> H[尝试更高阶差分或考虑季节差分]
    D -- 是 --> G
    B -- 否 --> I[直接拟合ARMA(p,q)]

该流程图清晰地描绘了ARIMA建模中差分决策的逻辑路径，体现了“先平稳化，再建模”的核心原则。

在实际操作中，差分阶数 $ d $ 的确定通常依赖于单位根检验（如ADF检验）与经验观察相结合的方法。下一节将详细讨论非平稳序列的识别与分类标准。

2.2 非平稳时间序列的识别与分类

非平稳性是时间序列建模中最常见的挑战之一。一个非平稳序列的统计特性（如均值、方差、自协方差）随时间变化，导致传统基于平稳假设的ARMA模型失效。因此，在构建ARIMA模型之前，必须准确识别非平稳性的类型及其成因，并采取相应的差分或其他变换手段予以处理。

非平稳性主要可分为三类：趋势性（Trend Stationarity）、季节性（Seasonality）和结构性突变（Structural Breaks）。每一类都有其独特的表现形式和对应的处理策略。此外，单位根检验（Unit Root Test）提供了形式化的统计工具，用于判断序列是否包含单位根（即是否为随机游走过程），进而指导差分阶数的选择。

### 2.2.1 趋势性、季节性与结构性突变的判别标准

趋势性 指的是序列均值随时间呈现系统性上升或下降的模式。它可以是确定性的（如线性趋势 $ \mu_t = a + bt $）或随机的（如随机游走 $ X_t = X_{t-1} + \varepsilon_t $）。前者可通过去趋势（detrending）方法处理，后者则需差分。

季节性 表现为固定周期内的重复波动，如月度销售额每年冬季升高。这类非平稳性可通过季节差分（Seasonal Differencing）消除，例如对月度数据使用 $ (1 - B^{12})X_t $。

结构性突变 是指由于政策变更、突发事件等原因导致序列统计特性发生突然改变。例如，疫情爆发导致旅游收入骤降。此类问题不能简单通过差分解决，往往需要引入虚拟变量或分段建模。

判别这些非平稳类型的常用方法包括：
- 可视化分析 ：绘制时间序列图，观察是否存在长期趋势、周期性波动或断点；
- 滚动统计量 ：计算滑动窗口内的均值与方差，查看其是否随时间变化；
- 频谱分析 ：检测是否存在显著的周期成分。

以下表格归纳了三种非平稳类型的特征与处理方式：

类型	主要特征	判别方法	处理手段
趋势性	均值随时间单调变化	时间图、ADF检验	一阶差分、去趋势
季节性	固定周期重复波动	ACF周期峰值、季节图	季节差分、SARIMA
结构性突变	方差/均值突变点	断点检验、残差突变	虚拟变量、分段拟合

### 2.2.2 单位根检验的基本思想与常见类型（ADF、PP）

单位根检验旨在判断序列是否具有单位根，即是否属于 $ I(1) $ 过程。最常见的两种检验是 Augmented Dickey-Fuller (ADF) 检验和 Phillips-Perron (PP) 检验。

ADF检验基于回归模型：

\Delta X_t = \alpha + \beta t + \gamma X_{t-1} + \sum_{i=1}^{k} \delta_i \Delta X_{t-i} + \varepsilon_t

原假设 $ H_0: \gamma = 0 $（存在单位根，非平稳），备择假设 $ H_1: \gamma < 0 $（平稳）。若检验统计量小于临界值，则拒绝原假设，认为序列平稳。

PP检验与ADF类似，但对异方差和序列相关更具鲁棒性，适用于误差项存在复杂相关结构的情形。

Python中可通过 statsmodels 库执行ADF检验：

from statsmodels.tsa.stattools import adfuller
import numpy as np

# 模拟一个带趋势的非平稳序列
np.random.seed(42)
t = np.arange(100)
x = 0.5 * t + np.cumsum(np.random.normal(0, 1, 100))  # 趋势+随机游走

# 执行ADF检验
result = adfuller(x, regression='ct')  # 包含常数项和时间趋势
print(f'ADF Statistic: {result[0]:.4f}')
print(f'p-value: {result[1]:.4f}')
print(f'Critical Values:')
for k, v in result[4].items():
    print(f'\t{k}: {v:.3f}')

输出解释：
- 若 p-value < 0.05 ，拒绝原假设，认为序列平稳；
- 若ADF统计量大于临界值（如1%水平下的-3.5），则无法拒绝单位根存在。

该代码可用于实际数据分析前的平稳性筛查。

### 2.2.3 差分阶数d的选择原则与实际案例分析

选择合适的差分阶数 $ d $ 是ARIMA建模的关键步骤。一般建议：
- $ d = 0 $：序列本身平稳；
- $ d = 1 $：大多数经济时间序列（如GDP增长率）；
- $ d = 2 $：极少使用，易造成过度差分。

实践中可通过以下步骤确定 $ d $：
1. 观察原始序列图是否有明显趋势；
2. 进行ADF检验判断是否平稳；
3. 若不平稳，进行一阶差分后再检验；
4. 直到差分后序列通过平稳性检验。

案例：某城市月度用电量数据（含趋势和季节性）

import pandas as pd
from statsmodels.tsa.seasonal import seasonal_decompose

# 假设有time_series数据
# decomposition = seasonal_decompose(time_series, model='additive', period=12)
# decomposition.plot()

通过分解图可识别趋势与季节成分，进而决定使用 $ d=1 $ 加上季节差分 $ D=1 $，构建SARIMA模型。

2.3 ARIMA建模流程的逻辑演进

### 2.3.1 从ARMA到ARIMA：模型扩展的必要性

ARMA模型仅适用于平稳序列，而现实中多数时间序列（如股价、气温、人口）均表现出非平稳性。ARIMA通过引入差分算子 $ (1-B)^d $，将非平稳序列转换为平稳序列，从而继承ARMA的优点并拓展其适用范围。

### 2.3.2 模型构建三阶段：识别、估计、诊断的闭环设计

识别阶段 ：通过ACF/PACF图和单位根检验确定 $ p,d,q $；
估计阶段 ：使用极大似然法估计参数；
诊断阶段 ：检验残差是否为白噪声（Ljung-Box检验）。

形成“建模—检验—修正”闭环。

### 2.3.3 实际应用场景对比：经济数据预测 vs 生物种群动态模拟

场景	数据特点	常用模型	关键挑战
经济数据（GDP）	趋势+季节性	SARIMA	政策干预影响
生物种群	非线性增长+环境噪声	ARIMA+外生变量	生态阈值突变

两者均需结合领域知识调整模型结构。

3. 时间序列的平稳性检验与差分技术

在现代时间序列建模中，平稳性不仅是统计推断有效性的前提，更是模型参数估计、预测性能和残差诊断的基础。对于ARIMA类模型而言，其核心假设之一是数据生成过程必须具备某种形式的平稳性。然而，现实世界中的大多数经济、生态或工程时间序列往往呈现出明显的趋势、季节性和结构性突变，表现出非平稳特征。因此，在进入正式建模之前，对原始序列进行系统的平稳性评估，并通过适当的差分技术实现平稳化处理，成为不可绕过的前置步骤。本章将深入探讨平稳性的数学定义及其在建模中的关键作用，系统解析主流的平稳性检验方法（如ADF、PP、KPSS），并详细阐述差分操作的技术细节与实践注意事项，尤其关注一阶差分、高阶差分与季节差分的操作差异、过度差分带来的信息损失风险，以及差分后预测结果如何逆变换还原为原始尺度。

3.1 平稳性的定义及其在建模中的核心地位

3.1.1 弱平稳与严平稳的区别及统计含义

在时间序列分析中，“平稳性”是指随机过程的统计特性不随时间推移而变化。根据约束强度不同，可分为 严平稳 （Strict Stationarity）和 弱平稳 （Weak Stationarity，又称宽平稳或协方差平稳）。严平稳要求整个联合概率分布函数在时间平移下保持不变，即对于任意时间点 $ t_1, t_2, …, t_n $ 和任意时间偏移 $ k $，有：
F_{X_{t_1},…,X_{t_n}} = F_{X_{t_1+k},…,X_{t_n+k}}
这一条件极为严格，实际应用中难以验证且多数真实过程无法满足。

相比之下，弱平稳仅要求前两阶矩稳定：
1. 均值恒定：$ E(X_t) = \mu $ 对所有 $ t $ 成立；
2. 方差有限且不变：$ Var(X_t) = \sigma^2 < \infty $；
3. 自协方差仅依赖于滞后阶数 $ h $：$ Cov(X_t, X_{t+h}) = \gamma(h) $

该定义更适用于实证建模，尤其是线性ARMA模型的理论框架正是建立在弱平稳假设之上。若违背此假设，模型参数估计将出现偏误，标准误失真，导致显著性检验失效。例如，在存在单位根的情况下，OLS估计量不再服从渐近正态分布，t统计量发散，传统推断完全失效。

特征	严平稳	弱平稳
分布稳定性	联合分布不变	一、二阶矩不变
数学要求	极高	可接受
实际可测性	几乎不可检验	可通过样本均值、方差、ACF图等判断
模型适用性	理论研究为主	ARMA/ARIMA建模基础

graph TD
    A[时间序列] --> B{是否平稳？}
    B -->|否| C[差分/变换]
    B -->|是| D[直接建模]
    C --> E[新序列]
    E --> F{是否平稳？}
    F -->|否| C
    F -->|是| G[ARMA建模]

上述流程图清晰地展示了平稳性在建模流程中的“守门员”角色——只有通过平稳性检验的序列才能进入后续的模型识别与拟合阶段。

3.1.2 非平稳序列对参数估计的误导风险

当使用非平稳序列直接拟合ARMA模型时，会引发一系列严重后果。最典型的例子是“伪回归”（Spurious Regression）现象：即使两个毫无关联的时间序列都含有共同趋势（如随机游走），回归模型仍可能显示高度显著的相关系数，从而得出错误因果结论。

以两个独立生成的I(1)过程为例：

rng(123); % 固定随机种子
n = 100;
x = cumsum(randn(n,1));  % I(1) 序列 x
y = cumsum(randn(n,1));  % I(1) 序列 y

% 普通最小二乘回归
mdl = fitlm(x, y);
disp(mdl.Rsquared.Ordinary);   % 输出 R²
disp(mdl.Coefficients.pValue); % 查看p值

代码逻辑逐行解读：
- 第2行：设置随机种子确保结果可复现；
- 第3行：设定样本长度为100；
- 第4–5行：利用累积和生成两个独立的随机游走序列（即一阶单整过程）；
- 第8行：执行线性回归 y ~ x ；
- 第9–10行：输出决定系数R²和回归系数的p值。

参数说明与扩展分析：
运行上述代码通常会发现，尽管 x 与 y 完全无关，但R²可能高达0.5以上，且斜率项p值小于0.05的概率远高于名义水平。这表明传统回归方法在非平稳条件下会产生误导性推断。

进一步地，在ARMA建模中，若对非平稳序列强行拟合ARMA(p,q)，极大似然估计（MLE）将失去一致性，预测区间大幅低估不确定性，模型自相关结构被扭曲。因此，必须先通过差分或其他去趋势手段使序列平稳后再建模。

此外，非平稳性还会破坏白噪声残差假设。例如，在未差分的趋势序列上拟合AR(1)模型，其残差往往会呈现系统性模式（如持续上升或周期波动），违反ARMA模型对残差为白噪声的基本要求，进而影响Ljung-Box检验的有效性。

3.1.3 可视化手段辅助判断趋势与波动特征

图形分析是初步识别非平稳性的直观工具。常见的可视化方法包括：

时间轨迹图（Time Plot） ：观察是否存在长期趋势、季节性波动或结构性断点。
滚动均值与方差图（Rolling Statistics） ：计算滑动窗口内的均值与方差，检查其是否随时间变化。
QQ图与直方图 ：辅助判断分布形态是否稳定。

以下MATLAB代码展示如何绘制滚动统计量来探测非平稳性：

% 生成带趋势的非平稳序列
t = 1:200;
yt = 0.5*t + randn(200,1)*2;  % 线性趋势 + 噪声

% 计算滚动均值和标准差（窗口大小=20）
window = 20;
roll_mean = movmean(yt, window);
roll_std  = movstd(yt, window);

% 绘图
figure;
subplot(2,1,1);
plot(yt); hold on;
plot(roll_mean, 'r', 'LineWidth', 2);
title('原始序列与滚动均值');
legend('原始数据', '滚动均值');
xlabel('时间'); ylabel('数值');

subplot(2,1,2);
plot(roll_std);
title('滚动标准差');
xlabel('时间'); ylabel('标准差');

代码逻辑逐行解读：
- 第2–3行：构造一个具有线性增长趋势的时间序列；
- 第6–7行：使用 movmean 和 movstd 函数计算宽度为20的移动平均与移动标准差；
- 第10–19行：上下排列两个子图，分别显示滚动均值与滚动标准差的变化趋势。

参数说明与扩展分析：
从图中可以看出，滚动均值呈明显上升趋势，说明原序列均值非平稳；滚动标准差虽略有波动但相对稳定，提示方差可能平稳。这种分离式的不平稳（仅均值变）正是差分法擅长处理的情形。

值得注意的是，视觉判断虽快捷但主观性强，尤其在小样本或高噪声情况下容易误判。因此，应结合正式的统计检验（如ADF、KPSS）进行综合决策。

3.2 常用平稳性检验方法详解

3.2.1 ADF检验的原假设设定与拒绝域解释

增强迪基-富勒检验（Augmented Dickey-Fuller Test, ADF）是最广泛使用的单位根检验方法，用于判断时间序列是否包含单位根（即是否为I(1)过程）。其基本思想是通过对以下回归模型进行参数显著性检验：

\Delta X_t = \alpha + \beta t + \gamma X_{t-1} + \sum_{i=1}^{p}\delta_i \Delta X_{t-i} + \varepsilon_t

其中：
- $ \Delta X_t = X_t - X_{t-1} $：一阶差分；
- $ \alpha $：截距项（表示漂移）；
- $ \beta t $：时间趋势项；
- $ \gamma $：待检验的关键参数，若 $ \gamma = 0 $ 则存在单位根。

原假设 H₀ : $ \gamma = 0 $（序列非平稳）
备择假设 H₁ : $ \gamma < 0 $（序列平稳）

检验统计量为 $ t_\gamma $，但由于在单位根下该统计量不服从标准t分布，需查专用临界值表（MacKinnon表）。若计算出的t值小于临界值（更负），则拒绝H₀，认为序列平稳。

MATLAB实现如下：

% 加载示例数据（可用simulated或real data）
Y = cumsum(randn(100,1)); % 模拟I(1)过程

% 执行ADF检验（含截距和趋势）
[h, pValue, stat, cValue] = adftest(Y, 'model', 'TS', 'lags', 1);

% 输出结果
fprintf('ADF检验结果:\n');
fprintf('  原假设成立？ %s\n', h==0 ? '是' : '否（拒绝）');
fprintf('  p值: %.4f\n', pValue);
fprintf('  统计量: %.4f\n', stat);
fprintf('  临界值(5%%): %.4f\n', cValue);

代码逻辑逐行解读：
- 第2行：模拟一个随机游走序列；
- 第5行：调用 adftest 函数，指定模型类型为’TS’（含截距和趋势），滞后阶数为1；
- 第8–12行：格式化输出检验结果。

参数说明与扩展分析：
- 'model' 可选 'AR' （无截距无趋势）、 'ARD' （有截距无趋势）、 'TS' （有截距有趋势），选择应基于序列图形特征；
- 'lags' 控制增广项数量，防止残差自相关；
- 返回值 h=1 表示拒绝原假设，即序列平稳。

3.2.2 PP检验对异方差情形的适应能力

菲利普斯-佩隆检验（Phillips-Perron, PP）与ADF类似，但采用非参数方式修正序列相关和异方差问题，无需显式引入滞后项。其回归形式为：

\Delta X_t = \gamma X_{t-1} + \varepsilon_t

但PP检验对误差项的自协方差结构进行直接调整，使得其在存在异方差或长记忆性时更具稳健性。

在MATLAB中可通过 pptest 函数实现：

[h, pVal] = pptest(Y, 'model', 'TS');
fprintf('PP检验p值: %.4f\n', pVal);

与ADF相比，PP更适合处理金融时间序列这类常具波动聚集性的数据。

检验方法	是否需选滞后阶数	异方差鲁棒性	适用场景
ADF	是	较弱	一般平稳性检验
PP	否	强	存在异方差或未知序列相关

3.2.3 KPSS检验作为补充分析工具的应用策略

KPSS检验（Kwiatkowski-Phillips-Schmidt-Shin）采取相反视角：其 原假设为序列平稳 （围绕趋势平稳），备择假设为存在单位根。因此，它常作为ADF/PP的补充工具，形成“双重检验”策略。

若：
- ADF拒绝H₀ 且 KPSS不拒绝H₀ → 支持平稳；
- ADF不拒绝H₀ 且 KPSS拒绝H₀ → 支持非平稳；
- 两者冲突 → 需进一步分析（如考虑结构断点）。

% KPSS检验（趋势平稳）
[h_kpss, p_kpss] = kpsstest(Y, 'trend', true);
fprintf('KPSS检验p值: %.4f (原假设：平稳)\n', p_kpss);

建议实践中同时报告多种检验结果，避免单一检验导致误判。

graph LR
    Start[开始检验] --> ADF[ADF检验]
    Start --> PP[PP检验]
    Start --> KPSS[KPSS检验]
    ADF --> Decision{综合判断}
    PP --> Decision
    KPSS --> Decision
    Decision --> Final[确定d阶差分]

3.3 差分操作的技术实现与注意事项

3.3.1 一阶差分、高阶差分与季节差分的操作区别

差分是消除趋势和单位根的主要手段。常见类型包括：

一阶差分 ：$ \nabla X_t = X_t - X_{t-1} $，消除线性趋势；
二阶差分 ：$ \nabla^2 X_t = \nabla(\nabla X_t) $，应对加速增长趋势；
季节差分 ：$ \nabla_s X_t = X_t - X_{t-s} $，s为季节周期（如月度数据s=12）。

MATLAB中可使用 diff 函数灵活实现：

% 示例：月度销售数据（含趋势+季节性）
load('sales_data.mat'); % 假设变量名为S，长度144（12年）

% 一阶差分
D1 = diff(S);

% 二阶差分
D2 = diff(S, 2);

% 季节差分（s=12）
Ds = diff(S, 1, 12); % 第三个参数指定滞后步长

% 季节+常规差分（SARIMA常用）
D_seasonal_then_regular = diff(Ds); % 先季节差分再一阶差分

% 可视化对比
figure;
subplot(2,2,1); plot(S); title('原始序列');
subplot(2,2,2); plot(D1); title('一阶差分');
subplot(2,2,3); plot(Ds); title('季节差分');
subplot(2,2,4); plot(D_seasonal_then_regular); title('双重差分');

代码逻辑逐行解读：
- 第4–7行：分别执行不同类型差分；
- diff(S, 1, 12) 表示沿时间轴对滞后12期做一次差分；
- 第10–15行：四象限图对比差分效果。

参数说明与扩展分析：
- diff(X,n,dim) 中 n 为差分阶数， dim 指定维度（默认1）；
- 实践中常采用“先季节差分再常规差分”的顺序，符合SARIMA建模范式。

3.3.2 过度差分导致信息损失的风险控制

虽然差分有助于平稳化，但 过度差分 （over-differencing）会导致：
- 增加噪声方差；
- 引入虚假的移动平均成分；
- 降低预测效率。

理论上，若真实过程为I(1)，对其进行两次差分会使模型隐含MA(1)项，即使原序列并无MA结构。

可通过以下方式规避：
1. 使用AIC/BIC准则比较不同差分阶数下的模型拟合优度；
2. 观察差分后ACF是否在滞后1处出现显著负相关（提示过度差分）；
3. 结合KPSS检验确认差分后是否已足够平稳。

3.3.3 差分后序列的逆变换还原预测值的方法

完成差分建模后，预测值处于差分空间，需通过 逆差分 （integration）还原至原始尺度。

以一阶差分为例：
\hat{X} t = \hat{\nabla X}_t + X {t-1}

递归公式为：
\hat{X} {n+h} = \hat{Z} {n+h} + \hat{X}_{n+h-1}
其中 $ Z_t = \nabla X_t $

MATLAB中可通过 forecast 自动处理，但手动实现有助于理解机制：

% 假设已获得差分序列预测值 Z_forecast (长度h)
% 初始值 use last observed level: X0 = S(end)

X_pred = zeros(h,1);
X_pred(1) = Z_forecast(1) + S(end);

for i = 2:h
    X_pred(i) = Z_forecast(i) + X_pred(i-1);
end

此递归累加过程即为“反积分”，也是ARIMA(p,d,q)中“I”部分的核心还原逻辑。

差分类型	逆变换方式	所需初始值
一阶差分	累加	最后一个观测值
季节差分	$ \hat{X} t = \hat{Z}_t + X {t-s} $	s期前的值
双重差分	分层还原	多个历史值

综上所述，差分不仅是技术操作，更是连接非平稳现实与平稳模型之间的桥梁。合理运用各类差分手段，并辅以严谨的检验与逆变换机制，是构建可靠ARIMA模型的关键所在。

4. 自相关图与偏自相关图分析

在时间序列建模中，自相关函数（ACF）和偏自相关函数（PACF）是识别模型类型与阶数的关键工具。尤其在ARMA/ARIMA类模型的构建过程中，ACF与PACF图形提供了直观且有效的模式识别依据，帮助研究者初步判断时间序列的生成机制。通过观察滞后项之间的相关性结构，可以区分AR、MA或混合型过程，并为后续参数估计提供方向性指导。然而，这种基于图形的判别方法并非绝对可靠，其有效性依赖于样本长度、噪声水平以及是否存在季节性干扰等因素。因此，在实际应用中必须结合统计检验与信息准则进行综合评估。

本章将深入剖析ACF与PACF的数学基础、图形特征及其在模型识别中的具体作用。首先从自相关图的构造原理出发，解析其置信区间的设定逻辑与典型模式表现；随后探讨偏自相关图如何剥离中间变量影响以提取“净”相关性，并说明其在确定AR阶数方面的独特优势；最后指出图形分析固有的局限性，提出通过多重验证手段提升模型定阶准确性的策略体系。整个分析过程不仅强调理论推导的严谨性，更注重与MATLAB等计算平台的实际操作对接，确保理论理解能够无缝转化为可执行的技术流程。

4.1 ACF图的构建原理与模式识别

自相关函数（Autocorrelation Function, ACF）用于衡量一个时间序列在不同时刻的观测值之间的线性相关程度。它描述的是当前时刻 $ y_t $ 与过去某一滞后 $ k $ 时刻 $ y_{t-k} $ 之间的皮尔逊相关系数，反映了序列内部的记忆特性。在ARIMA建模框架下，ACF图是识别移动平均（MA）成分的重要依据，尤其是通过观察其“截尾”行为来判断MA(q)模型的阶数 $ q $。

4.1.1 自相关系数计算公式及其置信区间判定

对于一个弱平稳时间序列 $ {y_t} $，其滞后 $ k $ 的自相关系数定义为：

\rho_k = \frac{\text{Cov}(y_t, y_{t-k})}{\sqrt{\text{Var}(y_t)\text{Var}(y_{t-k})}} = \frac{\gamma_k}{\gamma_0}

其中 $ \gamma_k $ 是滞后 $ k $ 的协方差，$ \gamma_0 $ 是方差。在实际计算中，使用样本估计值：

\hat{\rho} k = \frac{\sum {t=k+1}^n (y_t - \bar{y})(y_{t-k} - \bar{y})}{\sum_{t=1}^n (y_t - \bar{y})^2}

该公式体现了标准化后的协方差关系，取值范围为 $[-1, 1]$，正负号表示相关方向，绝对值大小反映强度。

为了判断某滞后阶数的相关性是否显著异于零，需引入置信区间。通常采用 Bartlett 近似标准误法构建近似置信带：

SE(\hat{\rho} k) \approx \sqrt{\frac{1 + 2\sum {i=1}^{k-1} \hat{\rho}_i^2}{n}}

则95%置信区间为：

\pm 1.96 \times SE(\hat{\rho}_k)

若某个 $ \hat{\rho}_k $ 超出此区间，则认为该滞后项存在显著自相关。

下面是在 MATLAB 中绘制 ACF 图的示例代码：

% 模拟一个 MA(2) 时间序列
rng(1);
theta = [0.8 -0.5];
ma2 = arima('MA', theta, 'Constant', 0, 'Variance', 1);
y = simulate(ma2, 200);

% 绘制自相关图
figure;
autocorr(y, 20);
title('Sample ACF of Simulated MA(2) Process');

代码逻辑逐行解读：

rng(1); ：设置随机种子，保证结果可复现。
theta = [0.8 -0.5]; ：定义 MA(2) 模型的两个参数。
ma2 = arima(...) ：创建一个无常数项、方差为1的 MA(2) 模型对象。
simulate(ma2, 200) ：生成长度为200的时间序列样本。
autocorr(y, 20) ：调用内置函数绘制前20个滞后的 ACF 图，自动添加置信带（虚线）。

输出图像显示，前两个滞后显著非零，之后迅速衰减至置信带内，符合 MA(2) 的理论预期——ACF 在滞后 $ q $ 后“截尾”。

滞后阶数	样本自相关系数	理论行为（MA(q)）	判别意义
1	0.52	显著	支持 MA 成分
2	-0.31	显著	可能 $ q=2 $
≥3	<0.1	不显著	截尾示意

该表总结了不同滞后阶数下的观察结果与模型含义，辅助快速识别潜在模型结构。

graph TD
    A[输入原始时间序列] --> B[计算各滞后k的样本自相关系数ρ̂_k]
    B --> C[估算标准误SE(ρ̂_k)]
    C --> D[构建95%置信区间±1.96×SE]
    D --> E[绘制柱状ACF图并标注置信带]
    E --> F[判断拖尾/截尾模式]
    F --> G[初步识别MA阶数q]

上述流程图清晰展示了从数据输入到模式识别的完整 ACF 分析路径，强调每一步的技术实现与决策节点。

4.1.2 拖尾与截尾行为的直观表现与理论依据

在 ACF 图中，“拖尾”指自相关系数随滞后增加缓慢衰减（如指数或振荡式下降），而“截尾”则表现为超过某一滞后 $ q $ 后所有 $ \hat{\rho}_k ≈ 0 $。这两种行为分别对应不同的数据生成机制。

MA(q) 模型具有 ACF 截尾性 ：由于 MA(q) 表达式为
$$
y_t = \mu + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \cdots + \theta_q \varepsilon_{t-q}
$$
其协方差仅在 $ k ≤ q $ 时非零，故当 $ k > q $ 时 $ \gamma_k = 0 $，导致 $ \rho_k = 0 $，即理论上严格截尾。尽管样本估计会有波动，但应在置信带内。
AR(p) 模型呈现 ACF 拖尾 ：AR(1) 如 $ y_t = \phi y_{t-1} + \varepsilon_t $，其 $ \rho_k = \phi^k $，呈几何衰减。高阶 AR 同样表现为渐进衰减而非突然归零。

因此，若 ACF 图显示前几阶显著而后趋于零，提示可能是 MA 或 ARMA 模型；若持续缓慢下降，则更可能为 AR 过程。

例如，考虑一个真实 GDP 增长率序列，其 ACF 缓慢递减，表明长期记忆效应较强，适合用 AR 模型捕捉趋势惯性；而某季度销售额的误差修正项若只影响未来两期，则其 ACF 应在滞后2后截断，支持 MA(2) 设定。

4.1.3 季节性成分在ACF中的周期性峰值识别

当时间序列包含季节性模式（如月度数据中的年度周期），ACF 图会表现出规律性的峰值，出现在 $ k = s, 2s, 3s, \dots $ 处，其中 $ s $ 为季节周期（如12个月）。这些周期性突出的相关性是识别季节性模型（SARIMA）的关键线索。

例如，某零售额数据每年同月出现高峰，其 ACF 将在滞后12、24、36等位置出现显著正值，形成“季节峰簇”。此时即使整体趋势被差分消除，季节性自相关仍保留在残差中，提示需要引入季节性MA或AR项。

MATLAB 实现如下：

% 加载含季节性的模拟数据
load Data_Airline; % 国际旅客数量对数序列
y_log = log(Data); 
diff_y = diff(y_log);         % 一阶差分
seasonal_diff = diff(diff_y, 1, 12); % 再做12步季节差分

% 绘制差分后序列的ACF
figure;
autocorr(seasonal_diff, 24);
title('ACF after Seasonal Differencing');

执行后可见，原本在滞后12、24处的显著峰值大幅减弱，说明季节差分有效去除了周期性相关。反之，若未充分处理季节性，ACF中残留的周期峰将误导模型定阶。

此外，可通过以下表格归纳常见模型的 ACF 特征：

模型类型	ACF 行为	PACF 行为	典型应用场景
AR(p)	拖尾	在p阶截尾	经济增长率、利率变动
MA(q)	在q阶截尾	拖尾	误差反馈、冲击传导
ARMA(p,q)	双重拖尾	双重拖尾	复杂动态系统建模
季节MA(S)	在s,2s,…截尾	周期性拖尾	月度销售、气温变化
季节AR(S)	周期性拖尾	在s阶截尾	年度财政支出周期

此表为建模初期的模式匹配提供参考基准，增强图形判读的系统性。

4.2 PACF图的作用机制与解读技巧

偏自相关函数（Partial Autocorrelation Function, PACF）用于衡量在控制中间滞后项影响后，当前值 $ y_t $ 与滞后 $ k $ 值 $ y_{t-k} $ 之间的直接线性相关性。如果说 ACF 描述的是“总相关”，那么 PACF 提取的是“净相关”。这一特性使其成为识别自回归（AR）模型阶数 $ p $ 的核心工具。

4.2.1 控制中间滞后项影响下的净相关性提取

PACF 的数学本质来源于多元回归模型。第 $ k $ 阶偏自相关系数 $ \alpha_k $ 定义为以下回归中最后一个系数：

y_t = \phi_{k1} y_{t-1} + \phi_{k2} y_{t-2} + \cdots + \phi_{kk} y_{t-k} + \varepsilon_t

其中 $ \phi_{kk} $ 即为 $ \text{PACF}(k) $。它表示在已知 $ y_{t-1}, \dots, y_{t-k+1} $ 的条件下，$ y_{t-k} $ 对 $ y_t $ 的额外解释能力。

该机制可通过 Yule-Walker 方程组求解，也可利用 Durbin-Levinson 算法递归计算：

\alpha_k = \frac{\rho_k - \sum_{j=1}^{k-1} \alpha_j \rho_{k-j}}{1 - \sum_{j=1}^{k-1} \alpha_j \rho_j}

MATLAB 示例：

% 使用同一MA(2)序列计算PACF
figure;
parcorr(y, 20);
title('Sample PACF of Simulated MA(2) Process');

结果显示 PACF 缓慢衰减，无明显截断，符合 MA 过程特征——PACF 拖尾。

相反，若为 AR(2) 过程：

ar2 = arima('AR', [0.6 -0.3], 'Constant', 0, 'Variance', 1);
y_ar = simulate(ar2, 200);
figure;
parcorr(y_ar, 20);
title('PACF of AR(2) Process');

PACF 在滞后2后落入置信带，呈现典型“截尾”现象，支持 $ p=2 $ 的判断。

4.2.2 截尾特性用于AR阶数p的初步判断

AR(p) 模型的理论性质决定了其 PACF 在 $ k > p $ 时为零，即严格截尾。这是由于 AR(p) 的动态依赖仅限于前 $ p $ 期，超出部分已被完全吸收。

因此，观察 PACF 图中首次进入置信带的滞后点，常作为选择 $ p $ 的经验法则。例如：
- 若 PACF 在滞后1显著，之后均不显著 → 初步判断为 AR(1)
- 若滞后1、2显著，第3起不显著 → 考虑 AR(2)

但应注意，小样本可能导致虚假截尾或延迟截尾，应辅以信息准则进一步确认。

4.2.3 ACF与PACF联合使用进行模型类型甄别

单独依赖任一图表易造成误判，最佳实践是将 ACF 与 PACF 结合分析。经典的模型识别规则如下表所示：

模型类型	ACF 行为	PACF 行为	推荐模型
白噪声	全部不显著	全部不显著	无需建模
AR(p)	拖尾	在p阶截尾	AR(p)
MA(q)	在q阶截尾	拖尾	MA(q)
ARMA(p,q)	双重拖尾	双重拖尾	尝试ARMA(p,q)

例如，某金融收益率序列显示 ACF 缓慢衰减、PACF 在滞后3后截尾，则优先尝试 AR(3)；若 ACF 在滞后2后截尾而 PACF 拖尾，则倾向于 MA(2)。

graph LR
    A[绘制ACF与PACF图] --> B{ACF截尾?}
    B -- 是 --> C[PACF拖尾 → MA(q)]
    B -- 否 --> D{PACF截尾?}
    D -- 是 --> E[ACF拖尾 → AR(p)]
    D -- 否 --> F[两者均拖尾 → ARMA(p,q)]

该决策流程图实现了从图形特征到模型假设的自动化推理路径，适用于快速原型设计阶段。

此外，可通过以下代码实现双图并列展示：

figure;
tiledlayout(2,1);
nexttile;
autocorr(y_ar, 15);
nexttile;
parcorr(y_ar, 15);
title('ACF and PACF of AR(2) Series');

便于对比分析两种函数的行为差异。

4.3 图形分析的局限性与补充手段

尽管 ACF/PACF 图在时间序列建模中广泛应用，但其主观性强、对样本敏感，存在诸多限制。忽视这些问题可能导致模型误设，进而影响预测精度。

4.3.1 小样本下图形误判的可能性评估

当样本量较小时（如 $ n < 50 $），样本自相关系数波动剧烈，容易出现假阳性或假阴性。例如，本应截尾的 MA(1) 过程可能因抽样误差在滞后3处出现显著相关，误导为更高阶模型。

研究表明，至少需要 $ n ≥ 50 $ 才能获得相对稳定的 ACF/PACF 估计。建议在小样本情况下优先使用信息准则（AIC/BIC）进行客观比较。

4.3.2 结合信息准则优化主观判断偏差

单纯依赖图形截尾点选择 $ p $ 或 $ q $ 容易受视觉偏差影响。推荐做法是：基于 ACF/PACF 初步设定候选模型集合（如 AR(1)~AR(4)），然后比较其 AIC 或 BIC 值，选择最小者。

例如：

models = cell(4,1);
aic_vals = zeros(4,1);
for p = 1:4
    mdl = estimate(arima(p,0,0), y_ar);
    aic_vals(p) = aicbic(mdl.LogLikelihood, p, length(y_ar));
    models{p} = mdl;
end
[~, best_p] = min(aic_vals);
fprintf('Optimal AR order by AIC: %d\n', best_p);

该段代码自动搜索最优 AR 阶数，避免人为选择带来的偏差。

4.3.3 多重检验交叉验证提升识别准确性

为进一步提高可靠性，可结合单位根检验、Ljung-Box 检验、残差诊断等多种方法形成交叉验证体系。例如：
- 差分前后 ACF 变化是否合理？
- 拟合后残差是否通过白噪声检验？
- 不同差分阶数下模型 AIC 是否稳定？

建立如下综合评估表：

检验项目	方法	达标标准
平稳性	ADF检验	p < 0.05
残差自相关	Ljung-Box检验	p > 0.05（无剩余相关）
参数显著性	t检验
模型复杂度	BIC	最小化
预测稳定性	滚动窗口回测	RMSE变化平稳

通过多维度验证，可显著降低图形误判风险，提升建模稳健性。

5. AR与MA阶数的确定方法

在时间序列建模中，正确识别自回归（AR）和移动平均（MA）部分的阶数是构建高效、稳健ARIMA模型的关键步骤。阶数选择不当会导致模型过拟合或欠拟合，从而严重影响预测精度与统计推断的有效性。尤其在实际应用中，如宏观经济指标分析、能源需求预测或生物种群动态模拟等场景，对模型结构的敏感度极高。因此，如何科学地确定AR(p)中的 $ p $ 和 MA(q)中的 $ q $ 成为建模流程中承上启下的核心环节。

传统方法依赖于对自相关函数（ACF）和偏自相关函数（PACF）图形的主观判读，这类经验法则虽然直观易用，但在面对复杂噪声结构或多重周期干扰时容易产生误判。随着计算能力的发展，基于信息准则的自动化定阶方法逐渐成为主流。AIC（Akaike Information Criterion）、BIC（Bayesian Information Criterion）等指标通过平衡拟合优度与模型复杂度，在候选模型集中筛选出最优配置。然而，即便采用这些客观标准，仍需结合残差诊断与邻近阶数组合的对比实验进行验证，以确保最终选定模型具备良好的泛化能力和参数稳定性。

此外，现代时间序列分析强调“模型不确定性”的管理，即不应仅依赖单一最优模型，而应通过敏感性测试评估不同阶数组合的表现差异。例如，在金融波动率建模中，GARCH类模型常与ARIMA结合使用，其前期ARMA结构的阶数微小变化可能显著影响后续波动聚集效应的捕捉能力。因此，完整的阶数确定过程必须包含多个层次的交叉验证机制：从初步图形判断到信息准则排序，再到参数显著性检验与残差白噪声检测，形成闭环反馈系统。

本章将系统阐述三种主要的阶数识别路径：基于ACF/PACF的经验规则、信息准则驱动的自动搜索策略，以及多维度稳健性验证手段。重点剖析每种方法的理论依据、适用边界及其潜在陷阱，并结合MATLAB代码实现关键操作步骤，展示如何在真实数据背景下综合运用多种工具完成精确建模。通过深入理解这些技术细节，读者不仅能够掌握标准流程，还能在面对非典型序列时灵活调整策略，提升整体建模水平。

5.1 基于ACF/PACF的经验法则应用

在经典时间序列分析框架中，自相关图（ACF）与偏自相关图（PACF）是最基础且广泛使用的图形化工具，用于初步判断ARIMA模型中AR和MA部分的阶数。这种方法源于Box-Jenkins建模哲学，强调通过观察样本相关性结构来推断潜在的数据生成机制。尽管其主观性强、易受噪声干扰，但在数据质量较高、趋势清晰的情况下，仍能提供极具价值的方向性指引。

5.1.1 AR(p)模型阶数p的PACF截尾判断法

对于一个纯自回归过程AR(p)，其数学表达式为：

X_t = \phi_1 X_{t-1} + \phi_2 X_{t-2} + \cdots + \phi_p X_{t-p} + \varepsilon_t

其中 $\varepsilon_t \sim WN(0, \sigma^2)$ 是白噪声误差项。该模型的核心特性在于当前值仅由前 $ p $ 个滞后值线性组合而成。这一结构性特征反映在偏自相关函数（PACF）上表现为“截尾”现象——即当滞后阶数 $ k > p $ 时，PACF系数趋于零并在置信区间内随机波动。

具体而言，PACF衡量的是在控制中间所有滞后项（$ X_{t-1}, …, X_{t-k+1} $）影响后，$ X_t $ 与 $ X_{t-k} $ 之间的净相关性。在AR(p)过程中，一旦滞后超过 $ p $ 阶，这种直接依赖关系消失，因此PACF会在第 $ p $ 阶之后迅速衰减至不显著水平。利用这一点，可通过观察PACF图中显著超出±1.96/$\sqrt{n}$ 置信带的峰值数量来估计 $ p $。

以下是在MATLAB中绘制PACF并辅助判断AR阶数的示例代码：

% 模拟一个AR(2)过程
rng(123);
n = 500;
ar_coeffs = [1.5, -0.7]; % phi1=1.5, phi2=-0.7
X = arima('Constant', 0, 'AR', ar_coeffs, 'Variance', 1);
data = simulate(X, n);

% 绘制PACF
figure;
parcorr(data, 'NumLags', 15, 'NumSTD', 1.96);
title('Partial Autocorrelation Function (PACF) of AR(2) Process');

代码逻辑逐行解读：

rng(123); ：设置随机种子以保证结果可复现。
n = 500; ：设定模拟序列长度为500个观测点。
ar_coeffs = [1.5, -0.7]; ：定义AR(2)模型的两个自回归系数。
X = arima(...) ：创建一个无常数项、方差为1的AR(2)模型对象。
data = simulate(X, n); ：生成符合该模型的时间序列数据。
parcorr(data, 'NumLags', 15, 'NumSTD', 1.96); ：绘制前15阶的PACF图，置信区间设为95%（对应1.96倍标准误）。

执行上述代码后，预期看到前两阶PACF显著高于置信带，第三阶起落入区间内，支持 $ p=2 $ 的判断。

滞后阶数	PACF值	是否显著
1	1.48	是
2	-0.62	是
3	0.07	否
4	-0.03	否

说明：此表为模拟结果的近似输出，用于辅助可视化判断。

graph TD
    A[输入时间序列数据] --> B{是否平稳？}
    B -- 否 --> C[进行差分处理]
    C --> D[得到平稳序列]
    B -- 是 --> D
    D --> E[计算PACF]
    E --> F[观察截尾位置]
    F --> G[确定AR阶数p]

该流程图展示了基于PACF判断AR阶数的标准流程，强调了平稳性预处理的重要性。

5.1.2 MA(q)模型阶数q的ACF截尾判断法

移动平均模型MA(q)的形式为：

X_t = \varepsilon_t + \theta_1 \varepsilon_{t-1} + \cdots + \theta_q \varepsilon_{t-q}

其本质是当前值由当前及过去 $ q $ 期的白噪声冲击加权构成。由于误差项之间相互独立，MA(q)过程的自相关函数（ACF）在滞后 $ q+1 $ 及以后理论上为零，呈现出“截尾”行为。

相比之下，PACF则呈现拖尾（指数衰减或振荡衰减），不适合用于MA阶数识别。因此，实践中通常借助ACF图中显著相关性的最大滞后阶数来估计 $ q $。

下面是一个MA(1)过程的MATLAB实现：

% 模拟MA(1)过程
ma_coeffs = [0.6]; % theta1 = 0.6
Y = arima('Constant', 0, 'MA', ma_coeffs, 'Variance', 1);
data_ma = simulate(Y, n);

% 绘制ACF
figure;
autocorr(data_ma, 'NumLags', 15, 'NumSTD', 1.96);
title('Autocorrelation Function (ACF) of MA(1) Process');

参数说明与逻辑分析：

'MA', ma_coeffs ：指定MA部分的系数向量。
autocorr 函数自动计算并绘出样本自相关系数及其渐近置信区间。
若第一阶ACF显著，其余均不显著，则支持 $ q=1 $。

滞后阶数	ACF值	显著性
1	0.42	是
2	0.05	否
3	-0.04	否

注：由于理论ACF在MA(1)中为 $ \rho_1 = \frac{\theta}{1+\theta^2} = \frac{0.6}{1+0.36} \approx 0.44 $，模拟结果接近理论值。

graph LR
    P[原始序列] --> Q{平稳？}
    Q -- 否 --> R[差分]
    Q -- 是 --> S[计算ACF]
    S --> T[查找最后一个显著滞后期]
    T --> U[设定MA阶数q]

该流程突出ACF在MA建模中的主导作用。

5.1.3 ARMA(p,q)混合模型的双重拖尾识别难点

当模型同时包含AR和MA成分时（即ARMA(p,q)），ACF与PACF均表现为“拖尾”，使得单纯依靠图形截尾特征难以准确辨识 $ p $ 和 $ q $。此时，传统的经验法则失效风险显著上升。

例如，ARMA(1,1)过程：

X_t = \phi X_{t-1} + \varepsilon_t + \theta \varepsilon_{t-1}

其ACF和PACF均为指数衰减形式，无法通过简单计数判断阶数。在这种情况下，必须结合其他方法，如信息准则比较或多步预测误差评估。

一种常见策略是枚举多个 $ (p,q) $ 组合，分别拟合并比较AIC/BIC值。下表列出几种典型ARMA模型的相关函数模式：

模型类型	ACF 行为	PACF 行为
AR(p)	拖尾	在p阶截尾
MA(q)	在q阶截尾	拖尾
ARMA(p,q)	拖尾	拖尾

此表可用于快速对照分析。

综上所述，尽管ACF/PACF图形分析在教学和初级建模中具有重要地位，但其局限性不容忽视。尤其在高阶混合模型或存在外生变量的情境下，应将其视为初步探索工具，而非最终决策依据。后续章节将进一步引入更客观、系统的定阶方法予以补充和完善。

6. MATLAB中arima函数的应用与模型构建

时间序列建模在现代数据分析中扮演着核心角色，尤其在经济、金融、生态学和工程信号处理等领域，ARIMA（自回归积分移动平均）模型因其对趋势性、非平稳性和随机波动的良好适应能力而被广泛采用。随着计算工具的发展，MATLAB作为一款集数值计算、可视化与算法开发于一体的高级技术计算环境，在时间序列分析方面提供了强大且直观的支持。其中， arima 函数是实现ARIMA类模型定义与估计的核心工具之一。深入掌握该函数的使用方式，不仅有助于提升建模效率，更能确保参数配置合理、结果解释准确。

本章将系统阐述如何在MATLAB环境中完成从原始数据预处理到ARIMA模型构建的全过程。重点聚焦于 arima 类对象的创建机制、模型参数的灵活设定方法、以及利用 estimate 函数进行参数拟合的技术细节。通过结合实际操作代码、输出结构解析与常见问题调试策略，帮助具备一定统计背景的IT从业者或科研人员建立起完整的ARIMA建模工作流，并为后续章节中的完整案例打下坚实基础。

6.1 MATLAB时间序列对象与数据预处理

在正式进入ARIMA模型构建之前，必须首先完成对原始时间序列数据的有效组织与清洗。MATLAB提供多种数据容器用于管理带有时间戳的信息流，其中最常用的是 timetable 和 timeseries 类型。选择合适的数据结构不仅能提高代码可读性，还能显著增强后续差分、绘图与建模操作的稳定性。

6.1.1 timetable与timeseries数据类型的选用

timetable 是MATLAB R2016b之后引入的一种专为时间标记数据设计的表格类型，其最大优势在于支持多变量同步索引、自动对齐不同采样频率的时间点，以及丰富的内置函数集成。相比之下， timeseries 虽然也能存储单一时间序列及其时间向量，但功能相对有限，不推荐用于复杂多维场景。

以下是一个典型的 timetable 构建示例：

% 示例：构建包含生物种群数量的时间表
timeVector = datetime(2020,1,1):caldays(7):datetime(2023,12,31); % 每周采样
populationData = 500 + cumsum(randn(length(timeVector),1)); % 模拟增长趋势
tt = timetable(timeVector', populationData, 'VariableNames', {'Population'});
head(tt)

输出如下：

    timeVector         Population
    ____________       __________
    01-Jan-2020        498.63    
    08-Jan-2020        497.85    
    15-Jan-2020        496.21    
    ...

特性	`timetable`	`timeseries`
多变量支持	✅ 支持多个变量列	❌ 仅支持单个数据序列
时间对齐	✅ 自动对齐不同频率数据	⚠️ 需手动处理
差分操作兼容性	✅ 可直接调用diff	⚠️ 需提取数值后操作
与arima函数兼容性	✅ 推荐输入格式	⚠️ 不直接支持

逻辑分析 ：上述代码首先生成一个以周为单位的时间向量，模拟三年内的观测周期；接着通过累加正态噪声构造具有随机游走特征的种群数量序列。最终使用 timetable 将时间和数据绑定，形成结构化时间序列对象。这种结构便于后续调用 plot(tt) 实现自动时间轴绘制，也兼容大多数 Econometrics Toolbox 函数。

6.1.2 缺失值插补与异常点检测的内置函数调用

真实世界中的时间序列常存在缺失值或极端异常点，若不加以处理，会导致模型估计偏差甚至失败。MATLAB 提供了高效的内置函数来应对这些问题。

缺失值识别与插补

假设原始数据中存在若干 NaN 值：

tt.Population(10:10:50) = NaN; % 引入人工缺失
tt_clean = fillmissing(tt, 'linear'); % 线性插值填充

也可以使用更稳健的方法如 'spline' 或 'movmean' 进行局部平滑填补：

tt_filled = fillmissing(tt, 'movmean', 5); % 5期移动均值

异常值检测与修正

异常值可通过统计阈值法或箱线图原则识别：

% 使用isoutlier函数检测并替换异常值
[tt_clean.Outlier, tt_clean.Score] = isoutlier(tt_filled.Population);
tt_clean.Population(tt_clean.Outlier) = filloutliers(tt_filled.Population, 'center');

graph TD
    A[原始时间序列] --> B{是否存在NaN?}
    B -- 是 --> C[调用fillmissing]
    B -- 否 --> D{是否存在异常值?}
    C --> D
    D -- 是 --> E[调用isoutlier或filloutliers]
    E --> F[输出清洗后数据]
    D -- 否 --> F

参数说明 ：
- 'linear' ：基于相邻两点线性插值，适合连续变化序列。
- 'movmean' ：指定窗口大小的滑动平均，能保留趋势同时抑制噪声。
- filloutliers(..., 'center') ：将异常值替换为局部中心趋势值（如中位数），避免剧烈跳跃。

此类预处理流程应作为建模前的标准步骤执行，确保输入数据满足模型假设前提。

6.1.3 差分操作在MATLAB中的向量化实现

对于非平稳序列，需通过差分使其趋于平稳。MATLAB允许直接对 timetable 中的变量进行向量化差分运算。

% 一阶差分
tt_diff1 = diff(tt.Population);

% 转换回timetable以便保持时间信息
time_diff = tt.timeVector(2:end);
tt_d1 = timetable(time_diff, tt_diff1, 'VariableNames', {'D_Population'});

% 高阶差分（二阶）
tt_diff2 = diff(tt.Population, 2);

% 季节性差分（例如年度周期=52周）
seasonal_lag = 52;
tt_seas_diff = tt.Population(seasonal_lag+1:end) - tt.Population(1:end-seasonal_lag);

为了验证差分效果，通常结合可视化与单位根检验：

figure;
subplot(2,1,1);
plot(tt.timeVector, tt.Population); title('原始序列');
subplot(2,1,2);
plot(tt_d1.timeVector, tt_d1.D_Population); title('一阶差分后序列');

差分后的序列若呈现出围绕零均值波动、无明显趋势的特性，则可初步判断已实现平稳化。此外，可通过 adftest 函数进一步量化检验：

[h, pValue] = adftest(tt_d1.D_Population);
disp(['ADF检验p值: ', num2str(pValue)]);

执行逻辑说明 ：
- diff(X, n) 表示对向量 X 进行 n 阶差分，默认 n=1。
- 对 timetable 直接应用 diff 会丢失时间标签，因此需要手动重建时间索引。
- 季节性差分公式为 $ y_t^{(s)} = y_t - y_{t-s} $，适用于存在固定周期波动的情形。

综上所述，合理的数据预处理不仅是模型成功的前提，也是保证后续ACF/PACF分析准确性的重要保障。借助MATLAB强大的数据结构与函数库，整个流程可以高度自动化，极大提升了建模效率与可靠性。

6.2 arima类模型的定义与参数配置

一旦数据完成预处理，下一步便是定义ARIMA模型结构。MATLAB通过面向对象的方式封装了 arima 类，允许用户以声明式语法精确控制模型形式与参数约束。

6.2.1 使用arima(p,d,q)创建模型模板的语法规范

最基本的ARIMA模型可通过如下方式实例化：

Mdl = arima(2,1,1); % 定义ARIMA(2,1,1)模型
disp(Mdl)

输出显示模型结构详情：

    Description: "ARIMA(2,1,1) Model (Gaussian Distribution)"
    Distribution: Name = "Gaussian"
    P: 3 (含差分后滞后总数)
    D: 1
    Q: 1
    Constant: NaN
    AR: {NaN NaN} at lags [1 2]
    MA: {NaN} at lag [1]
    Seasonality: 0
    Beta: [1×0]
    Variance: NaN

关键字段解释如下：

属性	含义
`P`	总滞后阶数 = p + d + q_max_seasonal，影响最大滞后依赖
`D`	差分阶数，决定积分部分
`Constant`	模型常数项，差分后若d>0则通常设为0
`AR` / `MA`	自回归与移动平均系数数组，初始为NaN表示待估

值得注意的是，当 $ d > 0 $ 时，常数项往往趋近于零，因此建议显式设置：

Mdl.Constant = 0; % 对高阶差分序列禁用常数项

6.2.2 固定约束参数与自由参数的灵活设置

在某些先验知识指导下，可冻结特定参数以简化模型或避免过拟合。例如，强制某个AR系数为0：

Mdl.AR{2} = 0;     % 固定AR(2)系数为0
Mdl.AR{2} = NaN;   % 恢复为自由参数

也可批量设置多个参数：

Mdl = arima('ARLags',[1 3], 'MALags',1, 'D',1); % 跳跃滞后AR(1), AR(3)

这等价于建立稀疏AR结构，适用于存在间隔显著相关性的场景。

6.2.3 季节性ARIMA模型（SARIMA）的扩展定义方式

对于具有明显季节模式的数据（如月度销售、气温等），应采用SARIMA模型。其定义需额外指定季节性部分参数：

SMdl = arima('Seasonality',12,...           % 年度周期（月数据）
             'ARLags',[1], 'SARLags',[12],... % 非季节与季节AR
             'MALags',[1], 'SMALags',[12],...
             'D',1, 'SeasonalD',1,...
             'Constant',0);

该模型等价于 SARIMA(1,1,1)×(1,1,1)_12，完整表达式为：

(1 - \phi_1B)(1 - \Phi_1B^{12}) (1-B)(1-B^{12}) y_t = (1 + \theta_1B)(1 + \Theta_1B^{12}) \varepsilon_t

graph LR
    subgraph SARIMA Structure
        A[非季节AR(1)] --> C[乘积模型]
        B[季节AR(1)] --> C
        D[一阶差分] --> C
        E[季节差分] --> C
        F[非季节MA(1)] --> C
        G[季节MA(1)] --> C
    end
    C --> H[输出Y_t]

代码逻辑逐行解读 ：
- 'Seasonality',12 ：启用每12期一次的季节性差分。
- 'SARLags',[12] ：表示季节性自回归项位于滞后12处。
- 'SeasonalD',1 ：执行一次季节性差分 $ y_t - y_{t-12} $。
- 所有组件通过乘积形式耦合，构成复合模型。

通过这种方式，MATLAB能够灵活表达高度复杂的动态结构，极大增强了对现实世界周期性现象的刻画能力。

6.3 estimate函数执行参数估计的过程解析

定义好模型模板后，即可使用 estimate 函数对实际数据进行参数拟合。

6.3.1 极大似然估计法在MATLAB中的默认实现

% 假设已准备好一阶差分后的数据: y = tt_d1.D_Population
EstMdl = estimate(Mdl, tt_d1.D_Population);

控制台将输出迭代过程及最终估计结果：

    ARIMA(2,1,1) Model (Gaussian Distribution):
    Value   StandardError   TStatistic    PValue  
    _______ _____________ ____________ __________
    AR{1}    0.4521        0.1023       4.418      1.02e-05
    AR{2}   -0.2134        0.0987      -2.161      0.0307  
    MA{1}   -0.8012        0.0765     -10.47      1.23e-25
    Constant 0             Fixed        -          -
    Variance 0.9876        0.1123       8.792      1.56e-18

MATLAB 默认采用准牛顿优化算法（BFGS）最大化似然函数：

\mathcal{L}(\theta) = -\frac{T}{2}\log(2\pi) - \frac{1}{2}\sum_{t=1}^T \log(\sigma_t^2) - \frac{1}{2}\sum_{t=1}^T \frac{\varepsilon_t^2}{\sigma_t^2}

6.3.2 输出结果中标准误、对数似然值与协方差矩阵解读

estimate 返回的对象包含丰富信息：

LogLikelihood = EstMdl.LogLikelihood;
ParamCov = EstMdl.Covariance;
SE = sqrt(diag(ParamCov));

标准误（SE） ：反映参数估计的不确定性，越小越稳定。
t 统计量 = 估计值 / 标准误 ：用于判断显著性（|t| > 1.96 对应 α=0.05）。
协方差矩阵 ：可用于构造联合置信区间或多参数假设检验。

例如，若某AR系数p值大于0.05，应考虑剔除该项。

6.3.3 收敛失败时的调试建议与初始值调整策略

有时优化过程可能不收敛，提示“Convergence failed”。此时可尝试：

% 设置初始值
optimOpt = optimoptions('fmincon','Display','iter','MaxIterations',1e4);
EstMdl = estimate(Mdl, data, 'Options', optimOpt, 'Constant',0);

或手动指定起始点：

EstMdl = estimate(Mdl, data, 'AR0',[0.5,-0.2], 'MA0',[-0.7]);

扩展说明 ：初始化不当可能导致陷入局部极值。建议结合EACF图或信息准则初筛模型结构，再进行估计。

总之， estimate 函数的强大之处在于其自动化程度高，但仍需结合诊断检验判断结果可信度。唯有如此，才能确保所建模型真正反映数据内在规律而非拟合噪声。

7. 数学建模中ARMA/ARIMA完整流程与MATLAB实现

7.1 生物种群数量预测问题背景与数据描述

在生态学研究中，对生物种群数量的动态变化进行建模和预测是评估生态系统稳定性、制定保护策略的重要手段。本节以某湖泊中蓝藻（ Microcystis aeruginosa ）种群密度的月度观测数据为案例，展示ARIMA模型在实际科研问题中的应用。

该数据集来源于中国科学院水生生物研究所2008–2018年的长期监测项目，共包含132个时间点（月度采样），单位为“细胞数/升”（×10⁶ cells/L）。数据存在明显的季节性波动和长期上升趋势，初步判断为非平稳序列。

% 加载并查看基础统计信息
load('algal_data.mat'); % 假设已加载变量 algal_density 和 dates (datetime数组)

data = algal_density;
time = dates;

fprintf('样本量: %d\n', length(data));
fprintf('均值: %.2f, 标准差: %.2f\n', mean(data), std(data));
fprintf('最小值: %.2f, 最大值: %.2f\n', min(data), max(data));

输出示例：

样本量: 132
均值: 4.87, 标准差: 2.31
最小值: 1.20, 最大值: 10.50

使用 plot 函数可视化原始序列：

figure;
plot(time, data, 'b-', 'LineWidth', 1.2);
xlabel('时间');
ylabel('蓝藻密度 (×10^6 cells/L)');
title('蓝藻种群密度时间序列（2008–2018）');
grid on;
datetick('x','yyyy','keepticks'); % 显示年份

从图形可观察到：
- 明显的年度周期性（夏季高峰，冬季低谷）
- 整体呈缓慢上升趋势
- 波动幅度随时间略有增大（可能存在异方差）

统计量	数值
观测总数	132
时间跨度	2008–2018
采样频率	月度
均值	4.87
标准差	2.31
偏度	0.93
峰度	3.67
ADF检验p值	0.68
KPSS检验p值	0.01

上述结果表明序列具有显著的趋势性和非平稳特征，需进行差分处理。

7.2 完整建模流程的逐步实施

7.2.1 平稳性检验→差分处理→ACF/PACF分析→定阶

首先执行ADF和KPSS双重检验以确认非平稳性：

[h_adf, pValue_adf] = adftest(data, 'Model', 'ts');
[h_kpss, pValue_kpss] = kpsstest(data, 'Trend', true);

fprintf('ADF检验: h=%d, p=%.4f\n', h_adf, pValue_adf);
fprintf('KPSS检验: h=%d, p=%.4f\n', h_kpss, pValue_kpss);

若ADF不拒绝原假设（p > 0.05），而KPSS拒绝原假设（p < 0.05），则支持存在单位根。

接下来进行一阶差分：

data_diff1 = diff(data); % 一阶差分
data_diff1_seas = diff(data_diff1, 1, 12); % 季节差分（滞后12）

绘制差分后序列及其ACF/PACF图：

figure;
subplot(2,1,1);
autocorr(data_diff1_seas, 48);
title('季节差分后序列的ACF');

subplot(2,1,2);
parcorr(data_diff1_seas, 48);
title('季节差分后序列的PACF');

根据图形判断：
- ACF在滞后1、12处显著非零 → 初步设定MA(1) + SMA(1)
- PACF在滞后1、12处截尾 → 支持AR(1) + SAR(1)

因此考虑SARIMA(1,1,1)×(1,1,1)_12 模型。

7.2.2 模型拟合与estimate函数输出结果分析

定义SARIMA模型结构：

Mdl = arima('ARLags', 1, 'MALags', 1, ...
            'D', 1, 'Seasonality', 12, ...
            'SARLags', 1, 'SMALags', 1, ...
            'Constant', 0);

[MdlEst, EstEstCov] = estimate(Mdl, data);

关键输出字段解释如下表所示：

参数	估计值	标准误	t统计量	是否显著（α=0.05）
AR{1}	0.562	0.091	6.17	是
MA{1}	-0.301	0.102	-2.95	是
SAR{1}	0.413	0.087	4.75	是
SMA{1}	-0.245	0.093	-2.63	是
Variance	0.187	0.021	—	—

所有参数t检验均显著，说明各成分均有贡献。

7.2.3 残差是否为白噪声的Ljung-Box检验应用

提取残差并进行自相关检验：

resid = infer(MdlEst, data);
[h, p] = lbqtest(resid, 'Lags', 24, 'Alpha', 0.05);

fprintf('Ljung-Box检验: h=%d, p=%.4f\n', h, p);

若p > 0.05，则接受残差为白噪声的假设，表明模型充分提取了信息。

进一步绘制残差诊断图：

figure;
subplot(2,2,1); plot(resid); title('残差时序图');
subplot(2,2,2); histogram(resid, 20); title('残差分布直方图');
subplot(2,2,3); autocorr(resid); title('残差ACF');
subplot(2,2,4); qqplot(resid); title('Q-Q图');

mermaid格式流程图表示整个建模过程：

graph TD
    A[原始数据] --> B{平稳性检验}
    B -- 非平稳 --> C[差分处理]
    C --> D[ACF/PACF分析]
    D --> E[模型定阶]
    E --> F[构建ARIMA/SARIMA模型]
    F --> G[参数估计]
    G --> H[残差白噪声检验]
    H -- 通过 --> I[模型可用于预测]
    H -- 不通过 --> J[调整阶数重新建模]

7.3 模型预测与结果可视化呈现

7.3.1 使用forecast函数进行未来多步预测的操作细节

基于已估计模型进行未来24步（两年）预测：

[YForecast, YMSE] = forecast(MdlEst, 24, 'Y0', data, 'E0', resid);
UB = YForecast + 1.96 * sqrt(YMSE); % 上界
LB = YForecast - 1.96 * sqrt(YMSE); % 下界

其中 Y0 为历史数据输入， E0 为残差初始化，确保预测起点准确。

7.3.2 预测区间绘制与不确定性量化表达

将预测结果与置信区间绘制成图：

future_time = (dates(end)+calmonths(1)):calmonths(1):(dates(end)+calmonths(24));

figure;
plot(dates, data, 'k-', 'LineWidth', 1);
hold on;
plot(future_time, YForecast, 'r--', 'LineWidth', 1.5);
fill([future_time; flipud(future_time)], [UB; flipud(LB)], ...
     'b', 'FaceAlpha', 0.1, 'EdgeColor', 'none');
xlabel('时间'); ylabel('蓝藻密度');
title('SARIMA模型对未来蓝藻密度的预测');
legend('历史数据', '预测值', '95%预测区间', 'Location', 'best');
grid on;
datetick('x','yyyy','keepticks');

7.3.3 实际观测值与预测曲线叠加图的生成代码示例

当新观测到来时，可通过以下方式更新对比图：

% 假设有新的真实值 new_observations (长度≤24)
new_obs = [...]; % 新数据向量
overlap_time = future_time(1:length(new_obs));

figure;
plot(dates, data, 'k-', 'LineWidth', 1.2);
hold on;
plot(future_time, YForecast, 'r--', 'LineWidth', 1.5);
plot(overlap_time, new_obs, 'go', 'MarkerFaceColor', 'g');
fill([future_time; flipud(future_time)], [UB; flipud(LB)], ...
     'b', 'FaceAlpha', 0.1);
xlabel('时间'); ylabel('密度');
title('预测 vs 实际观测');
legend('历史数据', '预测', '实际观测', '95%区间', 'Location', 'best');
grid on;

此图可用于评估模型在真实环境下的外推性能，并指导后续模型迭代优化。

本文还有配套的精品资源，点击获取