Python实战：用ACF和PACF分析时间序列数据（附完整代码）

最近在做一个销售预测项目时，我遇到了一个典型问题：面对月度销售额数据，如何判断它适合用哪种时间序列模型？是简单的移动平均，还是复杂的ARIMA？团队里有人凭感觉选了滞后3期，结果模型效果时好时坏。这让我意识到，很多数据分析师和量化研究员虽然知道ACF（自相关函数）和PACF（偏自相关函数）这两个名词，但在实际项目中，面对一张张由statsmodels生成的图表，往往感到无从下手——哪些滞后阶数是真正重要的？置信区间那条线到底该怎么看？如何从这些统计图形里，提炼出对模型选择有直接指导意义的结论？

这篇文章，就是为你解决这些实操痛点而写的。我不会过多纠缠于公式推导（那是教科书的事），而是聚焦于如何用Python工具，从真实数据中计算出ACF和PACF，并像一位经验丰富的数据侦探一样，解读图形背后的故事，最终为模型定阶做出可靠决策。无论你是金融领域的量化研究员，还是电商、物流等行业的数据分析师，只要你需要预测未来，这篇文章中的代码和思路都能直接套用。

1. 理解核心概念：超越公式的直觉

在打开Jupyter Notebook写代码之前，我们需要建立对ACF和PACF的“手感”。你可以暂时忘掉那些求和符号，用更直观的方式来理解它们。

想象你正在观察一条河流每日的水位记录。自相关函数（ACF） 回答的问题是：“今天的水位，与昨天、前天、乃至一周前的水位，有多相似？” 它衡量的是时间序列与其自身滞后版本之间的总体相关性。如果今天的值高度依赖于昨天的值，那么滞后1阶的ACF值就会接近1或-1。ACF图能帮助我们判断序列是否具有趋势或季节性。例如，一个具有长期上升趋势的序列，其ACF值会缓慢衰减，因为很久以前的值与当前值仍然有某种关联。

而偏自相关函数（PACF） 则更“挑剔”一些。它问的是：“在已经知道了昨天、前天…直到k-1天前的水位信息后，今天的水位与k天前的水位还有额外的、直接的相关性吗？” 它试图剥离中间滞后项的影响，找出最纯粹的直接关系。这在识别自回归（AR）模型的阶数时至关重要。一个AR(p)模型，本质上就是说，当前值只直接依赖于前p个历史值，PACF图能清晰地指出这个p值大概在哪里。

为了让你一目了然地抓住关键区别，我整理了下面这个对比表格：

特性	自相关函数 (ACF)	偏自相关函数 (PACF)
核心问题	当前值与过去值总体相关性如何？	在控制中间滞后项后，当前值与某一特定滞后值直接相关性如何？
主要用途	识别移动平均（MA）模型的阶数(q)、检测趋势与季节性	识别自回归（AR）模型的阶数(p)
图形特征	可能缓慢衰减或呈现周期性波动	通常在滞后p阶后出现截尾（急剧降至接近0）
在ARIMA中的作用	帮助确定MA(q) 部分	帮助确定AR(p) 部分

提示：很多初学者容易混淆两者的角色。一个简单的记忆法是：PACF看“头”（AR的p），ACF看“尾”（MA的q）。PACF在滞后p阶后突然“砍头”式地截断，ACF则在滞后q阶后“拖尾”式地逐渐衰减至零。

2. 环境搭建与数据准备

工欲善其事，必先利其器。我们首先确保有一个干净、可复现的Python分析环境。我强烈建议使用conda或venv创建独立的虚拟环境，避免包版本冲突。以下是核心依赖库及其作用：

# 创建并激活虚拟环境（以conda为例）
conda create -n timeseries-analysis python=3.9
conda activate timeseries-analysis

# 安装核心数据分析库
pip install numpy pandas matplotlib seaborn

# 安装时间序列分析专用库
pip install statsmodels scikit-learn

