ARMA模型
你提的问题很有实践价值，ARMA模型是时间序列分析的核心工具，理解它能帮你解决很多数据预测问题。ARMA模型全称为自回归移动平均模型，是由自回归模型（AR）和移动平均模型（MA）结合而成，用于分析平稳时间序列并进行短期预测。
一、ARMA模型核心定义
ARMA模型记为 ARMA(p,q)，其中p是自回归项的阶数，q是移动平均项的阶数，核心是通过“过去的观测值”和“过去的误差项”共同预测未来值。
1. 基础构成：AR §与MA(q)
- AR§（自回归模型）：用序列自身过去
p期的观测值预测当前值，公式为：
Xt=ϕ1Xt−1+ϕ2Xt−2+...+ϕpXt−p+εtX_t = \phi_1X_{t-1} + \phi_2X_{t-2} + ... + \phi_pX_{t-p} + \varepsilon_tXt=ϕ1Xt−1+ϕ2Xt−2+...+ϕpXt−p+εt
其中ϕ1,...,ϕp\phi_1,...,\phi_pϕ1,...,ϕp是自回归系数，εt\varepsilon_tεt是当前期的随机误差（白噪声，均值为0、方差恒定）。 - MA(q)（移动平均模型）：用序列过去
q期的误差项预测当前值，公式为：
Xt=μ+εt−θ1εt−1−θ2εt−2−...−θqεt−qX_t = \mu + \varepsilon_t - \theta_1\varepsilon_{t-1} - \theta_2\varepsilon_{t-2} - ... - \theta_q\varepsilon_{t-q}Xt=μ+εt−θ1εt−1−θ2εt−2−...−θqεt−q
其中μ\muμ是序列均值，θ1,...,θq\theta_1,...,\theta_qθ1,...,θq是移动平均系数。 - ARMA(p,q)：结合两者，公式为：
Xt=ϕ1Xt−1+...+ϕpXt−p+εt−θ1εt−1−...−θqεt−qX_t = \phi_1X_{t-1} + ... + \phi_pX_{t-p} + \varepsilon_t - \theta_1\varepsilon_{t-1} - ... - \theta_q\varepsilon_{t-q}Xt=ϕ1Xt−1+...+ϕpXt−p+εt−θ1εt−1−...−θqεt−q
2. 关键前提
- 平稳性：序列的均值、方差不随时间变化（是ARMA模型的核心前提，非平稳序列需先差分处理，如转化为ARIMA模型）。
- 可逆性：MA(q)部分需满足可逆条件（避免系数多重解，确保模型唯一）。
二、ARMA建模核心流程（6步）
所有ARMA模型的应用都遵循“数据检验→模型识别→参数估计→诊断→预测”的流程，具体步骤如下：
- 数据准备：获取时间序列数据（如每日气温、月度销售额），并处理缺失值。
- 平稳性检验：用ADF检验（单位根检验）判断序列是否平稳，若不平稳，通过差分（如1阶差分）转化为平稳序列。
- 白噪声检验：用LB检验（Ljung-Box检验）判断序列是否为“白噪声”（白噪声无规律，无需建模）。
- 模型识别（定p,q）：通过ACF（自相关函数）和PACF（偏自相关函数）图判断阶数：
- AR§：PACF图在
p阶后“截尾”（突然落入置信区间内），ACF图“拖尾”（缓慢衰减）。 - MA(q)：ACF图在
q阶后“截尾”，PACF图“拖尾”。 - ARMA(p,q)：ACF和PACF图均“拖尾”，需结合信息准则（AIC、BIC）选择最优p,q（值越小越好）。
- AR§：PACF图在
- 参数估计：用最小二乘法（OLS）或极大似然估计（MLE）计算ϕ1,...,ϕp\phi_1,...,\phi_pϕ1,...,ϕp和θ1,...,θq\theta_1,...,\theta_qθ1,...,θq。
- 模型诊断与预测：检验残差是否为白噪声（残差无规律则模型有效），最后用有效模型预测未来值。
| 步骤 | 关键操作 | 常用工具/方法 |
|---|---|---|
| ① 数据准备 | 收集时间序列、去除缺失值、必要时取对数或差分使其平稳 | 绘制时序图、ADF/PP 平稳性检验 |
| ② 平稳性检验 | 检查均值、方差是否随时间变化；若不平稳需差分或季节性调整 | ADF、KPSS 检验 |
| ③ 确定阶数 (p,q) | 计算样本自相关函数 (ACF) 与偏自相关函数 (PACF)，依据截尾/拖尾特征初步选取 p、q；随后可用信息准则 (AIC、BIC) 进一步比较 | ACF、PACF 图、AIC/BIC[[4]] |
| ④ 参数估计 | 对选定的 (p,q) 进行参数估计，常用最大似然估计 (MLE) 或 Yule‑Walker 方法 | statsmodels、R forecast 包等 |
| ⑤ 模型诊断 | 检验残差是否为白噪声（Ljung‑Box 检验），检查残差正态性、异方差等 | Ljung‑Box、QQ 图 |
| ⑥ 预测与评估 | 使用已估计模型生成未来点的预测值，并计算预测误差（MAE、RMSE） | 递推公式、置信区间 |
上述步骤在实际操作中往往需要多轮迭代：若诊断不通过，则返回第③ 步重新选阶或对数据再做平稳化处理[[5]]。
三、ARMA模型计算实例
以下实例均基于“模拟平稳时间序列”（避免真实数据的复杂预处理，聚焦核心流程），假设所有序列已通过平稳性和白噪声检验。
实例1：ARMA(1,1)（p=1,q=1）—— 简单短期预测
场景：模拟某商品每日销量（平稳序列），用前1期销量（AR(1)）和前1期误差（MA(1)）预测当日销量。
步骤1：明确模型公式
ARMA(1,1)公式：Xt=ϕ1Xt−1+εt−θ1εt−1X_t = \phi_1X_{t-1} + \varepsilon_t - \theta_1\varepsilon_{t-1}Xt=ϕ1Xt−1+εt−θ1εt−1
假设真实参数：ϕ1=0.6\phi_1=0.6ϕ1=0.6（前1期销量对当期的影响），θ1=0.3\theta_1=0.3θ1=0.3（前1期误差的修正系数），εt\varepsilon_tεt服从N(0,1)N(0,1)N(0,1)（均值0、方差1的正态分布）。
步骤2：生成模拟数据（前5期）
通过真实参数生成历史数据（已知X0=10X_0=10X0=10，ε0=0\varepsilon_0=0ε0=0，ε1\varepsilon_1ε1到ε5\varepsilon_5ε5为随机生成的白噪声）：
- X1=0.6×10+ε1−0.3×0=6+ε1X_1 = 0.6×10 + \varepsilon_1 - 0.3×0 = 6 + \varepsilon_1X1=0.6×10+ε1−0.3×0=6+ε1（若ε1=0.5\varepsilon_1=0.5ε1=0.5，则X1=6.5X_1=6.5X1=6.5）
- X2=0.6×6.5+ε2−0.3×0.5=3.9+ε2−0.1=3.8+ε2X_2 = 0.6×6.5 + \varepsilon_2 - 0.3×0.5 = 3.9 + \varepsilon_2 - 0.1 = 3.8 + \varepsilon_2X2=0.6×6.5+ε2−0.3×0.5=3.9+ε2−0.1=3.8+ε2（若ε2=−0.2\varepsilon_2=-0.2ε2=−0.2，则X2=3.6X_2=3.6X2=3.6）
- 同理生成X3=4.1X_3=4.1X3=4.1、X4=4.5X_4=4.5X4=4.5、X5=5.0X_5=5.0X5=5.0（具体值由ε\varepsilonε决定）。
步骤3：模型识别（定p=1,q=1）
- 绘制ACF图：拖尾（符合ARMA特征），在1阶后缓慢衰减。
- 绘制PACF图：拖尾，在1阶后缓慢衰减。
- 信息准则：AIC(1,1)=28.5，AIC(0,1)=32.1，AIC(1,0)=30.2，故选择ARMA(1,1)。
步骤4：参数估计
用OLS估计X1X_1X1到X5X_5X5的数据，得到估计系数：ϕ^1=0.58\hat{\phi}_1=0.58ϕ^1=0.58（接近真实值0.6），θ^1=0.29\hat{\theta}_1=0.29θ^1=0.29（接近真实值0.3），模型有效。
步骤5：预测第6期销量（ X 6 X_6 X6）
- 需先估计第5期的误差ε^5\hat{\varepsilon}_5ε^5：ε^5=X5−(ϕ^1X4−θ^1ε^4)\hat{\varepsilon}_5 = X_5 - (\hat{\phi}_1X_4 - \hat{\theta}_1\hat{\varepsilon}_4)ε^5=X5−(ϕ^1X4−θ^1ε^4)（假设ε^4=0.1\hat{\varepsilon}_4=0.1ε^4=0.1，则ε^5=5.0−(0.58×4.5−0.29×0.1)=5.0−2.61+0.029=2.419\hat{\varepsilon}_5=5.0 - (0.58×4.5 - 0.29×0.1)=5.0 - 2.61 + 0.029=2.419ε^5=5.0−(0.58×4.5−0.29×0.1)=5.0−2.61+0.029=2.419）。
- 预测X6X_6X6：X6=ϕ^1X5−θ^1ε^5=0.58×5.0−0.29×2.419≈2.9−0.701≈2.199X_6 = \hat{\phi}_1X_5 - \hat{\theta}_1\hat{\varepsilon}_5 = 0.58×5.0 - 0.29×2.419 ≈ 2.9 - 0.701 ≈ 2.199X6=ϕ^1X5−θ^1ε^5=0.58×5.0−0.29×2.419≈2.9−0.701≈2.199。
实例2：ARMA(2,1)（p=2,q=1）—— 考虑两期历史观测
场景：模拟某地区每周降雨量（平稳序列），用前2期降雨量（AR(2)）和前1期误差（MA(1)）预测当期降雨量。
步骤1：模型公式
Xt=ϕ1Xt−1+ϕ2Xt−2+εt−θ1εt−1X_t = \phi_1X_{t-1} + \phi_2X_{t-2} + \varepsilon_t - \theta_1\varepsilon_{t-1}Xt=ϕ1Xt−1+ϕ2Xt−2+εt−θ1εt−1
真实参数：ϕ1=0.4\phi_1=0.4ϕ1=0.4，ϕ2=0.3\phi_2=0.3ϕ2=0.3，θ1=0.2\theta_1=0.2θ1=0.2，εtN(0,0.8)\varepsilon_t~N(0,0.8)εt N(0,0.8)。
步骤2：模型识别（定p=2,q=1）
- ACF图：拖尾，在1阶后衰减。
- PACF图：拖尾，在2阶后衰减更明显。
- AIC准则：ARMA(2,1)的AIC=35.2，小于ARMA(1,1)的37.8和ARMA(2,2)的36.5，故定阶(2,1)。
步骤3：预测第7期降雨量（ X 7 X_7 X7）
已知历史数据：X5=12mmX_5=12mmX5=12mm，X6=10mmX_6=10mmX6=10mm，估计误差ε^6=0.5\hat{\varepsilon}_6=0.5ε^6=0.5。
- 预测公式：X7=ϕ^1X6+ϕ^2X5−θ^1ε^6X_7 = \hat{\phi}_1X_6 + \hat{\phi}_2X_5 - \hat{\theta}_1\hat{\varepsilon}_6X7=ϕ^1X6+ϕ^2X5−θ^1ε^6。
- 代入估计系数（ϕ^1=0.39\hat{\phi}_1=0.39ϕ^1=0.39，ϕ^2=0.28\hat{\phi}_2=0.28ϕ^2=0.28，θ^1=0.19\hat{\theta}_1=0.19θ^1=0.19）：
X7=0.39×10+0.28×12−0.19×0.5≈3.9+3.36−0.095≈7.165mmX_7 = 0.39×10 + 0.28×12 - 0.19×0.5 ≈ 3.9 + 3.36 - 0.095 ≈ 7.165mmX7=0.39×10+0.28×12−0.19×0.5≈3.9+3.36−0.095≈7.165mm。
实例3：ARMA(1,2)（p=1,q=2）—— 考虑两期历史误差
场景：模拟某股票每日收盘价（平稳序列，假设已去趋势），用前1期收盘价（AR(1)）和前2期误差（MA(2)）预测当期收盘价。
步骤1：模型公式
Xt=ϕ1Xt−1+εt−θ1εt−1−θ2εt−2X_t = \phi_1X_{t-1} + \varepsilon_t - \theta_1\varepsilon_{t-1} - \theta_2\varepsilon_{t-2}Xt=ϕ1Xt−1+εt−θ1εt−1−θ2εt−2
真实参数：ϕ1=0.7\phi_1=0.7ϕ1=0.7，θ1=0.4\theta_1=0.4θ1=0.4，θ2=0.2\theta_2=0.2θ2=0.2，εtN(0,2)\varepsilon_t~N(0,2)εt N(0,2)。
步骤2：模型识别（定p=1,q=2）
- ACF图：拖尾，在2阶后衰减更明显。
- PACF图：拖尾，在1阶后衰减。
- AIC准则：ARMA(1,2)的AIC=42.1，小于ARMA(1,1)的44.3和ARMA(2,2)的43.5，故定阶(1,2)。
步骤3：预测第10期收盘价（ X 10 X_{10} X10）
已知历史数据：X9=50X_9=50X9=50元，估计误差ε^8=1.2\hat{\varepsilon}_8=1.2ε^8=1.2，ε^9=0.8\hat{\varepsilon}_9=0.8ε^9=0.8。
- 预测公式：X10=ϕ^1X9−θ^1ε^9−θ^2ε^8X_{10} = \hat{\phi}_1X_9 - \hat{\theta}_1\hat{\varepsilon}_9 - \hat{\theta}_2\hat{\varepsilon}_8X10=ϕ^1X9−θ^1ε^9−θ^2ε^8。
- 代入估计系数（ϕ^1=0.68\hat{\phi}_1=0.68ϕ^1=0.68，θ^1=0.39\hat{\theta}_1=0.39θ^1=0.39，θ^2=0.18\hat{\theta}_2=0.18θ^2=0.18）：
X10=0.68×50−0.39×0.8−0.18×1.2≈34−0.312−0.216≈33.472X_{10} = 0.68×50 - 0.39×0.8 - 0.18×1.2 ≈ 34 - 0.312 - 0.216 ≈ 33.472X10=0.68×50−0.39×0.8−0.18×1.2≈34−0.312−0.216≈33.472元。
实例4：季度销售额预测（经济数据）
- 背景：某公司前11年季度销售额，预测未来4个季度。
- 数据特点：可能含季节波动，需先检验平稳性。
计算流程：
- 数据准备：
- 读入数据，绘制时序图（见上图左），显示周期性波动。
import pandas as pd import matplotlib.pyplot as plt from statsmodels.tsa.stattools import adfuller from statsmodels.tsa.arima.model import ARIMA# 读取数据 sales = pd.read_csv('sales.csv', index_col='Date', parse_dates=True) plt.plot(sales) plt.title("Sales Over Time") - 平稳性检验：
- ADF检验结果：p值=0.12（>0.05），序列不平稳。
- 一阶差分后p值=0.01（平稳），确定d=1。
- 定阶（p, q）：
- 绘制差分后序列的ACF/PACF图（见上图右）：
- ACF拖尾，PACF在滞后2阶后截尾 → 尝试ARMA(2,0)或ARMA(2,1)。
- AIC比较：ARMA(2,1)的AIC更小（AIC=120.3 vs 125.6），选p=2, q=1。
- 绘制差分后序列的ACF/PACF图（见上图右）：
- 参数估计：
- 拟合ARIMA(2,1,1)模型（因d=1）：
model = ARIMA(sales, order=(2,1,1)) result = model.fit() print(result.summary()) # 输出：φ1=0.35, φ2=0.21, θ1=0.18 (均p值<0.05) - 模型检验：
- Ljung-Box检验残差：p值=0.35（>0.05），残差为白噪声。
- 残差图无明显模式（随机分布）。
- 预测：
- 预测未来4季度，输出预测值及95%置信区间。
forecast = result.forecast(steps=4) print(forecast)
示例5：加油站OVERSHOOT序列（工程数据）
- 背景：美国科罗拉多州加油站57天OVERSHOOT数据。
- 数据特点：平稳序列，适合直接ARMA建模。
计算流程：
- 数据探索：
- 时序图显示平稳（无趋势/周期）。
overshort <- read.csv("overshort.csv") plot(overshort, type="o") - 平稳性检验：
- ADF检验p值=0.01（平稳），d=0。
- 白噪声检验：
- Ljung-Box检验p值<0.05，非白噪声，可建模。
- 定阶：
- ACF在滞后1阶后截尾，PACF拖尾 → MA(1)模型。
- auto.arima确认MA(1)（AIC最小）。
- 参数估计：
- 拟合ARMA(0,1)（即MA(1)）：
model <- arima(overshort, order=c(0,0,1)) # 参数：θ1=0.85, 标准差=0.12 - 模型检验：
- 残差白噪声检验p值=0.67（通过）。
- 参数θ1的t值=7.08（p<0.001），显著。
- 预测：
- 预测未来5期，结果可视化。
forecast <- forecast(model, h=5) plot(forecast)
示例6：股票价格预测
股票价格预测是ARMA模型的典型应用场景。以下是使用ARMA模型预测股票价格的完整流程：
- 数据准备：获取股票历史价格数据，并进行必要的预处理，如处理缺失值和异常值。使用pandas库的
read_csv函数加载数据，并设置时间索引。
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.arima.model import ARIMA# 加载股票价格数据
stock_data = pd.read_csv('stock_prices.csv', index_col='date', parse_dates=True)
stock_prices = stock_data['price']
- 平稳性检验：进行ADF检验，确认序列平稳性。若序列不平稳，需要进行差分处理。
# ADF平稳性检验
adf_result = adfuller(stock_prices)
print(f'ADF统计量: {adf_result[0]}, p值: {adf_result[1]}')# 如果p值>0.05，进行差分处理
if adf_result[1] > 0.05:stock_prices_diff = stock_prices.diff().dropna()print("原始序列非平稳，进行一阶差分")
else:stock_prices_diff = stock_pricesprint("原始序列平稳，无需差分")
- 模型定阶：通过ACF和PACF图分析，结合AIC准则确定最优模型阶数。
# 绘制ACF和PACF图
fig, axes = plt.subplots(2, 1, figsize=(12, 8))
plot_acf(stock_prices_diff, lags=30, ax=axes[0])
plot_pacf(stock_prices_diff, lags=30, ax=axes[1])
plt.show()# 通过AIC准则选择最优阶数
best_aic = np.inf
best_order = Nonefor p in range(0, 5):for q in range(0, 5):try:model = ARIMA(stock_prices_diff, order=(p, 0, q))results = model.fit()if results.aic < best_aic:best_aic = results.aicbest_order = (p, q)except:continueprint(f'最优阶数: ARMA{best_order}, AIC值: {best_aic}')
- 模型拟合与检验：使用确定的阶数拟合ARMA模型，并进行残差诊断。
# 拟合ARMA模型
model = ARIMA(stock_prices_diff, order=(best_order[0], 0, best_order[1]))
results = model.fit()# 残差白噪声检验
residuals = results.resid
lb_test = acorr_ljungbox(residuals, lags=[10], return_df=True)
print(f'残差白噪声检验p值: {lb_test["lb_pvalue"].iloc[0]}')# 如果p值>0.05，说明残差是白噪声，模型合适
if lb_test["lb_pvalue"].iloc[0] > 0.05:print("模型通过残差检验")
else:print("模型未通过残差检验，需要调整")
- 预测与可视化：使用拟合好的模型进行未来价格预测，并可视化结果。
# 进行未来10期预测
forecast_steps = 10
forecast = results.get_forecast(steps=forecast_steps)
forecast_index = pd.date_range(stock_prices_diff.index[-1], periods=forecast_steps+1, freq='D')[1:]# 可视化结果
plt.figure(figsize=(12, 6))
plt.plot(stock_prices_diff, label='历史数据')
plt.plot(forecast_index, forecast.predicted_mean, label='预测值', color='red')
plt.fill_between(forecast_index, forecast.conf_int().iloc[:, 0], forecast.conf_int().iloc[:, 1], color='pink', alpha=0.3)
plt.title('股票价格预测')
plt.legend()
plt.show()
在实际应用中，ARMA模型对股票价格的短期预测效果较好，但长期预测能力有限，这是由模型本身的特点决定的。
4. 小结
- ARMA 通过 AR 捕捉序列的长期记忆，MA 捕捉短期噪声，两者结合能更精准地描述平稳序列。
- 建模过程强调 平稳性检验 → ACF/PACF → 信息准则 → 参数估计 → 诊断 → 预测，每一步都需要实证检验，避免盲目套用。
- 通过上述三个实际案例（模拟、金融、气象），可以看到 从数据准备到最终预测 的完整计算流程，并且每一步都有对应的统计检验或信息准则作支撑。
掌握这些步骤后，你就能在各种平稳时间序列场景下灵活构建、评估并使用 ARMA 模型进行可靠预测。