本文目录导读：
我来详细介绍如何使用Python进行ARIMA时间序列分析。
安装必要的库
pip install pandas numpy matplotlib statsmodels
基础ARIMA模型实现
1 基础数据准备和建模
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import adfuller
import warnings
warnings.filterwarnings('ignore')
# 创建示例时间序列数据
np.random.seed(42)
dates = pd.date_range('2020-01-01', periods=200, freq='D')
# 生成带有趋势和季节性的数据
trend = np.linspace(100, 150, 200)
seasonal = 10 * np.sin(2 * np.pi * np.arange(200) / 30)
noise = np.random.normal(0, 2, 200)
data = trend + seasonal + noise
# 创建DataFrame
df = pd.DataFrame({'value': data}, index=dates)
print("数据预览：")
print(df.head())
2 数据平稳性检验
def check_stationarity(series):
"""检查时间序列的平稳性"""
result = adfuller(series.dropna())
print('ADF检验结果:')
print(f'ADF统计量: {result[0]:.6f}')
print(f'p值: {result[1]:.6f}')
print(f'临界值:')
for key, value in result[4].items():
print(f'\t{key}: {value:.6f}')
if result[1] <= 0.05:
print(" 序列是平稳的 (拒绝原假设)")
return True
else:
print(" 序列是非平稳的 (无法拒绝原假设)")
return False
# 检查原始序列平稳性
print("=== 原始序列平稳性检验 ===")
is_stationary = check_stationarity(df['value'])
# 如果非平稳，进行差分
if not is_stationary:
df['value_diff'] = df['value'].diff()
print("\n=== 一阶差分后序列平稳性检验 ===")
check_stationarity(df['value_diff'].dropna())
3 确定ARIMA参数(p,d,q)
def plot_acf_pacf(series, lags=40):
"""绘制ACF和PACF图"""
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
plot_acf(series.dropna(), lags=lags, ax=ax1)
plot_pacf(series.dropna(), lags=lags, ax=ax2)
plt.tight_layout()
plt.show()
# 绘制差分后的ACF和PACF图
print("绘制ACF和PACF图...")
plot_acf_pacf(df['value_diff'].dropna())
# 根据ACF和PACF图选择参数
# ACF: 确定MA(q) - 看截尾
# PACF: 确定AR(p) - 看截尾
# 一般先尝试简单参数: p=1, d=1, q=1
4 自动选择最优参数
from statsmodels.tsa.arima.model import ARIMA
from itertools import product
def find_best_arima_params(series, p_range=range(0, 4), d_range=range(0, 2), q_range=range(0, 4)):
"""自动搜索最优ARIMA参数"""
best_aic = float('inf')
best_order = None
results = []
for p, d, q in product(p_range, d_range, q_range):
try:
model = ARIMA(series.dropna(), order=(p, d, q))
result = model.fit()
aic = result.aic
results.append({
'order': (p, d, q),
'aic': aic
})
if aic < best_aic:
best_aic = aic
best_order = (p, d, q)
except:
continue
return best_order, best_aic, results
# 寻找最优参数
print("开始搜索最优ARIMA参数...")
best_order, best_aic, all_results = find_best_arima_params(df['value'])
print(f"最优参数: p={best_order[0]}, d={best_order[1]}, q={best_order[2]}")
print(f"最优AIC: {best_aic:.4f}")
# 显示前5个最优参数组合
all_results_sorted = sorted(all_results, key=lambda x: x['aic'])[:5]
print("\n前5个最优参数组合：")
for result in all_results_sorted:
print(f" ARIMA{result['order']}: AIC={result['aic']:.4f}")
5 训练ARIMA模型
# 使用最优参数训练模型
model = ARIMA(df['value'], order=best_order)
model_fit = model.fit()
print("\n=== 模型摘要 ===")
print(model_fit.summary())
6 模型诊断
from statsmodels.graphics.tsaplots import plot_acf
import seaborn as sns
def diagnose_model(model_fit):
"""模型诊断"""
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
# 残差图
residuals = model_fit.resid
axes[0, 0].plot(residuals)
axes[0, 0].set_title('残差图')
axes[0, 0].axhline(y=0, color='r', linestyle='--')
# 残差直方图
axes[0, 1].hist(residuals, bins=30, alpha=0.7, edgecolor='black')
axes[0, 1].set_title('残差分布')
# QQ图
from scipy import stats
stats.probplot(residuals, dist="norm", plot=axes[1, 0])
axes[1, 0].set_title('Q-Q图')
# 残差ACF
plot_acf(residuals.dropna(), lags=40, ax=axes[1, 1])
axes[1, 1].set_title('残差ACF')
plt.tight_layout()
plt.show()
# 正态性检验
stat, p_value = stats.normaltest(residuals.dropna())
print(f"正态性检验: 统计量={stat:.4f}, p值={p_value:.4f}")
# Durbin-Watson检验
from statsmodels.stats.stattools import durbin_watson
dw = durbin_watson(residuals.dropna())
print(f"Durbin-Watson统计量: {dw:.4f} (接近2表示无自相关)")
diagnose_model(model_fit)
7 预测与可视化
def forecast_and_plot(model_fit, df, steps=30):
"""预测并可视化"""
# 进行预测
forecast = model_fit.forecast(steps=steps)
forecast_index = pd.date_range(df.index[-1], periods=steps+1, freq='D')[1:]
# 获取置信区间
forecast_result = model_fit.get_forecast(steps=steps)
conf_int = forecast_result.conf_int()
# 可视化
plt.figure(figsize=(15, 6))
# 历史数据
plt.plot(df.index, df['value'], label='历史数据', color='blue')
# 预测数据
plt.plot(forecast_index, forecast, label='预测值', color='red', linewidth=2)
# 置信区间
plt.fill_between(
forecast_index,
conf_int.iloc[:, 0],
conf_int.iloc[:, 1],
color='red',
alpha=0.2,
label='95% 置信区间'
)
plt.title(f'ARIMA{best_order} 时间序列预测', fontsize=14)
plt.xlabel('日期')
plt.ylabel('值')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
return forecast, conf_int
# 进行预测
forecast_values, confidence_intervals = forecast_and_plot(model_fit, df, steps=30)
# 打印预测结果
print("\n=== 预测结果 ===")
forecast_df = pd.DataFrame({
'日期': forecast_index,
'预测值': forecast_values,
'下限(95%)': confidence_intervals.iloc[:, 0],
'上限(95%)': confidence_intervals.iloc[:, 1]
})
print(forecast_df.head(10))
实际应用示例：销售预测
# 生成模拟销售数据
np.random.seed(123)
dates_sales = pd.date_range('2023-01-01', periods=365, freq='D')
# 添加趋势、季节性和随机波动
base = 1000
weekly_pattern = 100 * np.sin(2 * np.pi * np.arange(365) / 7)
monthly_pattern = 200 * np.sin(2 * np.pi * np.arange(365) / 30)
growth = np.linspace(0, 300, 365)
noise = np.random.normal(0, 50, 365)
sales = base + weekly_pattern + monthly_pattern + growth + noise
sales = np.maximum(sales, 0) # 确保非负
df_sales = pd.DataFrame({'sales': sales}, index=dates_sales)
print("=== 销售数据ARIMA分析 ===")
# 1. 数据可视化
plt.figure(figsize=(15, 5))
plt.plot(df_sales.index, df_sales['sales'], label='日销售额', alpha=0.7)'每日销售数据')
plt.xlabel('日期')
plt.ylabel('销售额')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
# 2. 确定差分阶数
d_value = 1 # 通常销售额数据需要一阶差分
# 3. 搜索最优p和q
best_order, _, _ = find_best_arima_params(
df_sales['sales'],
p_range=range(0, 5),
d_range=[1],
q_range=range(0, 5)
)
print(f"最优ARIMA参数: {best_order}")
# 4. 训练模型
model_sales = ARIMA(df_sales['sales'], order=best_order)
model_sales_fitted = model_sales.fit()
# 5. 预测未来30天
forecast_sales = model_sales_fitted.forecast(steps=30)
forecast_index = pd.date_range(df_sales.index[-1], periods=31, freq='D')[1:]
# 6. 可视化结果
plt.figure(figsize=(15, 6))
plt.plot(df_sales.index[-90:], df_sales['sales'][-90:], label='历史数据(最近90天)', color='blue')
plt.plot(forecast_index, forecast_sales, label='30天预测', color='red', linewidth=2)f'销售额预测 (ARIMA{best_order})', fontsize=14)
plt.xlabel('日期')
plt.ylabel('销售额')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
# 7. 预测性能评估
from sklearn.metrics import mean_absolute_error, mean_squared_error
# 使用最后30天作为测试集
test_size = 30
train_data = df_sales['sales'][:-test_size]
test_data = df_sales['sales'][-test_size:]
model_temp = ARIMA(train_data, order=best_order)
model_temp_fit = model_temp.fit()
predictions = model_temp_fit.forecast(steps=test_size)
mae = mean_absolute_error(test_data, predictions)
rmse = np.sqrt(mean_squared_error(test_data, predictions))
print(f"\n=== 模型评估 ===")
print(f"MAE (平均绝对误差): {mae:.2f}")
print(f"RMSE (均方根误差): {rmse:.2f}")
print(f"MAPE (平均绝对百分比误差): {(mae/test_data.mean()*100):.2f}%")
高级技巧和注意事项
1 季节性ARIMA (SARIMA)
from statsmodels.tsa.statespace.sarimax import SARIMAX
# SARIMA模型(带季节性)
seasonal_order = (1, 1, 1, 7) # P, D, Q, S (S=7表示周季节性)
sarima_model = SARIMAX(
df_sales['sales'],
order=(best_order[0], best_order[1], best_order[2]),
seasonal_order=seasonal_order,
enforce_stationarity=False,
enforce_invertibility=False
)
sarima_result = sarima_model.fit()
print("SARIMA模型摘要：")
print(sarima_result.summary())
2 异常值处理
def detect_and_handle_outliers(series, threshold=3):
"""检测并处理异常值"""
mean = series.mean()
std = series.std()
# 检测异常值
outliers = np.abs(series - mean) > threshold * std
print(f"检测到 {outliers.sum()} 个异常值")
# 用中位数替换异常值
series_clean = series.copy()
median_val = series.median()
series_clean[outliers] = median_val
return series_clean, outliers
# 处理异常值
clean_data, outliers = detect_and_handle_outliers(df_sales['sales'])
df_sales['sales_clean'] = clean_data
3 模型比较
def compare_models(series, orders_list):
"""比较多个ARIMA模型"""
results = []
for order in orders_list:
try:
model = ARIMA(series.dropna(), order=order)
result = model.fit()
results.append({
'order': order,
'aic': result.aic,
'bic': result.bic,
'hqic': result.hqic
})
except:
print(f"ARIMA{order} 模型无法拟合")
results_df = pd.DataFrame(results)
return results_df.sort_values('aic')
# 比较不同模型
orders_to_compare = [(1,1,1), (2,1,2), (1,1,2), (2,1,1), (3,1,2)]
comparison_results = compare_models(df_sales['sales'], orders_to_compare)
print("\n模型比较结果：")
print(comparison_results)
4 滚动预测
def rolling_forecast(series, order, window_size=30):
"""滚动预测"""
predictions = []
actual_values = []
for i in range(window_size, len(series) - window_size):
train = series[:i]
test = series[i:i+window_size]
# 训练模型并预测
model = ARIMA(train, order=order)
fitted = model.fit()
forecast = fitted.forecast(steps=window_size)
predictions.extend(forecast)
actual_values.extend(test)
return np.array(predictions), np.array(actual_values)
# 执行滚动预测
window = 30
preds, actuals = rolling_forecast(
df_sales['sales'][:200],
best_order,
window_size=window
)
# 计算性能指标
mae_rolling = mean_absolute_error(actuals, preds)
rmse_rolling = np.sqrt(mean_squared_error(actuals, preds))
print(f"\n滚动预测评估:")
print(f"样本数: {len(preds)}")
print(f"MAE: {mae_rolling:.2f}")
print(f"RMSE: {rmse_rolling:.2f}")
完整的使用流程总结
- 数据准备: 确保数据是时间序列格式(pandas DatetimeIndex)
- 平稳性检验: 使用ADF检验，必要时进行差分
- 参数确定:
- 使用ACF/PACF图
- 或自动搜索最优参数
- 模型训练: 使用ARIMA或SARIMAX类
- 模型诊断: 检查残差是否满足白噪声假设
- 预测: 使用forecast()方法进行预测
- 评估: 使用MAE、RMSE等指标评估模型性能
这样你就可以用Python进行完整的ARIMA时间序列分析了！
