用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_colmonth) # 检查数据前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], modeladditive, period12) 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], orderparam, seasonal_orderparam_seasonal, enforce_stationarityFalse, enforce_invertibilityFalse) results mod.fit() if results.aic best_aic: best_aic results.aic best_params (param, param_seasonal) print(fSARIMA{param}x{param_seasonal} - AIC:{results.aic:.2f}) except: continue print(f\n最优参数组合: {best_params} - AIC: {best_aic:.2f})参数选择经验法则观察ACF/PACF图确定初步参数范围优先尝试dD≤2的组合季节性参数通常不超过1阶权衡模型复杂度(AIC)与过拟合风险4. 模型训练与验证确定最优参数后我们训练最终模型# 使用最优参数训练模型 best_order, best_seasonal_order best_params model SARIMAX(df[emission], orderbest_order, seasonal_orderbest_seasonal_order, enforce_stationarityFalse) 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], orderbest_order, seasonal_orderbest_seasonal_order) fitted model.fit() # 预测测试集 forecast fitted.get_forecast(steps24) 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], colorgray, alpha0.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(fMAE: {mae:.2f}) print(fRMSE: {rmse:.2f})5. 模型部署与生产应用将训练好的模型应用于实际业务预测# 全量数据重新训练 final_model SARIMAX(df[emission], orderbest_order, seasonal_orderbest_seasonal_order) final_results final_model.fit() # 预测未来24个月 forecast final_results.get_forecast(steps24) 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], periods25, freqM)[1:], forecast.predicted_mean, label未来预测) plt.fill_between(pd.date_range(df.index[-1], periods25, freqM)[1:], forecast_ci.iloc[:,0], forecast_ci.iloc[:,1], colorgray, alpha0.2) plt.title(未来两年碳排放预测) plt.ylabel(百万吨CO2) plt.legend() plt.grid(True) plt.show()生产环境部署建议使用Joblib或Pickle保存训练好的模型设置定期(如每月)模型重训练机制实现自动化预测结果推送建立模型性能监控体系# 模型保存示例 import joblib joblib.dump(final_results, sarima_emission_model.pkl) # 模型加载示例 loaded_model joblib.load(sarima_emission_model.pkl) new_forecast loaded_model.get_forecast(steps12)6. 模型优化与高级技巧基SARIMA模型可以进一步优化提升预测精度1. 外生变量引入当有其他影响因素数据时可以使用SARIMAX模型# 假设有温度数据作为外生变量 exog pd.read_csv(temperature.csv, index_colmonth, parse_datesTrue) model SARIMAX(df[emission], exogexog, 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], seasonalTrue, m12, traceTrue, error_actionignore, suppress_warningsTrue) print(model.summary())3. 预测区间调整根据业务需求调整置信区间# 获取不同置信水平的预测区间 forecast_95 final_results.get_forecast(steps24).conf_int(alpha0.05) forecast_80 final_results.get_forecast(steps24).conf_int(alpha0.2)4. 多周期预测比较评估不同预测周期下的模型表现预测周期(月)MAERMSE训练时间(s)62.12.815123.24.118245.77.3227. 业务应用与决策支持将模型预测结果转化为业务洞察是关键。以下是典型应用场景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, orientindex, columns[值])实际项目中我们曾遇到一个典型案例某电厂通过SARIMA模型预测发现如果不采取改进措施明年三季度将超出碳配额7.2%。基于这一预警他们提前实施了能效提升计划最终避免了约280万元的超额排放罚款。