Python时间序列分析：ACF与PACF计算逻辑与ARIMA定阶实战

简介：Python实现时间序列自相关图（ACF）与偏自相关图（PACF）的PDF教程，面向数据分析、统计建模及金融经济领域从业者，帮助读者理解时间序列模式并通过Python工具完成可视化。教程从ACF和PACF的基本概念讲起，清晰解释自相关系数与偏自相关系数的含义，结合statsmodels库的plot_acf和plot_pacf函数演示具体代码，并补充seaborn热力图对多变量相关性进行可视化，为ARIMA模型定阶与特征相关性分析提供实操参考。资源共1个PDF文件，大小77KB，内容精炼，适合快速阅读与随查随用。截至目前已有11434人浏览学习，口碑验证了其实用价值。通过学习，读者能快速掌握绘制ACF/PACF图的核心步骤，理解截尾与拖尾特征在模型识别中的作用，并借由相关性热力图辅助探索数据关联，是一份紧凑高效的时间序列可视化入门资料。

1. 从一张图画错模型阶数说起

很多人在做时间序列分析时，把 ACF 和 PACF 当成"画两张图，看一眼拖尾还是截尾"的步骤。但实际落地时你会发现，statsmodels默认参数画出来的图，和教科书上的示意图差距很大：样本量不够时置信区间宽得离谱，非平稳序列的 ACF 拖尾拖到天边，差分阶数没确认就急着看图定阶，最后 ARIMA 模型参数估计出来全是负数。更隐蔽的问题是plot_pacf默认用 Yule-Walker 方法，和小样本下用 OLS 回归估计的结果能差出一倍。

这篇文章讲清楚 ACF 和 PACF 背后的计算逻辑，以及用 Python 实现时从数据准备、平稳性检验、画图参数到定阶判断的完整路径。读者如果是刚接触时间序列的工程师，可以照着一行行跑通；如果已经做过几轮 ARIMA 建模，后面关于置信区间修正、残差白噪声验证和季节性判断的部分，也能补上一些平时容易忽略的细节。

2. ACF 与 PACF 的计算逻辑：为什么不能只看图

2.1 自相关函数到底在算什么

自相关函数度量的是同一个序列在不同滞后阶数下的线性相关性。给定时间序列 $x_1, x_2, ..., x_T$，滞后 $k$ 阶的样本自相关系数定义为：

$$\hat{\rho}k = \frac{\sum{t=k+1}^{T}(x_t - \bar{x})(x_{t-k} - \bar{x})}{\sum_{t=1}^{T}(x_t - \bar{x})^2}$$

分子是 $x_t$ 和 $x_{t-k}$ 的协方差，分母是方差。注意这里用的是同一个序列的均值 $\bar{x}$，而不是分别计算两组数据的均值——因为理论上平稳序列的均值是常数，这个假设是所有后续判断的前提。

实现上通常用两种方式：一种是直接按公式遍历计算，复杂度 $O(T \times k)$；另一种是先做 FFT 变换计算互相关，再把结果归一化。statsmodels底层用的acovf函数在样本量较大的时候会走 FFT 路径，但处理缺失值时退化到逐项计算。

import numpy as np def acf_by_hand(x, nlags=20): x = np.asarray(x, dtype=float) x = x - x.mean() # 中心化，保证零均值 n = len(x) # 先算方差（滞后0阶的自协方差） c0 = np.sum(x ** 2) / n acf_vals = [1.0] # lag 0 的 ACF 恒为 1 for k in range(1, nlags + 1): ck = np.sum(x[k:] * x[:-k]) / n acf_vals.append(ck / c0) return np.array(acf_vals)

这段代码里除以n而不是n - k，是有意为之。教科书里有两种分母的写法，除以 $n$ 得到的自相关矩阵保证正定性，在后续拟合 AR 模型时不会出现奇异矩阵；除以 $n-k$ 是无偏估计但可能破坏正定性。statsmodels默认也是除以n，这点和 R 的acf()函数保持一致。

使用这段手工实现时要注意两个参数：nlags控制计算到多少阶，经验上取 $\min(10 \log_{10}(n), n-1)$ 是一个相对稳妥的选择，样本量 100 时大约算到 20 阶；x.mean()这一步不能省，如果直接拿原始数据算，滞后阶数大的时候分子分母都会被均值项帯偏。

2.2 偏自相关：排除中间变量的干扰

