Python脚本数据时间序列ARIMA如何用

wen 实用脚本 2026-07-18 253

本文目录导读：

安装必要的库
基础ARIMA模型实现
实际应用示例：销售预测
高级技巧和注意事项
完整的使用流程总结

我来详细介绍如何使用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}")