• statsmodels： 本文的绝对主角，提供了计算和绘制ACF/PACF的成熟函数（plot_acf, plot_pacf）。
• pandas： 用于数据读取、清洗和结构化，其DataFrame和Series是处理时间序列的绝佳容器。
• matplotlib & seaborn： 可视化黄金组合，用于定制化图表，让结果更美观。

接下来，我们需要一份有代表性的时间序列数据。为了贴近真实场景，我们不使用过于完美的模拟数据，而是从公开数据集中选取。这里我使用statsmodels自带的“美国航空乘客数量”数据集，它包含了1949年至1960年的月度数据，兼具趋势性和明显的季节性，非常适合演示。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.datasets import get_rdataset

# 设置绘图风格，让图表更专业
plt.style.use('seaborn-v0_8-darkgrid')
sns.set_palette("husl")

# 加载经典数据集：AirPassengers
data = get_rdataset('AirPassengers')
df = data.data
df['date'] = pd.to_datetime(df['time'].astype(str) + '-' + (df['time']%1 * 12 + 1).astype(int).astype(str).str.zfill(2) + '-01')
df.set_index('date', inplace=True)
ts_series = df['value']

# 快速查看数据
print(f"数据时间范围: {ts_series.index.min()} 到 {ts_series.index.max()}")
print(f"数据形状: {ts_series.shape}")
print(ts_series.head())

# 绘制原始序列图
fig, ax = plt.subplots(figsize=(12, 5))
ts_series.plot(ax=ax, linewidth=2)
ax.set_title('美国航空月度乘客数量 (1949-1960)', fontsize=14, fontweight='bold')
ax.set_xlabel('日期')
ax.set_ylabel('乘客数 (千)')
plt.tight_layout()
plt.show()

运行这段代码，你会看到一条典型的非平稳时间序列曲线：长期上升趋势和每年冬季的周期性高峰。我们的分析将基于这个数据集展开。

3. 计算与可视化：statsmodels实战详解

现在进入核心环节。statsmodels.tsa.stattools模块提供了acf和pacf函数用于计算，而statsmodels.graphics.tsaplots中的plot_acf和plot_pacf则是快速可视化的利器。但直接调用绘图函数得到的结果往往需要进一步解读。

3.1 绘制基础ACF与PACF图

让我们先生成最基础的图形，看看数据“长什么样”。

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

fig, axes = plt.subplots(2, 1, figsize=(12, 8))

# 绘制ACF图，设置40阶滞后
plot_acf(ts_series, lags=40, ax=axes[0], alpha=0.05)
axes[0].set_title('自相关函数 (ACF) 图', fontsize=13, fontweight='bold')
axes[0].set_ylabel('自相关系数')

# 绘制PACF图
plot_pacf(ts_series, lags=40, ax=axes[1], alpha=0.05, method='ols') # 使用OLS回归法计算PACF
axes[1].set_title('偏自相关函数 (PACF) 图', fontsize=13, fontweight='bold')
axes[1].set_ylabel('偏自相关系数')
axes[1].set_xlabel('滞后阶数')

plt.tight_layout()
plt.show()

生成的图形中，你会看到：

• 每个柱状图代表对应滞后阶数的相关系数值。
• 蓝色阴影区域是95%的置信区间。通常，完全在阴影区域内的柱状图被认为在统计上不显著（即相关系数与0无显著差异）。
• ACF图呈现缓慢衰减且具有明显的周期性（大约每12个月一个高峰），这是趋势+季节性的典型标志。
• PACF图在滞后1阶和13阶处有显著峰值，之后基本落在置信区间内。

注意：plot_pacf函数中的method参数默认为'ywm'（Yule-Walker方程，带偏差修正），对于较长的序列表现良好。'ols'（普通最小二乘法）是另一种常用方法，尤其在样本量较大时更为准确。如果结果差异不大，任选其一即可；如果差异显著，可以尝试两种方法并参考更稳健的那个。