偏自相关函数度量的是在剔除 $x_{t-1}, x_{t-2}, ..., x_{t-k+1}$ 对 $x_t$ 和 $x_{t-k}$ 的影响之后，两者之间剩余的线性关系。AR(1) 过程 $x_t = \phi x_{t-1} + \varepsilon_t$ 的 ACF 在滞后 1 阶之后仍然不为零（拖尾），因为 $x_t$ 通过 $x_{t-1}$ 间接和 $x_{t-2}$ 相关。但 PACF 在滞后 2 阶及以后应该接近零，因为直接关联已经被 1 阶滞后解释了。

PACF 的估计有三种常见路径。第一种是 Yule-Walker 方程，用样本 ACF 值代入：

$$\begin{bmatrix} 1 & \hat{\rho}1 & \cdots & \hat{\rho}{k-1} \ \hat{\rho}1 & 1 & \cdots & \hat{\rho}{k-2} \ \vdots & \vdots & \ddots & \vdots \ \hat{\rho}{k-1} & \hat{\rho}{k-2} & \cdots & 1 \end{bmatrix} \begin{bmatrix} \phi_{k1} \ \phi_{k2} \ \vdots \ \phi_{kk} \end{bmatrix} = \begin{bmatrix} \hat{\rho}_1 \ \hat{\rho}_2 \ \vdots \ \hat{\rho}_k \end{bmatrix}$$

解出 $\phi_{kk}$ 就是滞后 $k$ 阶的 PACF 值。第二种是 OLS 回归：把 $x_t$ 对 $x_{t-1}, ..., x_{t-k}$ 做回归，最后一个回归系数的估计值就是 PACF。第三种是 Levinson-Durbin 递推，利用 Toeplitz 矩阵结构把复杂度压到 $O(k^2)$。

from statsmodels.regression.linear_model import OLS from statsmodels.tools.tools import add_constant def pacf_via_ols(x, nlags=20): x = np.asarray(x, dtype=float) n = len(x) pacf_vals = [1.0] # 从滞后0开始逐个构造滞后矩阵并回归 for k in range(1, nlags + 1): # 滞后矩阵：第 t 行是 [x_{t-1}, x_{t-2}, ..., x_{t-k}] X = np.column_stack([x[k - i - 1:n - i - 1] for i in range(k)]) y = x[k:] # 加截距项后拟合 X_design = add_constant(X, has_constant='add') model = OLS(y, X_design).fit() pacf_vals.append(model.params[-1]) # 最后一个系数为 k 阶 PACF return np.array(pacf_vals)

OLS 方法的优势在于能同时拿到系数的标准误，对后续判断"PACF 是否显著非零"有帮助。但它对样本量敏感，当 $k$ 接近 $n/10$ 时，设计矩阵接近奇异，回归系数方差暴涨。实际操作中statsmodels的plot_pacf默认用的就是 Yule-Walker，参数method='ols'可以切换到回归法。

这里有个值得注意的点：Yule-Walker 和 OLS 在小样本下的结果差异不小。模拟一组 $n=50$ 的 AR(1) 数据，$\phi=0.7$，滞后 5 阶的 PACF 用 Yule-Walker 估计约在 0.03 附近，OLS 可能到 0.08——都在置信区间内，但数值分布完全不同。别因为两张图长得不一样就怀疑代码写错了。

2.3 置信区间：判断显著性的基准线

plot_acf和plot_pacf画出的阴影区域代表的是"若真实自相关为零，样本估计值的抽样分布"的近似置信区间。默认用的是 $\pm 1.96 / \sqrt{n}$，即正态近似下 95% 的置信带。

这个公式隐含的假设是：样本量足够大，且序列本身是白噪声。当你分析的是残差序列时这个假设基本成立；但分析原始序列时，ACF 的各阶估计值之间存在强相关，使用统一带宽会低估联合显著性——多个点同时在区间外，未必代表真的显著。

import numpy as np import matplotlib.pyplot as plt from statsmodels.graphics.tsaplots import plot_acf, plot_pacf # 模拟一个 AR(1) 过程验证置信区间 np.random.seed(42) n = 200 phi = 0.75 x = np.zeros(n) for t in range(1, n): x[t] = phi * x[t-1] + np.random.randn() fig, axes = plt.subplots(1, 2, figsize=(14, 5)) plot_acf(x, lags=20, ax=axes[0], alpha=0.05, title='ACF of AR(1)') plot_pacf(x, lags=20, ax=axes[1], alpha=0.05, method='ywm', title='PACF of AR(1)') plt.tight_layout() plt.show()

