尧图网站建设 
简介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{tk1}^{T}(x_t - \bar{x})(x_{t-k} - \bar{x})}{\sum_{t1}^{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, nlags20): x np.asarray(x, dtypefloat) 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-k1}$ 对 $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, nlags20): x np.asarray(x, dtypefloat) 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_constantadd) 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参数methodols可以切换到回归法。这里有个值得注意的点Yule-Walker 和 OLS 在小样本下的结果差异不小。模拟一组 $n50$ 的 AR(1) 数据$\phi0.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, lags20, axaxes[0], alpha0.05, titleACF of AR(1)) plot_pacf(x, lags20, axaxes[1], alpha0.05, methodywm, titlePACF of AR(1)) plt.tight_layout() plt.show()参数alpha控制置信区间的显著性水平默认 0.05 对应 95% 区间methodywm是 Yule-Walker 的修正版本在小样本下偏差比原始 YW 更小。改这两个参数图形会明显变化。做研究或线上报告需要给人看的时候建议把alpha0.01或alpha0.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, periods365, freqD) 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, indexdates, namedaily_visits) # ADF 检验p 值 0.05 则不能拒绝单位根序列非平稳 adf_result adfuller(series, autolagAIC) print(fADF p-value: {adf_result[1]:.4f}) # 一阶差分后再检验 diff_series series.diff().dropna() adf_result_diff adfuller(diff_series, autolagAIC) print(fADF p-value after diff: {adf_result_diff[1]:.4f})autolagAIC表示回归滞后阶数由信息准则自动选择默认值也是这个。如果这里 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(nsample500, burnin50) # 画图核心参数 fig, axes plt.subplots(2, 2, figsize(16, 10)) plot_acf(sample_data, lags30, axaxes[0][0], titleACF (default)) plot_pacf(sample_data, lags30, axaxes[0][1], titlePACF (default YW)) plot_acf(sample_data, lags30, axaxes[1][0], titleACF (zeroFalse), zeroFalse) plot_pacf(sample_data, lags30, axaxes[1][1], titlePACF (OLS), methodols) plt.tight_layout() plt.show()参数逐一说明lags30计算并展示 0 到 30 阶的相关性。样本量 500 时30 阶偏大常规建议 $\sqrt{n} \approx 22$ 左右但 ARMA 结构可能隐藏在更高阶多画几阶没坏处。判断时重点看前 5-10 阶即可。zeroFalse跳过滞后 0 阶。滞后 0 阶 ACF 恒为 1画出来会拉伸纵轴把后面有信息量的柱状图压扁。zeroFalse在可视化上往往更实用。methodolsPACF 用回归法估计。小样本下 OLS 给出的标准误可以用于构造置信区间ywm则更快但区间计算依赖正态近似。# 如果环境里还没有这些库按顺序安装 pip install numpy pandas matplotlib statsmodels运行上面的代码后你能看到 ACF 在第 1 阶显著、之后拖尾PACF 在第 1 阶显著、第 2 阶后截尾——这对应 ARMA(1,1) 的典型特征ACF 拖尾而 PACF 在 1 阶截尾。实际数据分析中很少遇到教科书级完美的截尾图噪声和样本截断都会让柱状图在边界附近摆动。3.3 滞后阶数怎么选三个经验法则对比方法公式适用场景固定阈值$\min(10 \log_{10}(n), n-1)$快速预览通用样本量根号$\lfloor \sqrt{n} \rfloor$样本量 100-500 的常规建模信息准则AIC/BIC 自动选择需要严谨定阶时配合 AR 模型滞后阶数选少了会漏掉季节性周期处的尖峰选多了容易把噪声点误读为结构性相关。我的做法是先用根号法则确定一个基线再按数据的最小业务周期比如日数据看 7、14、21 阶额外标注几根竖线快速检查是否存在周期相关。# 在图中标注季节性滞后 lags_baseline int(np.sqrt(len(sample_data))) plot_acf(sample_data, lagsmax(lags_baseline, 20), axaxes[0][0]) for lag in [5, 10, 15, 20]: axes[0][0].axvline(xlag, colorred, linestyle--, alpha0.5)红线和蓝色柱状图的交叉情况能直观判断这个峰值是否可能来自周期。比如日粒度数据在 7 和 14 阶出现显著峰值那就要调整建模方向从普通 ARMA 转向 SARIMA而不是强行提高 AR 阶数去拟合周期效应——那样做会让参数数量爆炸且可解释性几乎为零。4. 从图形特征判断模型阶数拖尾、截尾与实战定阶4.1 四种典型形态对应到 AR/MA/ARMA 的映射表模型ACF 形态PACF 形态判读要点AR(p)拖尾指数衰减或振荡衰减在 p 阶后截尾看 PACF 最有效MA(q)在 q 阶后截尾拖尾指数衰减看 ACF 最有效ARMA(p, q)拖尾拖尾无法直接从图确定 p 和 q白噪声所有阶都不显著所有阶都不显著无建模必要这张表所有时间序列教材都有但实际看图时有个坑样本量不足时即使真实模型是 AR(1)PACF 也可能在第 3 阶出现看似显著的柱——因为 95% 置信区间意味着平均 20 根柱里有一根会误报。滞后总阶数越多误报的绝对概率越大。from statsmodels.tsa.stattools import acf, pacf # 用数值方式输出各阶系数辅助图形判读 acf_vals acf(sample_data, nlags10) pacf_vals pacf(sample_data, nlags10, methodywm) for i, (a, p) in enumerate(zip(acf_vals, pacf_vals)): print(flag {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, keyresults.get) print(fBest (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, lags10, return_dfTrue) print(lb_test) # 残差 ACF 图理论上所有柱都在置信区间内 fig, ax plt.subplots(figsize(10, 4)) plot_acf(residuals, lags20, axax, zeroFalse, titleResidual 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(nsample300, burnin50) group2 np.concatenate([arma_process.generate_sample(nsample150), arma_process.generate_sample(nsample150)]) fig, axes plt.subplots(2, 2, figsize(14, 8)) for idx, data in enumerate([group1, group2]): plot_acf(data, lags20, axaxes[idx][0], zeroFalse, titlefGroup {idx1} ACF) plot_pacf(data, lags20, axaxes[idx][1], zeroFalse, titlefGroup {idx1} 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, indexpd.date_range(2020-01-01, periods120, freqMS) ) fig, axes plt.subplots(2, 2, figsize(14, 8)) # 原始序列的 ACF/PACF plot_acf(monthly, lags30, axaxes[0][0], titleOriginal ACF (seasonal)) plot_pacf(monthly, lags30, axaxes[0][1], titleOriginal PACF (seasonal)) # 季节差分后的 ACF/PACF plot_acf(monthly.diff(12).dropna(), lags30, axaxes[1][0], titleSeasonal Diff ACF, zeroFalse) plot_pacf(monthly.diff(12).dropna(), lags30, axaxes[1][1], titleSeasonal Diff PACF, zeroFalse) plt.tight_layout() plt.show()判断季节性时有个常被忽略的参数lags至少要大于两个完整周期。比如月粒度数据要看到第 24 阶日粒度带周周期的数据要至少看到第 21 阶。如果只取默认的 20 阶部分季节性结构会被截断在视野之外。5.3 三个实践技巧第一个技巧是调整vlines_kwargs和marker参数改善图形可读性。默认的柱状图条幅较宽样本量大时相邻柱子会互相遮挡人眼很难分辨细微的边界情况。plot_acf(series, lags30, axax, zeroFalse, vlines_kwargs{colors: steelblue, linewidth: 1.2}, markero, markersize4)柱状变细、点上加圆点标记之后目光能更快扫出哪些点落在置信区间之外。这个改动不影响任何计算逻辑只是把默认的误差条画法换成点线画法。第二个技巧是用matplotlib的交互模式动态调整置信区间快速排查刚好压线的情况。因为alpha参数会重新计算阴影区域的上下界这是应对边界情况最直接的办法。注意不要只调一个方向——alpha0.05下显著的点在alpha0.01下往往不再显著这提示该点很可能只是噪声。第三个技巧是在建模完成之后单独跑一次 PACF 的残差检验。即使 ACF 图和 Ljung-Box 检验都通过PACF 在某个滞后阶数上仍可能出现孤立峰值。遇到这种情况优先怀疑数据录入错误或者存在缺失值插补的人为痕迹而不是急着修改模型阶数。缺失值用线性插值处理过的序列会在插值区间的边界产生微弱的伪相关这种伪相关在 ACF 上不明显但在 PACF 上会暴露出来。本文还有配套的精品资源点击获取 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{tk1}^{T}(x_t - \bar{x})(x_{t-k} - \bar{x})}{\sum_{t1}^{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, nlags20): x np.asarray(x, dtypefloat) 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-k1}$ 对 $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, nlags20): x np.asarray(x, dtypefloat) 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_constantadd) 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参数methodols可以切换到回归法。这里有个值得注意的点Yule-Walker 和 OLS 在小样本下的结果差异不小。模拟一组 $n50$ 的 AR(1) 数据$\phi0.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, lags20, axaxes[0], alpha0.05, titleACF of AR(1)) plot_pacf(x, lags20, axaxes[1], alpha0.05, methodywm, titlePACF of AR(1)) plt.tight_layout() plt.show()参数alpha控制置信区间的显著性水平默认 0.05 对应 95% 区间methodywm是 Yule-Walker 的修正版本在小样本下偏差比原始 YW 更小。改这两个参数图形会明显变化。做研究或线上报告需要给人看的时候建议把alpha0.01或alpha0.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, periods365, freqD) 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, indexdates, namedaily_visits) # ADF 检验p 值 0.05 则不能拒绝单位根序列非平稳 adf_result adfuller(series, autolagAIC) print(fADF p-value: {adf_result[1]:.4f}) # 一阶差分后再检验 diff_series series.diff().dropna() adf_result_diff adfuller(diff_series, autolagAIC) print(fADF p-value after diff: {adf_result_diff[1]:.4f})autolagAIC表示回归滞后阶数由信息准则自动选择默认值也是这个。如果这里 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(nsample500, burnin50) # 画图核心参数 fig, axes plt.subplots(2, 2, figsize(16, 10)) plot_acf(sample_data, lags30, axaxes[0][0], titleACF (default)) plot_pacf(sample_data, lags30, axaxes[0][1], titlePACF (default YW)) plot_acf(sample_data, lags30, axaxes[1][0], titleACF (zeroFalse), zeroFalse) plot_pacf(sample_data, lags30, axaxes[1][1], titlePACF (OLS), methodols) plt.tight_layout() plt.show()参数逐一说明lags30计算并展示 0 到 30 阶的相关性。样本量 500 时30 阶偏大常规建议 $\sqrt{n} \approx 22$ 左右但 ARMA 结构可能隐藏在更高阶多画几阶没坏处。判断时重点看前 5-10 阶即可。zeroFalse跳过滞后 0 阶。滞后 0 阶 ACF 恒为 1画出来会拉伸纵轴把后面有信息量的柱状图压扁。zeroFalse在可视化上往往更实用。methodolsPACF 用回归法估计。小样本下 OLS 给出的标准误可以用于构造置信区间ywm则更快但区间计算依赖正态近似。# 如果环境里还没有这些库按顺序安装 pip install numpy pandas matplotlib statsmodels运行上面的代码后你能看到 ACF 在第 1 阶显著、之后拖尾PACF 在第 1 阶显著、第 2 阶后截尾——这对应 ARMA(1,1) 的典型特征ACF 拖尾而 PACF 在 1 阶截尾。实际数据分析中很少遇到教科书级完美的截尾图噪声和样本截断都会让柱状图在边界附近摆动。3.3 滞后阶数怎么选三个经验法则对比方法公式适用场景固定阈值$\min(10 \log_{10}(n), n-1)$快速预览通用样本量根号$\lfloor \sqrt{n} \rfloor$样本量 100-500 的常规建模信息准则AIC/BIC 自动选择需要严谨定阶时配合 AR 模型滞后阶数选少了会漏掉季节性周期处的尖峰选多了容易把噪声点误读为结构性相关。我的做法是先用根号法则确定一个基线再按数据的最小业务周期比如日数据看 7、14、21 阶额外标注几根竖线快速检查是否存在周期相关。# 在图中标注季节性滞后 lags_baseline int(np.sqrt(len(sample_data))) plot_acf(sample_data, lagsmax(lags_baseline, 20), axaxes[0][0]) for lag in [5, 10, 15, 20]: axes[0][0].axvline(xlag, colorred, linestyle--, alpha0.5)红线和蓝色柱状图的交叉情况能直观判断这个峰值是否可能来自周期。比如日粒度数据在 7 和 14 阶出现显著峰值那就要调整建模方向从普通 ARMA 转向 SARIMA而不是强行提高 AR 阶数去拟合周期效应——那样做会让参数数量爆炸且可解释性几乎为零。4. 从图形特征判断模型阶数拖尾、截尾与实战定阶4.1 四种典型形态对应到 AR/MA/ARMA 的映射表模型ACF 形态PACF 形态判读要点AR(p)拖尾指数衰减或振荡衰减在 p 阶后截尾看 PACF 最有效MA(q)在 q 阶后截尾拖尾指数衰减看 ACF 最有效ARMA(p, q)拖尾拖尾无法直接从图确定 p 和 q白噪声所有阶都不显著所有阶都不显著无建模必要这张表所有时间序列教材都有但实际看图时有个坑样本量不足时即使真实模型是 AR(1)PACF 也可能在第 3 阶出现看似显著的柱——因为 95% 置信区间意味着平均 20 根柱里有一根会误报。滞后总阶数越多误报的绝对概率越大。from statsmodels.tsa.stattools import acf, pacf # 用数值方式输出各阶系数辅助图形判读 acf_vals acf(sample_data, nlags10) pacf_vals pacf(sample_data, nlags10, methodywm) for i, (a, p) in enumerate(zip(acf_vals, pacf_vals)): print(flag {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, keyresults.get) print(fBest (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, lags10, return_dfTrue) print(lb_test) # 残差 ACF 图理论上所有柱都在置信区间内 fig, ax plt.subplots(figsize(10, 4)) plot_acf(residuals, lags20, axax, zeroFalse, titleResidual 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(nsample300, burnin50) group2 np.concatenate([arma_process.generate_sample(nsample150), arma_process.generate_sample(nsample150)]) fig, axes plt.subplots(2, 2, figsize(14, 8)) for idx, data in enumerate([group1, group2]): plot_acf(data, lags20, axaxes[idx][0], zeroFalse, titlefGroup {idx1} ACF) plot_pacf(data, lags20, axaxes[idx][1], zeroFalse, titlefGroup {idx1} 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, indexpd.date_range(2020-01-01, periods120, freqMS) ) fig, axes plt.subplots(2, 2, figsize(14, 8)) # 原始序列的 ACF/PACF plot_acf(monthly, lags30, axaxes[0][0], titleOriginal ACF (seasonal)) plot_pacf(monthly, lags30, axaxes[0][1], titleOriginal PACF (seasonal)) # 季节差分后的 ACF/PACF plot_acf(monthly.diff(12).dropna(), lags30, axaxes[1][0], titleSeasonal Diff ACF, zeroFalse) plot_pacf(monthly.diff(12).dropna(), lags30, axaxes[1][1], titleSeasonal Diff PACF, zeroFalse) plt.tight_layout() plt.show()判断季节性时有个常被忽略的参数lags至少要大于两个完整周期。比如月粒度数据要看到第 24 阶日粒度带周周期的数据要至少看到第 21 阶。如果只取默认的 20 阶部分季节性结构会被截断在视野之外。5.3 三个实践技巧第一个技巧是调整vlines_kwargs和marker参数改善图形可读性。默认的柱状图条幅较宽样本量大时相邻柱子会互相遮挡人眼很难分辨细微的边界情况。plot_acf(series, lags30, axax, zeroFalse, vlines_kwargs{colors: steelblue, linewidth: 1.2}, markero, markersize4)柱状变细、点上加圆点标记之后目光能更快扫出哪些点落在置信区间之外。这个改动不影响任何计算逻辑只是把默认的误差条画法换成点线画法。第二个技巧是用matplotlib的交互模式动态调整置信区间快速排查刚好压线的情况。因为alpha参数会重新计算阴影区域的上下界这是应对边界情况最直接的办法。注意不要只调一个方向——alpha0.05下显著的点在alpha0.01下往往不再显著这提示该点很可能只是噪声。第三个技巧是在建模完成之后单独跑一次 PACF 的残差检验。即使 ACF 图和 Ljung-Box 检验都通过PACF 在某个滞后阶数上仍可能出现孤立峰值。遇到这种情况优先怀疑数据录入错误或者存在缺失值插补的人为痕迹而不是急着修改模型阶数。缺失值用线性插值处理过的序列会在插值区间的边界产生微弱的伪相关这种伪相关在 ACF 上不明显但在 PACF 上会暴露出来。本文还有配套的精品资源点击获取