3.2 深入解读图形：从模式到决策

看图说话是门技术活。面对上面的图形，我们可以系统地拆解出以下信息：

1. 平稳性判断：一个平稳时间序列的ACF会相对快速地衰减至零。我们的ACF图衰减非常慢，这是非平稳的强信号，意味着数据需要进行差分处理。
2. 季节性识别：ACF在滞后12、24、36阶处出现显著高峰，这明确指出了12个月的季节性周期。这在月度数据中非常常见。
3. 模型阶数初探：
   • AR(p)阶数：观察PACF图。它在滞后1阶处有一个非常显著的尖峰，然后在滞后13阶处有另一个较小但超出置信区间的尖峰。滞后1阶的显著性表明可能存在AR(1) 成分。滞后13阶的显著性可能暗示着季节性自回归效应（即当前月与去年同月的直接关系）。一个常见的初步判断是 p=1。
   • MA(q)阶数：观察ACF图。它没有在某个滞后阶数后突然截断，而是拖尾衰减。仅从经典理论看，这暗示MA阶数可能为0。但请注意，对于非平稳数据，ACF的拖尾是常态，必须先做差分使其平稳后，再重新分析ACF。

基于此，我们下一步的行动路线就很清晰了：对原序列进行一阶差分和季节性差分，以消除趋势和季节性，得到一个平稳序列，然后再对其绘制ACF/PACF图。

# 进行一阶非季节性差分和一阶季节性差分（周期s=12）
ts_diff = ts_series.diff().dropna() # 一阶差分去趋势
ts_diff_seasonal = ts_diff.diff(periods=12).dropna() # 对差分后的序列再做12步差分去季节性

# 绘制差分后的序列
fig, axes = plt.subplots(3, 1, figsize=(14, 10))
ts_series.plot(ax=axes[0], title='原始序列')
ts_diff.plot(ax=axes[1], title='一阶差分后序列 (去趋势)')
ts_diff_seasonal.plot(ax=axes[2], title='一阶差分+季节性差分后序列 (去趋势&季节性)')
plt.tight_layout()
plt.show()

# 对平稳化后的序列重新绘制ACF/PACF
fig, axes = plt.subplots(2, 1, figsize=(12, 8))
plot_acf(ts_diff_seasonal, lags=40, ax=axes[0], alpha=0.05, title='平稳序列的ACF图')
plot_pacf(ts_diff_seasonal, lags=40, ax=axes[1], alpha=0.05, method='ols', title='平稳序列的PACF图')
plt.tight_layout()
plt.show()

观察新的ACF/PACF图，你会发现：

• ACF在滞后1阶和12阶处有显著的负相关，之后基本在置信区间内波动。这提示我们可能需要对差分后的序列考虑 MA(1) 和季节性MA(1) 成分。
• PACF在滞后1、2、3阶以及12、24阶处有若干超出置信区间的值，呈现拖尾衰减，这符合AR模型的特征，但不如ACF的截尾特征明显。

综合来看，对于原始序列AirPassengers，一个合理的ARIMA模型结构猜想是：(p,d,q) = (1,1,0) 或 (0,1,1)，并需要加上季节性部分 (P,D,Q,s) = (0,1,1,12) 或 (1,1,0,12)。这正是Box-Jenkins方法的核心：通过观察差分后平稳序列的ACF和PACF图，来初步识别ARIMA模型的(p,d,q)参数。

4. 高级技巧与避坑指南

掌握了基础操作后，我们来看看在实际项目中可能遇到的复杂情况和解决方案。

4.1 处理置信区间的陷阱

默认的95%置信区间（alpha=0.05）是基于一个假设：序列是白噪声。但在很多真实数据中，尤其是金融时间序列，存在波动率聚集现象，这会导致ACF/PACF的估计方差变大，使得默认的置信区间过窄，从而可能将本不显著的滞后误判为显著。