参数alpha控制置信区间的显著性水平，默认 0.05 对应 95% 区间；method='ywm'是 Yule-Walker 的修正版本，在小样本下偏差比原始 YW 更小。改这两个参数，图形会明显变化。做研究或线上报告需要给人看的时候，建议把alpha=0.01或alpha=0.1的结果都跑一遍，避免只展示一种设定下"恰好在边界"的图形。

3. 用 Python 画 ACF 和 PACF 的完整流程

3.1 数据准备：先做平稳性检验再画图

画 ACF/PACF 之前必须确认序列平稳——非平稳序列的自相关函数会以极慢的速度衰减，画出来的图往往是一大片滞后阶数都在置信区间外，看不出任何结构。实际工作流通常是先可视化原始序列，再跑 ADF 检验，必要时做差分，最后才进入相关图环节。

import pandas as pd from statsmodels.tsa.stattools import adfuller # 以某电商平台的日访问量数据为例（模拟结构） dates = pd.date_range('2024-01-01', periods=365, freq='D') trend = np.linspace(50, 120, 365) seasonal = 10 * np.sin(2 * np.pi * np.arange(365) / 7) noise = np.random.randn(365) * 3 series = pd.Series(trend + seasonal + noise, index=dates, name='daily_visits') # ADF 检验：p 值 > 0.05 则不能拒绝单位根，序列非平稳 adf_result = adfuller(series, autolag='AIC') print(f'ADF p-value: {adf_result[1]:.4f}') # 一阶差分后再检验 diff_series = series.diff().dropna() adf_result_diff = adfuller(diff_series, autolag='AIC') print(f'ADF p-value after diff: {adf_result_diff[1]:.4f}')

autolag='AIC'表示回归滞后阶数由信息准则自动选择，默认值也是这个。如果这里 p 值仍然大于 0.05，就需要考虑二阶差分或对数变换。有几个容易踩的坑：一是dropna()之后索引会断，画图没问题但后续建模要注意频率；二是diff()默认是一阶差分，对周期性数据应该先做季节差分，频率为 24 的小时数据做series.diff(24)。

序列平稳之后，ACF 图才能暴露真实结构。一个典型的现象是：趋势数据的一阶差分后 ACF 在滞后 1 阶显示明显负值，这是过度差分的信号，后面会展开讲。

3.2 核心代码：plot_acf 与 plot_pacf 的最小可运行示例

import pandas as pd import numpy as np import matplotlib.pyplot as plt from statsmodels.graphics.tsaplots import plot_acf, plot_pacf from statsmodels.tsa.arima_process import ArmaProcess # 构造一个 ARMA(1,1) 过程作为示例数据 ar = np.array([1, -0.6]) # AR 系数：x_t = 0.6*x_{t-1} + eps_t ma = np.array([1, 0.4]) # MA 系数：加上 0.4*eps_{t-1} arma_process = ArmaProcess(ar, ma) sample_data = arma_process.generate_sample(nsample=500, burnin=50) # 画图核心参数 fig, axes = plt.subplots(2, 2, figsize=(16, 10)) plot_acf(sample_data, lags=30, ax=axes[0][0], title='ACF (default)') plot_pacf(sample_data, lags=30, ax=axes[0][1], title='PACF (default YW)') plot_acf(sample_data, lags=30, ax=axes[1][0], title='ACF (zero=False)', zero=False) plot_pacf(sample_data, lags=30, ax=axes[1][1], title='PACF (OLS)', method='ols') plt.tight_layout() plt.show()

参数逐一说明：

• lags=30：计算并展示 0 到 30 阶的相关性。样本量 500 时，30 阶偏大，常规建议 $\sqrt{n} \approx 22$ 左右；但 ARMA 结构可能隐藏在更高阶，多画几阶没坏处。判断时重点看前 5-10 阶即可。
• zero=False：跳过滞后 0 阶。滞后 0 阶 ACF 恒为 1，画出来会拉伸纵轴，把后面有信息量的柱状图压扁。zero=False在可视化上往往更实用。
• method='ols'：PACF 用回归法估计。小样本下 OLS 给出的标准误可以用于构造置信区间，ywm则更快但区间计算依赖正态近似。

# 如果环境里还没有这些库，按顺序安装 pip install numpy pandas matplotlib statsmodels

