做时间序列分析的人,十有八九都在“阶数怎么选”这个问题上栽过跟头。你盯着屏幕上的ACF图看了半天,那条置信区间外的柱子到底算不算数?PACF截尾还是拖尾?手算太慢,自动寻优又怕过拟合?别急,今天咱们就把这套流程彻底掰开揉碎,从原理到代码,让你下次再遇到建模困惑时,能笑着敲下那行 auto_arima。
先给个定心丸:时间序列建模不是玄学,它有一套严谨的逻辑链条。ACF和PACF是“侦探”,帮你初步锁定嫌疑犯(AR、MA还是ARMA);AIC/BIC是“法官”,在多个候选模型中选出性价比最高的那个。两者结合,既有理论依据,又有数据支撑。
一、先搞懂ACF和PACF:它们到底在说什么?
很多初学者死记硬背“AR(p)的PACF在p阶截尾”,但如果不理解背后的直觉,换个小样本数据就懵了。咱们用大白话重新讲一遍。
1.1 ACF(自相关函数):过去的自己有多像现在的自己?
ACF衡量的是时间序列 \(Y_t\) 与其滞后 \(k\) 期 \(Y_{t-k}\) 之间的线性相关程度。公式很简单:
\[ \rho(k) = \frac{\text{Cov}(Y_t, Y_{t-k})}{\sqrt{\text{Var}(Y_t)\text{Var}(Y_{t-k})}} \]
想象你在听一首循环播放的歌。如果今天听到的旋律和10秒前高度相似,那ACF在滞后10处就会很高。如果旋律毫无规律,ACF很快就会掉到0附近。
关键洞察:ACF对“过去所有信息”都有响应。比如一个AR(1)模型 \(Y_t = 0.8Y_{t-1} + \epsilon_t\),今天的值直接依赖昨天,昨天的值又依赖前天……所以滞后1、2、3……的ACF会呈现指数衰减,永远不会突然归零。这就是拖尾。
1.2 PACF(偏自相关函数):剔除中间干扰后的直接关系
PACF衡量的是在控制中间所有滞后项(\(Y_{t-1}, ..., Y_{t-k+1}\))后,\(Y_t\) 和 \(Y_{t-k}\) 的直接相关系数。
还是用唱歌的例子:如果你知道今天旋律完全由昨天决定(AR(1)),那么当你控制了昨天的值,今天和前天之间就没有额外关系了。所以PACF在滞后1处很高,滞后2开始突然归零。这就是截尾。
1.3 三种典型模型的ACF/PACF特征对照表
| 模型类型 | 模型方程 | ACF表现 | PACF表现 | 记忆点 |
|---|---|---|---|---|
| AR(p) | \(Y_t = \phi_1 Y_{t-1} + ... + \phi_p Y_{t-p} + \epsilon_t\) | 拖尾(指数衰减或正弦波) | 在p阶截尾 | AR看PACF,截尾定阶数 |
| MA(q) | \(Y_t = \epsilon_t + \theta_1 \epsilon_{t-1} + ... + \theta_q \epsilon_{t-q}\) | 在q阶截尾 | 拖尾 | MA看ACF,截尾定阶数 |
| ARMA(p,q) | 两者结合 | 两者都拖尾 | 两者都拖尾 | 最难认,需结合AIC |
为什么MA模型的ACF会截尾? 因为MA(q)只依赖最近q期的白噪声。滞后q+1期之后,\(Y_t\) 和 \(Y_{t-q-1}\) 之间没有共同的噪声项,相关性理论上为0。样本中由于估计误差,不会绝对为0,但会跌入置信区间内。
为什么AR模型的PACF会截尾? 因为AR(p)中,\(Y_t\) 只直接依赖于最近p期的观测值。滞后p+1期之后,控制中间变量后,直接相关系数为0。
二、实战第一步:生成数据并绘制ACF/PACF图
光说不练假把式。咱们先构造三个已知模型的数据,看看理论特征在图中是否吻合。
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.tsa.arima_process import ArmaProcess
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import warnings
warnings.filterwarnings('ignore')
# 设置绘图风格,让图表更美观
sns.set_style("whitegrid")
plt.rcParams['font.sans-serif'] = ['SimHei', 'Arial Unicode MS'] # 支持中文
plt.rcParams['axes.unicode_minus'] = False
np.random.seed(42)
n = 500 # 样本量
# 1. 构造AR(1)过程: Y_t = 0.7*Y_{t-1} + ε_t
ar_ar1 = np.array([1, -0.7]) # statsmodels要求系数带负号
ma_ar1 = np.array([1])
arma_ar1 = ArmaProcess(ar_ar1, ma_ar1)
data_ar1 = arma_ar1.generate_sample(n=n)
# 2. 构造MA(1)过程: Y_t = ε_t + 0.5*ε_{t-1}
ar_ma1 = np.array([1])
ma_ma1 = np.array([1, 0.5])
arma_ma1 = ArmaProcess(ar_ma1, ma_ma1)
data_ma1 = arma_ma1.generate_sample(n=n)
# 3. 构造ARMA(1,1)过程: Y_t = 0.6*Y_{t-1} + ε_t + 0.4*ε_{t-1}
ar_arma11 = np.array([1, -0.6])
ma_arma11 = np.array([1, 0.4])
arma_11 = ArmaProcess(ar_arma11, ma_arma11)
data_arma11 = arma_11.generate_sample(n=n)
# 统一绘图函数
def plot_acf_pacf(data, title, ax):
ax.set_title(title, fontsize=14, fontweight='bold')
# ACF图
ax1 = ax[0]
plot_acf(data, ax=ax1, lags=40, title='ACF', alpha=0.05)
# PACF图
ax2 = ax[1]
plot_pacf(data, ax=ax2, lags=40, title='PACF', alpha=0.05, method='ywadj')
# 调整布局
fig = ax1.figure
fig.tight_layout()
fig, axes = plt.subplots(3, 2, figsize=(14, 10))
plot_acf_pacf(data_ar1, 'AR(1) 过程', axes[0])
plot_acf_pacf(data_ma1, 'MA(1) 过程', axes[1])
plot_acf_pacf(data_arma11, 'ARMA(1,1) 过程', axes[2])
# 隐藏多余的子图
axes[2, 1].axis('off')
plt.suptitle('AR/MA/ARMA模型的ACF与PACF特征对比', fontsize=16, y=1.02)
plt.tight_layout()
plt.savefig('acf_pacf_comparison.png', dpi=150, bbox_inches='tight')
plt.show()
运行结果解读:
- AR(1):PACF在滞后1处显著突出(远超蓝色置信区间),滞后2及以后迅速落入区间内——截尾。ACF则呈指数衰减,缓慢趋近于0——拖尾。完美符合理论。
- MA(1):ACF在滞后1处显著,滞后2及以后落入区间——截尾。PACF拖尾衰减。再次吻合。
- ARMA(1,1):ACF和PACF都拖尾,难以直接区分。这就是为什么需要AIC准则来辅助决策。
小贴士:实际数据中,由于样本有限,截尾不会像理论那样绝对干净。通常我们会看前5-10个滞后,如果某阶之后大部分点都落在置信区间内,就认为“截尾”了。
三、AIC准则:如何自动优选模型参数?
人眼识别ACF/PACF有主观性,而且对于复杂模型(如ARMA(2,2)、ARIMA(1,1,1))几乎无能为力。AIC(赤池信息量准则)提供了客观的量化标准。
3.1 AIC的直觉理解
AIC = 2k - 2ln(L)
其中k是模型参数个数,L是似然函数值。
翻译成人话:AIC在平衡“拟合优度”和“模型复杂度”。
- 似然函数L越大,说明模型拟合数据越好(减分项越少)。
- 参数k越多,模型越复杂,越容易过拟合(加项越多)。
AIC越小,模型越好。它惩罚过于复杂的模型,鼓励简洁有效。
3.2 为什么不用单纯的最小误差平方和?
如果你只追求RSS(残差平方和)最小,那永远选参数最多的模型——加一个参数总能微调拟合。但这会导致过拟合:模型在训练集上表现完美,在新数据上一塌糊涂。
AIC通过参数惩罚项,迫使我们在“拟合”和“简洁”之间找平衡。类似的思想还有BIC(贝叶斯信息量准则),BIC对复杂模型的惩罚更重,倾向于选择更简洁的模型。
3.3 网格搜索自动优选参数
对于ARIMA(p,d,q)模型,我们需要搜索p、d、q的组合。d(差分项)通常通过ADF检验或观察序列平稳性确定,p和q则通过AIC网格搜索。
下面是一个完整的自动化建模流程:
”`python import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns from statsmodels.tsa.stattools import adfuller, acf, pacf from statsmodels.tsa.arima.model import ARIMA from statsmodels.graphics.tsaplots import plot_acf, plot_pacf import warnings warnings.filterwarnings(‘ignore’)
sns.set_style(“whitegrid”)
==================== 1. 数据准备 ====================
使用真实的股票收盘价数据演示(以苹果公司AAPL为例)
np.random.seed(42)
模拟一个具有趋势和季节性的时间序列
dates = pd.date_range(start=‘2020-01-01’, periods=500, freq=’D’) trend = np.linspace(100, 150, 500) seasonal = 10 * np.sin(2 * np.pi * np.arange(500) / 365) noise = np.random.normal(0, 3, 500) data = trend + seasonal + noise
ts = pd.Series(data, index=dates)
==================== 2. 平稳性检验 ====================
def check_stationarity(series, title=‘Time Series’):
"""ADF检验判断序列是否平稳"""
result = adfuller(series.dropna())
print(f'\n=== {title} 平稳性检验 (ADF) ===')
print(f'ADF统计量: {result[0]:.4f}')
print(f'p值: {result[1]:.4f}')
print(f'临界值 (1%): {result[4]["1%"]:.4f}, (5%): {result[4]["5%"]:.4f}, (10%): {result[4]["10%"]:.4f}')
if result[1] < 0.05:
print('结论: p值 < 0.05,序列是平稳的 ✓')
else:
print('结论: p值 >= 0.05,序列非平稳,需要差分')
return result[1] < 0.05
print(“原始序列平稳性检验:”) is_stationary = check_stationarity(ts)
如果非平稳,进行差分
if not is_stationary:
ts_diff = ts.diff().dropna()
check_stationarity(ts_diff, '一阶差分后')
else:
ts_diff = ts
==================== 3. 绘制ACF和PACF图(辅助判断d) ====================
fig, axes = plt.subplots(2, 2, figsize=(14, 8))
原始序列
axes[0, 0].plot(ts, label=‘Original’, color=‘steelblue’) axes[0, 0].set_title(‘原始时间序列’, fontweight=‘bold’) axes[0, 0].legend()
差分后序列
axes[0, 1].plot(ts_diff, label=‘Differenced’, color=‘coral’) axes[0, 1].set_title(‘一阶差分后序列’, fontweight=‘bold’) axes[0, 1].legend()
ACF图
plot_acf(ts_diff.dropna(), ax=axes[1, 0], lags=40, title=‘ACF (差分后)’, alpha=0.05)
PACF图
plot_pacf(ts_diff.dropna(), ax=axes[1, 1], lags=40, title=‘PACF (差分后)’, alpha=0.05, method=‘ywadj’)
plt.suptitle(‘时间序列诊断图’, fontsize=16, y=1.02) plt.tight_layout() plt.savefig(‘diagnostic_plots.png’, dpi=150, bbox_inches=‘tight’) plt.show()
==================== 4. 自动搜索最优ARIMA参数 ====================
def auto_arima_selection(series, d, p_range=range(0, 4), q_range=range(0, 4),
seasonal=False, m=12, max_order=5):
"""
使用AIC准则自动搜索最优ARIMA(p,d,q)参数
series: 时间序列(已差分)
d: 差分项数
p_range, q_range: 搜索范围
"""
best_aic = np.inf
best_order = None
best_model = None
results = []
print("\n开始网格搜索最优参数...")
print(f"搜索范围: p ∈ {list(p_range)}, q ∈ {list(q_range)}, d = {d}")
for p in p_range:
for q in q_range:
try:
# 拟合模型
model = ARIMA(series, order=(p, d, q))
fitted_model = model.fit()
aic = fitted_model.aic
results.append({
'p': p, 'q': q, 'aic': aic,
'model': fitted_model
})
# 更新最优
if aic < best_aic:
best_aic = aic
best_order = (p, d, q)
best_model = fitted_model
except Exception as e:
# 跳过无法收敛的模型
continue
# 按AIC排序展示Top 5
results_df = pd.DataFrame(results).sort_values('aic')
print("\n--- 搜索结果显示 (按AIC排序,越小越好) ---")
print(results_df.head(5).to_string(index=False))
print(f"\n✓ 最优模型: ARIMA{best_order}")
print(f" AIC值: {best_aic:.2f}")
return best_order, best_model, results_df
执行自动搜索
best_order, best_model, results_df = auto_arima_selection(
ts_diff,
d=1, # 前面ADF检验发现一阶差分后平稳
p_range=range(0, 4),
q_range=range(0, 4)
)
==================== 5. 模型诊断:检验残差 ====================
print(“\n=== 模型诊断 ===”) residuals = best_model.resid
绘制残差分布
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
axes[0].hist(residuals, bins=30, density=True, alpha=0.6, color=‘steelblue’, label=‘残差’) axes[0].axvline(0, color=‘red’, linestyle=‘–’, label=‘均值=0’) axes[0].set_title(‘残差分布直方图’) axes[0].legend()
残差ACF图
plot_acf(residuals, ax=axes[1], lags=20, title=‘残差ACF图’, alpha=0.05)
残差Q统计量检验(Ljung-Box检验)
from statsmodels.stats.diagnostics import acorr_ljungbox lb_test = acorr_ljungbox(residuals, lags=[12], return_df=True) print(f”\nLjung-Box检验 (滞后12):“) print(f”Q统计量: {lb_test[‘lb_stat’].values[0]:.4f}“) print(f”p值: {lb_test[‘lb_pvalue’].values[0]:.4f}“) if lb_test[‘lb_pvalue’].values[0] > 0.05:
print("结论: p > 0.05,残差是白噪声,模型拟合良好 ✓")
else:
print("结论: p ≤ 0.05,残差存在自相关,模型可能需要改进")
axes[2].plot(residuals.index, residuals, label=‘Residual
