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']])

@浙大疏錦行