运行上面的代码后，你能看到 ACF 在第 1 阶显著、之后拖尾，PACF 在第 1 阶显著、第 2 阶后截尾——这对应 ARMA(1,1) 的典型特征：ACF 拖尾而 PACF 在 1 阶截尾。实际数据分析中，很少遇到教科书级完美的截尾图，噪声和样本截断都会让柱状图在边界附近摆动。

3.3 滞后阶数怎么选：三个经验法则对比

滞后阶数选少了会漏掉季节性周期处的尖峰，选多了容易把噪声点误读为结构性相关。我的做法是：先用根号法则确定一个基线，再按数据的最小业务周期（比如日数据看 7、14、21 阶）额外标注几根竖线，快速检查是否存在周期相关。

# 在图中标注季节性滞后 lags_baseline = int(np.sqrt(len(sample_data))) plot_acf(sample_data, lags=max(lags_baseline, 20), ax=axes[0][0]) for lag in [5, 10, 15, 20]: axes[0][0].axvline(x=lag, color='red', linestyle='--', alpha=0.5)

红线和蓝色柱状图的交叉情况能直观判断"这个峰值是否可能来自周期"。比如日粒度数据在 7 和 14 阶出现显著峰值，那就要调整建模方向，从普通 ARMA 转向 SARIMA，而不是强行提高 AR 阶数去拟合周期效应——那样做会让参数数量爆炸且可解释性几乎为零。

4. 从图形特征判断模型阶数：拖尾、截尾与实战定阶

4.1 四种典型形态对应到 AR/MA/ARMA 的映射表

这张表所有时间序列教材都有，但实际看图时有个坑：样本量不足时，即使真实模型是 AR(1)，PACF 也可能在第 3 阶出现"看似显著"的柱——因为 95% 置信区间意味着平均 20 根柱里有一根会误报。滞后总阶数越多，误报的绝对概率越大。

from statsmodels.tsa.stattools import acf, pacf # 用数值方式输出各阶系数，辅助图形判读 acf_vals = acf(sample_data, nlags=10) pacf_vals = pacf(sample_data, nlags=10, method='ywm') for i, (a, p) in enumerate(zip(acf_vals, pacf_vals)): print(f'lag {i}: ACF={a:.3f}, PACF={p:.3f}')

把数值打印出来看曲线摆动幅度比只盯着图更准确——图形渲染时纵轴的范围会自动缩放，某些小幅波动在图上看起来"很大"，数值上可能只有 0.02，完全在噪声范围内。

4.2 定阶不只看显著性：AIC/BIC 交叉验证

图形提供的只是候选阶数范围，最终需要信息准则来裁决。

from statsmodels.tsa.arima.model import ARIMA # 图形建议可能是 ARMA(1,1) 或 AR(1)，用 AIC 对候选模型逐个打分 results = {} for p in range(0, 4): for q in range(0, 4): try: model = ARIMA(sample_data, order=(p, 0, q)).fit() results[(p, q)] = model.aic except Exception: continue best = min(results, key=results.get) print(f'Best (p, q) by AIC: {best}, AIC={results[best]:.2f}')

图形定阶和信息准则冲突时怎么办？我的经验是：如果候选模型之间 AIC 差距小于 2，选阶数更少的那一个（简约原则）；差距在 2-7 之间，选 AIC 更小的；差距大于 7，基本可以判定更复杂模型更优。同时要检查残差是否白噪声——这一步经常被跳过，导致整个定阶链条在最后一环脱节。

4.3 残差白噪声验证：ACF 图的最终用途

拟合出模型后，残差的 ACF 图是最直接的质量检查手段。如果残差序列还有显著的自相关，说明信息没有被完全提取，模型设定有遗漏。

# 残差白噪声检验：Ljung-Box Q 统计量 from statsmodels.stats.diagnostic import acorr_ljungbox residuals = model.resid lb_test = acorr_ljungbox(residuals, lags=10, return_df=True) print(lb_test) # 残差 ACF 图：理论上所有柱都在置信区间内 fig, ax = plt.subplots(figsize=(10, 4)) plot_acf(residuals, lags=20, ax=ax, zero=False, title='Residual ACF') plt.show()

acorr_ljungbox返回的 p 值都大于 0.05 时，残差可以近似视为白噪声。一个细节：lags不能太大，当 $lags$ 接近样本量的 10% 时检验功效会退化。如果模型残差 ACF 在第 1 阶显著为负，多半是过度差分导致的；如果在某个商业周期的阶数上显著正相关，就要重新考虑季节性建模。

