用Python实战SARIMA模型：手把手教你预测月度用电碳排放（附完整代码）

Python实战SARIMA模型：从数据清洗到碳排放预测全流程解析

当企业需要制定碳中和战略时，准确预测未来碳排放量成为关键决策依据。某能源集团的数据分析师王敏最近就遇到了这样的挑战：管理层要求她基于历史数据，预测未来两年集团电力生产的月度碳排放趋势。传统方法难以捕捉季节性波动，而SARIMA模型恰好能解决这个问题。

1. 环境准备与数据加载

工欲善其事，必先利其器。我们首先配置Python环境，安装必要的库：

# 基础数据处理库 import pandas as pd import numpy as np # 统计分析库 import statsmodels.api as sm from statsmodels.tsa.statespace.sarimax import SARIMAX from statsmodels.tsa.seasonal import seasonal_decompose # 可视化库 import matplotlib.pyplot as plt import seaborn as sns # 模型评估 from sklearn.metrics import mean_absolute_error, mean_squared_error # 忽略警告信息 import warnings warnings.filterwarnings("ignore")

加载碳排放数据集时，需要特别注意数据质量。真实业务数据往往存在以下问题：

• 时间戳格式不统一
• 异常值或缺失值
• 计量单位不一致

# 加载数据示例 df = pd.read_csv('power_emission.csv', parse_dates=['month'], index_col='month') # 检查数据前5行 print(df.head()) # 检查缺失值 print(df.isnull().sum())

常见数据问题处理方案：

2. 数据探索与平稳性处理

高质量的数据可视化能帮助我们直观理解数据特征。以下是关键可视化步骤：

# 绘制原始数据趋势 plt.figure(figsize=(12,6)) df['emission'].plot(title='月度碳排放趋势(1973-2020)') plt.ylabel('百万吨CO2') plt.grid(True) plt.show()

通过STL分解观察数据的季节性、趋势和残差分量：

# 季节性分解 decomposition = seasonal_decompose(df['emission'], model='additive', period=12) decomposition.plot() plt.tight_layout() plt.show()

平稳性检验是时间序列分析的关键步骤。我们使用ADF检验：

from statsmodels.tsa.stattools import adfuller def adf_test(series): result = adfuller(series.dropna()) print('ADF统计量: %f' % result[0]) print('p值: %f' % result[1]) print('临界值:') for key, value in result[4].items(): print('\t%s: %.3f' % (key, value)) if result[1] < 0.05: print("拒绝原假设，数据平稳") else: print("无法拒绝原假设，数据非平稳") adf_test(df['emission'])

当数据不平稳时，我们需要进行差分处理：

# 一阶差分去趋势 df['diff_1'] = df['emission'].diff(1) # 季节性差分(周期12个月) df['diff_seasonal'] = df['diff_1'].diff(12) # 再次检验平稳性 adf_test(df['diff_seasonal'].dropna())

3. 模型构建与参数优化

SARIMA模型有7个关键参数：(p,d,q)(P,D,Q)m。确定这些参数的最佳组合是建模的核心挑战。

参数网格搜索实现：

# 定义参数搜索空间 p = d = q = range(0, 2) P = D = Q = range(0, 2) m = 12 # 月度数据的季节周期 # 生成所有参数组合 pdq = list(itertools.product(p, d, q)) seasonal_pdq = list(itertools.product(P, D, Q, [m])) # 网格搜索寻找最优参数 best_aic = float("inf") best_params = None for param in pdq: for param_seasonal in seasonal_pdq: try: mod = SARIMAX(df['emission'], order=param, seasonal_order=param_seasonal, enforce_stationarity=False, enforce_invertibility=False) results = mod.fit() if results.aic < best_aic: best_aic = results.aic best_params = (param, param_seasonal) print(f'SARIMA{param}x{param_seasonal} - AIC:{results.aic:.2f}') except: continue print(f'\n最优参数组合: {best_params} - AIC: {best_aic:.2f}')

参数选择经验法则：

1. 观察ACF/PACF图确定初步参数范围
2. 优先尝试d+D≤2的组合
3. 季节性参数通常不超过1阶
4. 权衡模型复杂度(AIC)与过拟合风险

4. 模型训练与验证

确定最优参数后，我们训练最终模型：

# 使用最优参数训练模型 best_order, best_seasonal_order = best_params model = SARIMAX(df['emission'], order=best_order, seasonal_order=best_seasonal_order, enforce_stationarity=False) results = model.fit() # 输出模型摘要 print(results.summary())

模型诊断要点：

• 残差应近似白噪声
• Q-Q图应接近直线
• 残差自相关函数(ACF)无显著相关性

# 模型诊断图 results.plot_diagnostics(figsize=(12,8)) plt.tight_layout() plt.show()

