经典时间序列建模：差分处理与 AR/MA/ARMA 模型识别

2. 构造非平稳序列并观察差分效果

随机游走（Random Walk）是时间序列里非常经典的'非平稳'例子，它会随时间漂移，没有稳定的均值。叠加一个线性趋势后，非平稳性会更明显。

# 生成非平稳数据：随机游走 + 线性趋势
np.random.seed(42)
n = 500
random_walk = np.random.randn(n).cumsum()  # 累积和 -> 随机游走
trend = np.linspace(0, 100, n)  # 线性趋势
data = pd.Series(random_walk + trend)
data.index = pd.date_range(start="2022-01-01", periods=n, freq="D")

plt.figure(figsize=(12, 4))
plt.plot(data)
plt.title("原始序列：随机游走 + 趋势（明显非平稳）")
plt.tight_layout()
plt.show()
adf_test(data, name="原始序列")

可以看到 ADF 检验的 p-value 接近 1，说明无法拒绝原假设（非平稳）。接下来我们尝试一阶差分。

# 一阶差分：尝试消除趋势
data_diff_1 = data.diff().dropna()

plt.figure(figsize=(12, 4))
plt.plot(data_diff_1)
plt.axhline(0, color="gray", linewidth=1, linestyle="--")
plt.title("一阶差分后的序列（通常更接近平稳）")
plt.tight_layout()
plt.show()
adf_test(data_diff_1, name="一阶差分后")

此时 p-value 降为 0，统计量显著小于临界值，说明序列已经变得平稳了。

二、处理季节性：季节性差分

很多序列（销量、客流、气温）会在固定周期内重复波动，这就是季节性。季节性差分的思想与普通差分相同，只是'减去的不是上一期'，而是'减去上一个周期的同位置'。

公式：$\Delta_s y_t = y_t - y_{t-s}$

• 月度数据的年度季节：s = 12
• 季度数据的年度季节：s = 4

如果序列同时存在趋势 + 季节性，常见做法是先做季节性差分（去周期），再对结果做一阶差分（去趋势）。

# 生成一段'趋势 + 季节性 + 噪声'的月度序列，并做季节性差分
np.random.seed(7)
# 4 年月度数据（48 个点），周期为 12
time_index = pd.date_range(start="2020-01-01", periods=48, freq="M")
seasonal = 10 * np.sin(np.arange(48) * (2 * np.pi / 12))  # 振幅 10，周期 12
trend = np.linspace(0, 20, 48)
noise = np.random.randn(48) * 2
seasonal_data = pd.Series(seasonal + trend + noise, index=time_index)
seasonal_diff_12 = seasonal_data.diff(periods=12).dropna()

fig, axes = plt.subplots(2, 1, figsize=(12, 6), sharex=True)
axes[0].plot(seasonal_data)
axes[0].set_title("原始季节性序列（趋势 + 季节 + 噪声）")
axes[1].plot(seasonal_diff_12)
axes[1].axhline(0, color="gray", linewidth=1, linestyle="--")
axes[1].set_title("季节性差分后（s=12）")
plt.tight_layout()
plt.show()
adf_test(seasonal_data, name="原始季节性序列")
adf_test(seasonal_diff_12, name="季节性差分后（s=12）")

从输出结果看，季节性差分显著降低了 p-value，消除了周期性带来的非平稳影响。

三、模型选择：AR / MA / ARMA 与 ACF/PACF

现在假设我们已经拿到了平稳序列（原序列平稳，或通过差分变得平稳）。接下来就是定阶的关键步骤。

1. ACF 与 PACF：分别在看什么？

• ACF（自相关）：$y_t$ 与 $y_{t-k}$ 的相关性（包含间接影响）。
• PACF（偏自相关）：控制了中间滞后项后，$y_t$ 与 $y_{t-k}$ 的'直接相关'。

2. '截尾'与'拖尾'经验规律

注意：这是非常常用的入门规律，但真实业务数据里常需要结合 AIC/BIC、残差白噪声检验等一起判断。

3. AR(p)：自回归（'惯性/记忆'）

AR(p) 的直观解释：当前值由过去 p 期的值线性决定（再加上随机扰动）。

$$y_t = \phi_1 y_{t-1} + \cdots + \phi_p y_{t-p} + \varepsilon_t$$

下面模拟一个 AR(2) 序列，并用 ACF/PACF 验证'PACF 在 2 阶后截尾'的现象。

# --- 案例一：AR(2) ---
print("--- 案例一：AR(2) 模型 ---")
# AR(2)：y_t = 0.7 y_{t-1} + 0.2 y_{t-2} + ε_t
ar_params = np.array([0.7, 0.2])
ma_params = np.array([])
ar_process = ArmaProcess.from_coeffs(arcoefs=ar_params, macoefs=ma_params)