5. 进阶用法：多序列对比画布、季节滞后标注与常见坑

5.1 多序列 ACF 对比：判断两组数据是否同构

在对比两个模型的残差结构，或者判断不同分组的序列是否共享相同的 ARMA 结构时，并排画多个 ACF 比逐个画再靠记忆对比可靠得多。

# 假设 group1 和 group2 是两个分组的聚合指标 group1 = arma_process.generate_sample(nsample=300, burnin=50) group2 = np.concatenate([arma_process.generate_sample(nsample=150), arma_process.generate_sample(nsample=150)]) fig, axes = plt.subplots(2, 2, figsize=(14, 8)) for idx, data in enumerate([group1, group2]): plot_acf(data, lags=20, ax=axes[idx][0], zero=False, title=f'Group {idx+1} ACF') plot_pacf(data, lags=20, ax=axes[idx][1], zero=False, title=f'Group {idx+1} PACF') plt.tight_layout() plt.show()

注意group2的构造方式：前 150 个点和后 150 个点各自独立生成再拼接，会导致连接处出现虚假的相关性。真实业务里类似情况会出现在"两个销售渠道数据简单合并"或"跨年数据未做季节调整"的场景。观察这种拼接序列的 ACF 图，你会发现滞后 1 阶的峰值被明显拉高，但滞后 2 阶以上仍然正常——这是结构性断点的典型信号。

5.2 季节性数据的 ACF 特征

周期性数据的 ACF 不会在第 7 阶突然截尾，而是在滞后 7, 14, 21 阶持续出现峰值，峰值大小随滞后阶数缓慢衰减。这种模式无法通过提高 AR 阶数来拟合，正确方向是做季节差分后画图对比。

# 月粒度含季节性的模拟数据 monthly = pd.Series( 100 + 20 * np.sin(2 * np.pi * np.arange(120) / 12) + np.random.randn(120) * 2, index=pd.date_range('2020-01-01', periods=120, freq='MS') ) fig, axes = plt.subplots(2, 2, figsize=(14, 8)) # 原始序列的 ACF/PACF plot_acf(monthly, lags=30, ax=axes[0][0], title='Original ACF (seasonal)') plot_pacf(monthly, lags=30, ax=axes[0][1], title='Original PACF (seasonal)') # 季节差分后的 ACF/PACF plot_acf(monthly.diff(12).dropna(), lags=30, ax=axes[1][0], title='Seasonal Diff ACF', zero=False) plot_pacf(monthly.diff(12).dropna(), lags=30, ax=axes[1][1], title='Seasonal Diff PACF', zero=False) plt.tight_layout() plt.show()

判断季节性时有个常被忽略的参数：lags至少要大于两个完整周期。比如月粒度数据要看到第 24 阶，日粒度带周周期的数据要至少看到第 21 阶。如果只取默认的 20 阶，部分季节性结构会被截断在视野之外。

5.3 三个实践技巧

第一个技巧是调整vlines_kwargs和marker参数改善图形可读性。默认的柱状图条幅较宽，样本量大时相邻柱子会互相遮挡，人眼很难分辨细微的边界情况。

plot_acf(series, lags=30, ax=ax, zero=False, vlines_kwargs={'colors': 'steelblue', 'linewidth': 1.2}, marker='o', markersize=4)

柱状变细、点上加圆点标记之后，目光能更快扫出哪些点落在置信区间之外。这个改动不影响任何计算逻辑，只是把默认的误差条画法换成点线画法。

第二个技巧是用matplotlib的交互模式动态调整置信区间，快速排查"刚好压线"的情况。因为alpha参数会重新计算阴影区域的上下界，这是应对边界情况最直接的办法。注意不要只调一个方向——alpha=0.05下显著的点，在alpha=0.01下往往不再显著，这提示该点很可能只是噪声。

第三个技巧是在建模完成之后，单独跑一次 PACF 的残差检验。即使 ACF 图和 Ljung-Box 检验都通过，PACF 在某个滞后阶数上仍可能出现孤立峰值。遇到这种情况优先怀疑数据录入错误或者存在缺失值插补的人为痕迹，而不是急着修改模型阶数。缺失值用线性插值处理过的序列，会在插值区间的边界产生微弱的伪相关，这种伪相关在 ACF 上不明显，但在 PACF 上会暴露出来。

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