DAY 58 经典时序预测模型2

知识点回顾：

时序建模的流程
时序任务经典单变量数据集
ARIMA（p，d，q）模型实战
SARIMA摘要图的理解
处理不平稳的2种差分
1. n阶差分---处理趋势
2. 季节性差分---处理季节性

建立一个ARIMA模型，通常遵循以下步骤：

1. 数据可视化：观察原始时间序列图，判断是否存在趋势或季节性。

2. 平稳性检验：

- 对原始序列进行ADF检验。

- 如果p值 > 0.05，说明序列非平稳，需要进行差分。

3. 确定差分次数 d:

- 进行一阶差分，然后再次进行ADF检验。

- 如果平稳了，则 d=1。否则，继续差分，直到平稳。

4. 确定 p 和 q:

- 对差分后的平稳序列绘制ACF和PACF图。

- 根据昨天学习的规则（PACF定p，ACF定q）来选择p和q的值。

5. 建立并训练ARIMA(p, d, q)模型。

模型评估与诊断：查看模型的摘要信息，检查残差是否为白噪声。
1. AIC用来对比不同模型选择，越小越好
2. 关注系数是否显著，通过p值或者置信区间均可
3. 残差的白噪声检验+正态分布检验

7. 进行预测（需要还原回差分前的结构）

作业：对太阳黑子数量数据集用arima完成流程

1.加载太阳黑子数据集，并绘制原始数据图像

import matplotlib.pyplot as plt
from statsmodels.datasets import sunspots
from statsmodels.tsa.stattools import adfuller
import pandas as pd
from statsmodels.tsa.arima.model import ARIMA
# 设置中文显示
plt.rcParams["font.family"] = ["SimHei", "WenQuanYi Micro Hei", "Heiti TC"]
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题# 加载太阳黑子数据集
df = sunspots.load_pandas().data
df.columns = ['YEAR', 'SUNACTIVITY']  # 设置列名
df['YEAR'] = pd.date_range(start='1700', end='2009', freq='A')  # 修正结束年份为2008
df.set_index('YEAR', inplace=True)  # 设置年份为索引# 创建画布和子图
plt.figure(figsize=(12, 6))# 绘制太阳黑子活动随时间变化的折线图
plt.plot(df.index, df['SUNACTIVITY'], 'b-', linewidth=1.5)# 添加标题和标签
plt.title('太阳黑子活动随时间变化趋势图 (1700-2008)', fontsize=16)
plt.xlabel('年份', fontsize=14)
plt.ylabel('太阳黑子数量', fontsize=14)# 美化图表
plt.tight_layout()  # 自动调整布局
plt.show()

无明显趋势，但有强周期性：大约每11年一个周期

2.检验平稳性

# 进行ADF检验
result = adfuller(df['SUNACTIVITY'])
print('ADF检验结果:')
print(f'ADF统计量: {result[0]}')
print(f'p值: {result[1]}')
print('临界值:')
for key, value in result[4].items():print(f'\t{key}: {value}')if result[1] <= 0.05:print("强烈拒绝原假设，时间序列是平稳的")
else:print("接受原假设，时间序列存在单位根，是非平稳的")

3.一阶差分

# 进行一阶差分
df['SUNACTIVITY_DIFF'] = df['SUNACTIVITY'].diff()
# 删除差分后产生的NaN值
df_diff = df.dropna()
# 绘制一阶差分后的时间序列图
plt.plot(df_diff.index, df_diff['SUNACTIVITY_DIFF'], 'r-', linewidth=1.5)
plt.title('一阶差分后的太阳黑子活动时间序列', fontsize=16)
plt.xlabel('年份', fontsize=14)
plt.ylabel('差分后的值', fontsize=14)# 美化图表
plt.tight_layout()  # 自动调整布局
plt.show()

4.一阶差分后验证平稳性

result_diff = adfuller(df_diff['SUNACTIVITY_DIFF'])
print('一阶差分后数据ADF检验结果:')
print(f'ADF统计量: {result_diff[0]}')
print(f'p值: {result_diff[1]}')
print('临界值:')
for key, value in result_diff[4].items():print(f'\t{key}: {value}')if result_diff[1] <= 0.05:print("强烈拒绝原假设，一阶差分后的时间序列是平稳的")
else:print("接受原假设，一阶差分后的时间序列存在单位根，是非平稳的")

5.绘制ACF和PACF图

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf# 一阶差分后的数据是平稳的，直接绘制ACF和PACF
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))plot_acf(df_diff['SUNACTIVITY_DIFF'], ax=ax1, lags=40)
ax1.set_title('自相关函数 (ACF)')plot_pacf(df_diff['SUNACTIVITY_DIFF'], ax=ax2, lags=40)
ax2.set_title('偏自相关函数 (PACF)')plt.tight_layout()
plt.show()

所以p=2或3,d=1,q=2或3

import warnings
warnings.filterwarnings("ignore")
# 建立ARIMA(2, 0, 0)模型
model_sun_1 = ARIMA(df['SUNACTIVITY'], order=(2, 1, 3))
arima_result_sun_1 = model_sun_1.fit()# 打印模型摘要
print(arima_result_sun_1.summary())

6.预测

# 预测部分代码
# 样本内预测（历史数据的拟合）
start_date = '1950-12-31'  # 从1950年开始显示预测结果
end_date = '2008-12-31'    # 预测到数据的最后# 获取样本内预测值和置信区间
pred = arima_result_sun_1.get_prediction(start=pd.to_datetime(start_date), end=pd.to_datetime(end_date), dynamic=False)
pred_ci = pred.conf_int()  # 置信区间# 未来预测（预测未来20年的数据）
forecast_steps = 20  # 预测未来20年
pred_uc = arima_result_sun_1.get_forecast(steps=forecast_steps)
pred_uc_ci = pred_uc.conf_int()  # 预测的置信区间# 为未来预测创建日期索引
last_date = df.index[-1]
future_dates = pd.date_range(start=last_date + pd.DateOffset(years=1), periods=forecast_steps, freq='A')# 创建预测结果的DataFrame
forecast_df = pd.DataFrame({'forecast': pred_uc.predicted_mean.values,'lower_ci': pred_uc_ci.iloc[:, 0].values,'upper_ci': pred_uc_ci.iloc[:, 1].values
}, index=future_dates)# 绘制预测结果
plt.figure(figsize=(14, 7))# 绘制原始数据
plt.plot(df.index, df['SUNACTIVITY'], 'b-', label='历史数据', linewidth=1.5)# 绘制样本内预测
plt.plot(pred.predicted_mean.index, pred.predicted_mean, 'r--', label='样本内预测', linewidth=1.5)
plt.fill_between(pred_ci.index,pred_ci.iloc[:, 0],pred_ci.iloc[:, 1], color='r', alpha=0.1)# 绘制未来预测
plt.plot(forecast_df.index, forecast_df['forecast'], 'g--', label='未来预测', linewidth=1.5)
plt.fill_between(forecast_df.index,forecast_df['lower_ci'],forecast_df['upper_ci'], color='g', alpha=0.1)# 添加标题和标签
plt.title('太阳黑子活动的预测结果', fontsize=16)
plt.xlabel('年份', fontsize=14)
plt.ylabel('太阳黑子数量', fontsize=14)
plt.legend(loc='upper left', fontsize=12)
plt.grid(True, linestyle='--', alpha=0.7)plt.tight_layout()
plt.show()# 打印未来20年的预测结果
print("未来20年太阳黑子活动的预测结果：")
print(forecast_df[['forecast', 'lower_ci', 'upper_ci']])