Python day58
@浙大疏锦行
ARIMA建模步骤：
- 数据可视化：观察原始时间序列图，判断是否存在趋势或季节性。
- 平稳性检验：
- 对原始序列进行ADF检验。
- 如果p值 > 0.05，说明序列非平稳，需要进行差分。
- 确定差分次数 d:
- 进行一阶差分，然后再次进行ADF检验。
- 如果平稳了，则 d=1。否则，继续差分，直到平稳。
- 确定 p 和 q:
- 对差分后的平稳序列绘制ACF和PACF图。
- 根据昨天学习的规则（PACF定p，ACF定q）来选择p和q的值。
- 建立并训练ARIMA(p, d, q)模型。
- 模型评估与诊断：查看模型的摘要信息，检查残差是否为白噪声。
- 进行预测
时间序列建模的迭代过程：初步选择 -> 拟合检验 -> 模型修正 -> 数据还原。
三个时序数据集大概图示
import pandas as pd
url = 'https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv'
df_air = pd.read_csv(url, header=0, index_col=0, parse_dates=True)
print(df_air.head())
df_air.plot()
from statsmodels.datasets import sunspots
df_sun = sunspots.load_pandas().data['SUNACTIVITY']
print(df_sun.head())
df_sun.plot()url = 'https://raw.githubusercontent.com/jbrownlee/Datasets/master/daily-total-female-births.csv'
df_births = pd.read_csv(url, header=0, index_col=0, parse_dates=True)
# df_births.plot()
df_births.plot()实战
import pandas as pd
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# 设置matplotlib以正确显示中文
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False# 加载数据
url = 'https://raw.githubusercontent.com/jbrownlee/Datasets/master/daily-total-female-births.csv'
df = pd.read_csv(url, header=0, index_col=0, parse_dates=True)
df.columns = ['Births']
ts_data = df['Births']print("--- 原始数据预览 ---")
print(ts_data.head())# 绘制时序图
plt.figure(figsize=(14, 7))
plt.plot(ts_data)
plt.title('1959年加州每日女性出生数量')
plt.xlabel('日期')
plt.ylabel('出生数量')
plt.show()def adf_test(timeseries):print('--- ADF检验结果 ---')# H0: 序列非平稳; H1: 序列平稳result = adfuller(timeseries)print(f'ADF Statistic: {result[0]}')print(f'p-value: {result[1]}') # 重点关注这个值if result[1] <= 0.05:print("结论: 成功拒绝原假设，序列是平稳的。")else:print("结论: 未能拒绝原假设，序列是非平稳的。")adf_test(ts_data)fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
plot_acf(ts_data, ax=ax1, lags=40)
ax1.set_title('自相关函数 (ACF)')
plot_pacf(ts_data, ax=ax2, lags=40)
ax2.set_title('偏自相关函数 (PACF)')
plt.tight_layout()
plt.show()以此确定p、q候选的ARIMA模型
import warnings
warnings.filterwarnings("ignore")
# 建立ARIMA(2, 0, 0)模型
model = ARIMA(ts_data, order=(2, 0, 0))
arima_result = model.fit()# 打印模型摘要
print(arima_result.summary())# 让我们预测未来30天
forecast_steps = 30
forecast = arima_result.get_forecast(steps=forecast_steps)
pred_mean = forecast.predicted_mean
conf_int = forecast.conf_int()# 绘制结果
plt.figure(figsize=(14, 7))
plt.plot(ts_data, label='原始数据')
# 绘制模型在历史数据上的拟合值
plt.plot(arima_result.fittedvalues, color='orange', label='模型拟合值')
# 绘制未来预测值
plt.plot(pred_mean, color='red', label='未来预测值')
# 绘制置信区间
plt.fill_between(conf_int.index,conf_int.iloc[:, 0],conf_int.iloc[:, 1], color='pink', alpha=0.5, label='95%置信区间')
plt.title('ARIMA(2,0,0)模型拟合与预测')
plt.legend()
plt.show()航空乘客数量预测实战
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import warningswarnings.filterwarnings("ignore") # 忽略一些模型拟合时产生的警告信息
# 加载数据
url = 'https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv'
df_air = pd.read_csv(url, header=0, index_col=0, parse_dates=True)
df_air.rename(columns={'#Passengers': 'Passengers'}, inplace=True) # 列名简化# 绘制原始数据
plt.figure(figsize=(12, 6))
plt.plot(df_air['Passengers'])
plt.title('Airline Passengers (1949-1960)')
plt.xlabel('Year')
plt.ylabel('Number of Passengers')
plt.grid(True)
plt.show()# 对数变换
df_air['log_passengers'] = np.log(df_air['Passengers'])plt.figure(figsize=(12, 5))
plt.plot(df_air['log_passengers'])
plt.title('Log-Transformed Airline Passengers')
plt.grid(True)
plt.show()# 季节性差分 (m=12)
df_air['log_seasonal_diff'] = df_air['log_passengers'].diff(12)# 丢掉因差分产生的NaN值
deseasonalized_data = df_air['log_seasonal_diff'].dropna()plt.figure(figsize=(12, 5))
plt.plot(deseasonalized_data)
plt.title('Deseasonalized and Log-Transformed Data')
plt.grid(True)
plt.show()# 对去季节性的数据再进行一阶差分，使其平稳
stationary_data = deseasonalized_data.diff(1).dropna()# 绘制ACF和PACF图
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
plot_acf(stationary_data, lags=24, ax=ax1)
plot_pacf(stationary_data, lags=24, ax=ax2)
plt.show()from statsmodels.tsa.arima.model import ARIMA
# 准备训练数据（经过预处理的）
# 注意：我们是对去季节性的数据进行建模
train_data = deseasonalized_data[:'1959']
# 定义并拟合ARIMA(1,1,1)模型
model = ARIMA(train_data, order=(1, 1, 1))
model_fit = model.fit()
print(model_fit.summary())
# 预测未来12个点（1960年全年）
# 模型预测的是经过变换后的数据
forecast_diff = model_fit.forecast(steps=12)发现这个模型的显著性还是不够，可以根据ma值推断我们的q应该是0，就重新构建ARIMA模型
model_simplified = ARIMA(train_data, order=(1, 1, 0))
model_fit_simplified = model_simplified.fit()
print(model_fit_simplified.summary())处理完后的数据预测，要重新把“残差”加回去得到原来的预测结果
# forecast_deseasonalized 是对“去季节性”序列的直接预测
forecast_deseasonalized = model_fit.forecast(steps=12)
# 逆向变换过程
predictions = []
# 逐步还原预测值
for i in range(len(forecast_deseasonalized)):# 1. 还原季节性差分 (直接使用预测值)# log(y_t) = y'_t + log(y_{t-12})# y'_t 就是 forecast_deseasonalized[i]# 获取12个月前的历史对数值last_year_log_val = df_air['log_passengers']['1959'].iloc[i]# 加上历史值，得到预测的对数值pred_log = forecast_deseasonalized.iloc[i] + last_year_log_val# 2. 还原对数变换 (exp)pred_original_scale = np.exp(pred_log)predictions.append(pred_original_scale)
# 将预测结果转换为Series，方便绘图
# 注意：forecast_deseasonalized 自带了正确的日期索引，可以直接使用
predictions_series = pd.Series(predictions, index=forecast_deseasonalized.index)
# --- 重新绘图 ---
plt.figure(figsize=(14, 7))
# 原始数据
plt.plot(df_air['Passengers'], label='Original Data')
# 我们的预测
plt.plot(predictions_series, color='red', linestyle='--', label='ARIMA Forecast (Corrected)')
plt.title('Airline Passengers Forecast using Manual Preprocessing + ARIMA')
plt.xlabel('Year')
plt.ylabel('Number of Passengers')
plt.legend()
plt.grid(True)
plt.show()