2023第十二届“认证杯”A题：太阳黑子周期|数学中国数学建模国际赛（小美赛）| 建模秘籍&文章代码思路大全

铛铛！小秘籍来咯！ 
小秘籍希望大家都能轻松建模呀，认证杯也会持续给大家放送思路滴~
抓紧小秘籍，我们出发吧~
 
来看看认证杯（A题）！

问题重述：
太阳黑子是太阳光球上出现的暂时比周围区域更暗的斑点现象。它们是由于磁通量的浓集而引起的表面温度降低区域，抑制对流而形成的。太阳黑子通常出现在活跃区域内，通常成对出现，具有相反的磁极性。它们的数量根据大约11年的太阳周期而变化。我们需要预测太阳黑子，通常我们需要将结果在月度基础上进行平均。
 
问题一、请预测当前太阳周期和下一个太阳周期的开始和结束时间； 
预测太阳周期的开始和结束时间通常是基于太阳黑子活动的观测数据。我们可以使用时间序列分析方法来预测太阳周期的开始和结束时间。在这里，我们将使用ARIMA（差分整合移动平均自回归模型）作为一种常用的时间序列分析方法。
ARIMA 模型是一种包含自回归（AR）、差分（I，表示整合）、和移动平均（MA）的时间序列分析模型。以下是 ARIMA 模型的一般形式：
 
ARIMA(p, d, q)
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
 
# 包含两列：'Date' 和 'Sunspot_Count'
# 日期列是 datetime 类型
data = pd.read_csv("sunspot_data.csv")
data['Date'] = pd.to_datetime(data['Date'])
data.set_index('Date', inplace=True)
sunspot_series = data['Sunspot_Count']
 
# 绘制时序图
plt.figure(figsize=(10, 6))
plt.plot(sunspot_series, label='Sunspot Count')
plt.title('Sunspot Activity Over Time')
plt.xlabel('Date')
plt.ylabel('Sunspot Count')
plt.legend()
plt.show()
 
# 平稳性检验
# 这里可以进行差分操作，观察是否能够使时序数据变得平稳
diff_sunspot_series = sunspot_series.diff().dropna()
 
# 绘制差分后的时序图
plt.figure(figsize=(10, 6))
plt.plot(diff_sunspot_series, label='Differenced Sunspot Count')
plt.title('Differenced Sunspot Activity Over Time')
plt.xlabel('Date')
plt.ylabel('Differenced Sunspot Count')
plt.legend()
plt.show()
 
# ACF 和 PACF 图形
plot_acf(sunspot_series, lags=20)
plt.title('Autocorrelation Function (ACF)')
plt.show()
 
plot_pacf(sunspot_series, lags=20)
plt.title('Partial Autocorrelation Function (PACF)')
plt.show()
 
# ARIMA 模型拟合
# 根据ACF和PACF的图形，选择合适的p和q值
p = 2 # 示例值，根据实际情况调整
d = 1 # 示例值，根据实际情况调整
q = 2 # 示例值，根据实际情况调整
 
model = ARIMA(sunspot_series, order=(p, d, q))
result = model.fit()
 
# 模型诊断
result.plot_diagnostics(figsize=(12, 8))
plt.show()
 
# 预测未来值
forecast_steps = 12 # 示例值，根据实际情况调整
forecast = result.get_forecast(steps=forecast_steps)
forecast_index = pd.date_range(start=sunspot_series.index[-1], periods=forecast_steps + 1, freq='M')[1:]
 
# 绘制预测结果
plt.figure(figsize=(10, 6))
plt.plot(sunspot_series, label='Historical Sunspot Count')
plt.plot(forecast_index, forecast.predicted_mean, color='red', label='Forecasted Sunspot Count')
plt.fill_between(forecast_index,
                 forecast.conf_int()['lower Sunspot_Count'],
                 forecast.conf_int()['upper Sunspot_Count'],
                 color='red', alpha=0.2)
plt.title('Sunspot Activity Forecast')
plt.xlabel('Date')
plt.ylabel('Sunspot Count')
plt.legend()
plt.show()