一个更稳健的做法是使用Ljung-Box检验或Bartlett公式计算的置信区间，但statsmodels的绘图函数默认不提供。我们可以手动计算并叠加到图上，作为交叉验证。

from statsmodels.stats.diagnostic import acorr_ljungbox

# 计算Ljung-Box检验的p值
lb_test = acorr_ljungbox(ts_diff_seasonal, lags=[10, 20, 30], return_df=True)
print("Ljung-Box检验结果 (检验差分后序列是否为白噪声):")
print(lb_test)

# 如果p值很小（如<0.05），则拒绝“是白噪声”的原假设，说明序列仍有自相关。
# 这时，仅凭图形判断阶数就需要更加谨慎，可能需要结合AIC/BIC等信息准则。

4.2 样本量不足时的策略

当你的时间序列数据点很少（比如少于50个）时，ACF/PACF的估计会非常不可靠，图形可能显得杂乱无章。此时，盲目依赖图形定阶风险极高。你可以采取以下策略：

• 减少滞后阶数：将lags参数设置为小于n/4的值，例如 lags=min(10, len(series)//4)。
• 优先使用信息准则：使用ARIMA模型的auto_arima函数（来自pmdarima库）或通过网格搜索，选择使AIC或BIC最小的(p,q)组合。
• 结合业务理解：例如，在销售数据中，滞后1（上月）、滞后12（去年同月）通常是最有意义的。可以优先考察这些特定滞后阶的显著性。

4.3 自动化探索与结果输出

在需要批量分析多个时间序列（如不同产品的销量、不同城市的指标）时，手动看图效率太低。我们可以编写一个函数，自动计算关键特征并生成报告。

def analyze_ts_autocorr(series, name="序列", max_lags=40):
    """
    自动化分析时间序列的自相关特征并生成文本报告。
    """
    from statsmodels.tsa.stattools import acf, pacf
    import textwrap

    # 计算ACF和PACF值及置信区间
    acf_vals, acf_confint = acf(series, nlags=max_lags, alpha=0.05, fft=False)
    pacf_vals, pacf_confint = pacf(series, nlags=max_lags, alpha=0.05)

    # 找出显著不为0的滞后阶 (绝对值超出置信区间)
    sig_acf_lags = np.where(np.abs(acf_vals[1:]) > np.abs(acf_confint[1:, 0] - acf_vals[1:]))[0] + 1
    sig_pacf_lags = np.where(np.abs(pacf_vals[1:]) > np.abs(pacf_confint[1:, 0] - pacf_vals[1:]))[0] + 1

    # 生成分析报告
    report = f"""
    {'='*60}
    时间序列分析报告: {name}
    样本量: {len(series)}
    {'='*60}
    1. 自相关 (ACF) 分析:
        - 显著滞后阶: {sig_acf_lags if len(sig_acf_lags)>0 else '无'}
        - 首个显著滞后后的模式: {'截尾' if len(sig_acf_lags)>0 and max(sig_acf_lags)<10 else '拖尾'}
    2. 偏自相关 (PACF) 分析:
        - 显著滞后阶: {sig_pacf_lags if len(sig_pacf_lags)>0 else '无'}
        - 首个显著滞后后的模式: {'截尾' if len(sig_pacf_lags)>0 and max(sig_pacf_lags)<10 else '拖尾'}
    3. 初步建模建议 (基于图形法):
        - 若ACF拖尾、PACF在p阶后截尾 -> 考虑 AR(p) 模型。
        - 若ACF在q阶后截尾、PACF拖尾 -> 考虑 MA(q) 模型。
        - 若两者均拖尾 -> 考虑 ARMA(p,q) 或 ARIMA(p,d,q) 模型。
        - 注意季节性模式 (如滞后12,24阶显著)。
    **注意**: 此建议仅为初步参考，务必结合差分后序列的ACF/PACF及信息准则最终确定模型。
    {'='*60}
    """
    print(textwrap.dedent(report))
    return sig_acf_lags, sig_pacf_lags

# 对差分后的平稳序列使用该函数
sig_acf, sig_pacf = analyze_ts_autocorr(ts_diff_seasonal, name="差分平稳化序列")

这个函数会输出一个简洁的报告，直接指出显著的滞后阶数，并给出初步的模型类型建议，极大提升了批量分析的效率。

5. 从分析到建模：完整工作流示例

最后，我们将整个流程串联起来，从原始数据出发，通过ACF/PACF分析初步确定模型结构，然后建立一个简单的ARIMA模型，并评估其效果。这形成了一个从诊断到建模的闭环。

from statsmodels.tsa.arima.model import ARIMA
from sklearn.metrics import mean_absolute_error, mean_squared_error

# 基于前面的分析，我们尝试一个季节性ARIMA模型: (0,1,1) x (0,1,1,12)
# 注意：这里仅为示例，最优参数需通过更系统的网格搜索确定。
order = (0, 1, 1) # 非季节性部分 (p,d,q)
seasonal_order = (0, 1, 1, 12) # 季节性部分 (P,D,Q,s)

# 划分训练集和测试集 (最后24个月作为测试)
train_size = len(ts_series) - 24
train, test = ts_series.iloc[:train_size], ts_series.iloc[train_size:]

# 拟合模型
model = ARIMA(train, order=order, seasonal_order=seasonal_order)
model_fit = model.fit()
print(model_fit.summary())

# 进行预测
forecast_steps = 24
forecast_result = model_fit.get_forecast(steps=forecast_steps)
forecast_mean = forecast_result.predicted_mean
forecast_conf_int = forecast_result.conf_int()

# 计算预测误差
mae = mean_absolute_error(test, forecast_mean)
rmse = np.sqrt(mean_squared_error(test, forecast_mean))
print(f"\n预测性能评估 (测试集最后24个月):")
print(f"平均绝对误差 (MAE): {mae:.2f}")
print(f"均方根误差 (RMSE): {rmse:.2f}")

# 可视化预测结果
fig, ax = plt.subplots(figsize=(14, 7))
train.iloc[-60:].plot(ax=ax, label='训练集 (后60期)', linewidth=2)
test.plot(ax=ax, label='测试集 (真实值)', linewidth=2)
forecast_mean.plot(ax=ax, label='预测值', style='--', linewidth=2)
ax.fill_between(forecast_conf_int.index,
                forecast_conf_int.iloc[:, 0],
                forecast_conf_int.iloc[:, 1], color='gray', alpha=0.2, label='95% 置信区间')
ax.set_title('季节性ARIMA模型预测效果', fontsize=14, fontweight='bold')
ax.set_xlabel('日期')
ax.set_ylabel('乘客数 (千)')
ax.legend()
plt.tight_layout()
plt.show()

# 最后，检查模型的残差是否符合白噪声假设（理想的模型其残差应为白噪声）
residuals = model_fit.resid
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
plot_acf(residuals, lags=40, ax=axes[0], alpha=0.05, title='模型残差ACF图')
plot_pacf(residuals, lags=40, ax=axes[1], alpha=0.05, title='模型残差PACF图')
plt.tight_layout()
plt.show()

如果残差的ACF/PACF图显示没有显著的自相关（所有柱状图基本都在置信区间内），说明模型已经较好地捕捉了数据中的规律，残差接近随机噪声，这是一个好迹象。如果残差仍存在显著的自相关，则说明模型还有改进空间，可能需要调整(p,d,q)参数或考虑更复杂的模型。

整个项目做下来，我的体会是，ACF和PACF分析更像是一门“看图艺术”而非精确科学。它为我们提供了强有力的初始线索，但绝不能教条化。尤其是在面对复杂的现实数据时，一定要将图形解读与业务知识、统计检验（如Ljung-Box检验）以及模型选择准则（如AIC）结合起来。多试几组参数，多看看残差图，模型的效果才会更稳健。