np.random.seed(100)
ar_data = ar_process.generate_sample(nsample=500)

fig, axes = plt.subplots(3, 1, figsize=(10, 9))
axes[0].plot(ar_data)
axes[0].set_title("模拟序列：AR(2)")
plot_acf(ar_data, ax=axes[1], lags=20, title="ACF（AR 通常拖尾）")
plot_pacf(ar_data, ax=axes[2], lags=20, title="PACF（AR(p) 通常在 p 阶截尾）")
plt.tight_layout()
plt.show()

# 使用 ARIMA(p,d,q) 形式拟合：AR(2) 等价于 ARIMA(2,0,0)
model_ar = ARIMA(ar_data, order=(2, 0, 0)).fit()
print(model_ar.summary())

从拟合结果可以看出，AR 系数显著不为 0，且残差诊断基本通过，说明模型有效。

4. MA(q)：移动平均（'短期冲击'）

MA(q) 的直观解释：当前值由当前与过去 q 期的随机冲击（误差项）叠加形成。

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

它很适合描述'突发冲击会影响一小段时间然后消失'的现象，例如生产线短暂故障、突发促销等。MA(q) 的典型经验规律是：ACF 在 q 阶后截尾，而 PACF 往往拖尾。

# --- 案例二：MA(2) ---
print("--- 案例二：MA(2) 模型 ---")
# MA(2)：y_t = ε_t + 0.8 ε_{t-1} + 0.4 ε_{t-2}
ar_params = np.array([])
ma_params = np.array([0.8, 0.4])
ma_process = ArmaProcess.from_coeffs(arcoefs=ar_params, macoefs=ma_params)

np.random.seed(200)
ma_data = ma_process.generate_sample(nsample=500)

fig, axes = plt.subplots(3, 1, figsize=(10, 9))
axes[0].plot(ma_data)
axes[0].set_title("模拟序列：MA(2)")
plot_acf(ma_data, ax=axes[1], lags=20, title="ACF（MA(q) 通常在 q 阶截尾）")
plot_pacf(ma_data, ax=axes[2], lags=20, title="PACF（MA 通常拖尾）")
plt.tight_layout()
plt.show()

# MA(2) 等价于 ARIMA(0,0,2)
model_ma = ARIMA(ma_data, order=(0, 0, 2)).fit()
print(model_ma.summary())

5. ARMA(p,q)：惯性 + 冲击的混合

真实世界经常同时存在惯性（AR）和短期冲击（MA）。ARMA(p,q) 的经验特征是：ACF 与 PACF 往往都会拖尾，因此仅凭图形精确定阶会更难。常见策略是从低阶（如 (1,1)）开始，结合 AIC/BIC 与残差诊断逐步选择。

# --- 案例三：ARMA(1,1) ---
print("--- 案例三：ARMA(1,1) 模型 ---")
# ARMA(1,1)：y_t = 0.8 y_{t-1} + ε_t + 0.5 ε_{t-1}
ar_params = np.array([0.8])
ma_params = np.array([0.5])
arma_process = ArmaProcess.from_coeffs(arcoefs=ar_params, macoefs=ma_params)

np.random.seed(300)
arma_data = arma_process.generate_sample(nsample=500)

fig, axes = plt.subplots(3, 1, figsize=(10, 9))
axes[0].plot(arma_data)
axes[0].set_title("模拟序列：ARMA(1,1)")
plot_acf(arma_data, ax=axes[1], lags=20, title="ACF（ARMA 常拖尾）")
plot_pacf(arma_data, ax=axes[2], lags=20, title="PACF（ARMA 常拖尾）")
plt.tight_layout()
plt.show()

# ARMA(1,1) 等价于 ARIMA(1,0,1)
model_arma = ARIMA(arma_data, order=(1, 0, 1)).fit()
print(model_arma.summary())

总结

今天我们梳理了时间序列建模的基础流程：

1. 趋势处理：遇到非平稳序列，先用一阶或二阶差分（Series.diff()）消除趋势。
2. 季节性处理：对于有固定周期的数据，使用季节性差分（Series.diff(periods=s)）。
3. 模型定阶：在平稳序列上观察 ACF/PACF 的截尾与拖尾特征，初步判断 AR、MA 或 ARMA 的阶数。

最后补充一下 ARIMA(p, d, q) 的参数含义：

• p：AR 阶数，常从 PACF 的截尾特征获取线索。
• d：差分次数，让序列更接近平稳。
• q：MA 阶数，常从 ACF 的截尾特征获取线索。

实际应用中，不要迷信单一指标，结合 AIC/BIC 与残差白噪声检验综合判断会更稳妥。