验证模型预测能力时，我们保留最后24个月作为测试集：

# 划分训练测试集 train = df.iloc[:-24] test = df.iloc[-24:] # 在训练集上重新训练模型 model = SARIMAX(train['emission'], order=best_order, seasonal_order=best_seasonal_order) fitted = model.fit() # 预测测试集 forecast = fitted.get_forecast(steps=24) forecast_ci = forecast.conf_int() # 可视化预测结果 plt.figure(figsize=(12,6)) plt.plot(train.index, train['emission'], label='训练数据') plt.plot(test.index, test['emission'], label='实际值') plt.plot(test.index, forecast.predicted_mean, label='预测值') plt.fill_between(test.index, forecast_ci.iloc[:,0], forecast_ci.iloc[:,1], color='gray', alpha=0.2) plt.title('SARIMA模型预测效果验证') plt.legend() plt.show()

评估指标计算：

# 计算评估指标 mae = mean_absolute_error(test['emission'], forecast.predicted_mean) rmse = np.sqrt(mean_squared_error(test['emission'], forecast.predicted_mean)) print(f'MAE: {mae:.2f}') print(f'RMSE: {rmse:.2f}')

5. 模型部署与生产应用

将训练好的模型应用于实际业务预测：

# 全量数据重新训练 final_model = SARIMAX(df['emission'], order=best_order, seasonal_order=best_seasonal_order) final_results = final_model.fit() # 预测未来24个月 forecast = final_results.get_forecast(steps=24) forecast_ci = forecast.conf_int() # 可视化长期预测 plt.figure(figsize=(12,6)) plt.plot(df.index, df['emission'], label='历史数据') plt.plot(pd.date_range(df.index[-1], periods=25, freq='M')[1:], forecast.predicted_mean, label='未来预测') plt.fill_between(pd.date_range(df.index[-1], periods=25, freq='M')[1:], forecast_ci.iloc[:,0], forecast_ci.iloc[:,1], color='gray', alpha=0.2) plt.title('未来两年碳排放预测') plt.ylabel('百万吨CO2') plt.legend() plt.grid(True) plt.show()

生产环境部署建议：

1. 使用Joblib或Pickle保存训练好的模型
2. 设置定期(如每月)模型重训练机制
3. 实现自动化预测结果推送
4. 建立模型性能监控体系

# 模型保存示例 import joblib joblib.dump(final_results, 'sarima_emission_model.pkl') # 模型加载示例 loaded_model = joblib.load('sarima_emission_model.pkl') new_forecast = loaded_model.get_forecast(steps=12)

6. 模型优化与高级技巧

基���SARIMA模型可以进一步优化提升预测精度：

1. 外生变量引入

当有其他影响因素数据时，可以使用SARIMAX模型：

# 假设有温度数据作为外生变量 exog = pd.read_csv('temperature.csv', index_col='month', parse_dates=True) model = SARIMAX(df['emission'], exog=exog, order=(1,1,1), seasonal_order=(1,1,1,12)) results = model.fit()

2. 参数自动优化

使用pmdarima库实现自动参数选择：

from pmdarima import auto_arima model = auto_arima(df['emission'], seasonal=True, m=12, trace=True, error_action='ignore', suppress_warnings=True) print(model.summary())

3. 预测区间调整

根据业务需求调整置信区间：

# 获取不同置信水平的预测区间 forecast_95 = final_results.get_forecast(steps=24).conf_int(alpha=0.05) forecast_80 = final_results.get_forecast(steps=24).conf_int(alpha=0.2)

4. 多周期预测比较

评估不同预测周期下的模型表现：

7. 业务应用与决策支持

将模型预测结果转化为业务洞察是关键。以下是典型应用场景：

1. 碳配额规划

基于预测结果制定碳配额采购计划，避免超额排放罚款或配额浪费。

2. 减排措施评估

模拟不同减排措施实施后的预测曲线变化，评估措施效果。

3. 能源结构调整

分析不同能源占比变化对碳排放的影响，优化能源结构。

4. 报告自动化

将预测结果自动生成可视化报告，支持管理层决策。

# 生成预测报告示例 report_data = { '当前排放水平': df['emission'][-1], '下季度预测': forecast.predicted_mean[:3].mean(), '明年同期变化率': (forecast.predicted_mean[12]/df['emission'][-12]-1)*100 } pd.DataFrame.from_dict(report_data, orient='index', columns=['值'])

实际项目中，我们曾遇到一个典型案例：某电厂通过SARIMA模型预测发现，如果不采取改进措施，明年三季度将超出碳配额7.2%。基于这一预警，他们提前实施了能效提升计划，最终避免了约280万元的超额排放罚款。