季节性时间序列预测之SARIMA模型

大家好，欢迎来到IT知识分享网。

该文主要讲解有关季节性时序模型——SARIMA模型，如果想了解ARIMA模型的小伙伴可以去我的这篇博客《时间序列预测之ARIMA模型》

一、基本概念

一些针对SARIMA模型的预备知识，若只需应用SARIMA模型求解可跳至第二节

1、季节性序列

在某些时间序列中，存在明显的周期性变化。这种周期是由于季节性变化（包括季度、月度、周度等变化）或其他一些固有因素（天气、节假日）引起的。这类序列可以统称为季节性序列。

2、SARIMA模型

SARIMA（Seansonal ARIMA）可以支持带有季节性成分特点的时间序列数据。它在ARIMA模型的基础上增加了4个季节性参数，分别是3个超参数（P,D,Q）和季节性周期参数s。SARIMA模型的参数概念与特点如下表所示：

$SARIMA(p,d,q)(P,D,Q{)}_S$模型具体公式如下
$${\varphi }_p(B){\Phi }_P(B_s)(1-B_s{)}^d(1-B_s{)}^D{x}_t={\theta }_q(B){\Theta }_Q(B_s){ϵ}_t$$是不是感觉很复杂（反正博主第一次学的时候是一头雾水），下面我们拆开讲讲公式中各个符号的含义与作用:

• 延迟算子$B$
  延迟算子类似于一个时间指针，当前序列值乘以一个延迟算子，就相当于把当前序列值的时间向过去拨了一个时刻$$x_{t-p}=B^p{x}_t$$特别地，$(1-B)x_t=x_t-x_{t-1}$表示对$x_t$做一阶差分，以此类推$(1-B{)}^d{x}_t$表示对$x_t$做d阶差分。
  
• 季节性延迟算子$B_s$
  该延迟算子相当于把当前序列值的时间向过去回拨了一个周期$$x_{t-s\ast P}=B_s^P{x}_t$$特别地，$(1-B_s)x_t=x_t-x_{t-s}$可以通俗理解为对$x_t$做了一次周期为s的季节性差分，以此类推$(1-B_s{)}^D{x}_t$表示对$x_t$做D次周期为s的季节性差分
  。
  
  
• 延迟多项式算子$\varphi (B)$与$\theta (B)$
  $$\varphi (B)=1-{\varphi }_{1}B-\dots -{\varphi }_p{B}^p$$$$\theta (B)=1-{\theta }_{1}B-\dots -{\theta }_q{B}^q$$如AR模型可以简化为$\varphi (B)x_t={ϵ}_t$，同理ARMA模型可以简化为$\varphi (B)x_t=\theta (B){ϵ}_t$，ARIMA模型可以简化为$\varphi (B)(1-B{)}^d{x}_t=\theta (B){ϵ}_t$
  
• 季节性延迟多项式算子$\Phi (B_s)$与$\Theta (B_s)$
  $$\Phi (B_s)=1-{\Phi }_1{B}_s-\dots -{\Phi }_p{B}_s^P$$$$\Theta (B_s)=1-{\Theta }_1{B}_s-\dots -{\Theta }_q{B}_s^Q$$
  

所以，这个看似复杂的SARIMA模型其实就是对一个具有季节性特点的时间序列{
$x_t$}做D次季节性差分（去周期）和d次常差分（去趋势）得到一个新的序列{
$y_t$}，这段话用数学公式表示为：$y_t=(1-B{)}^d(1-B_s{)}^D{x}_t$。然后对差分后的{
$y_t$}建立如下模型：
$$\begin{gathered} y_t={\varphi }_1{x}_{t-1}+{\varphi }_2{x}_{t-2}+\dots +{\varphi }_p{x}_{t-p} \\ +{\Phi }_1{x}_{t-s}+{\Phi }_2{x}_{t-2s}+\dots {\Phi }_P{x}_{t-P\ast s} \\ +{\theta }_1{ϵ}_{t-1}+{\theta }_2{ϵ}_{t-2}+\dots +{\theta }_q{ϵ}_{t-q} \\ +{\Theta }_1{ϵ}_{t-s}+{\Theta }_2{ϵ}_{t-2s}+\dots +{\Theta }_Q{ϵ}_{t-Q\ast s} \\ +{ϵ}_t \end{gathered}$$这里博主举个栗子给大家加深理解（如果栗子有质量问题，欢迎大家在评论区中留言讨论），假如需要对某品牌羽绒服12月份的销量进行预测（假定羽绒服的月销量是一个周期为12的时间序列），那么上式中：

• 第一行表示的是12月份羽绒服的销量与今年各个月份销量的联系
• 第二行表示今年12月份羽绒服的销量与去年、前年…的12月份销量的联系
• 第三行表示今年内某时间点上的突发事件（比如代言人偷税漏税被抓了或者产品质量被曝欺骗消费者等等）对12月份羽绒服销量的影响
• 第四行表示去年、前年…12月份的突发事件对今年12月份羽绒服的销量的影响

二、SARIMA模型求解

1、数据准备

1）导入库

#准备工作 import numpy as np import pandas as pd import statsmodels.api as sm import matplotlib.pyplot as plt from statsmodels.graphics.tsaplots import plot_acf,plot_pacf import itertools from sklearn.metrics import r2_score as rs import warnings warnings.filterwarnings("ignore")#忽略输出警告 plt.rcParams["font.sans-serif"]=["SimHei"]#用来正常显示中文标签 plt.rcParams["axes.unicode_minus"]=False#用来显示负号 %matplotlib inline #在jupyter中显示figure，该语句只在jupyter上有用 

2）读取数据
本文使用的是每月发电产生的二氧化碳排放量的公共数据集，数据时间跨度从1973年到2020年。

df=pd.read_excel("E:\\代码学习\\CSDN博客\\SARIMA模型\\Carbon_Dioxide_Emissions_From_Energy_Consumption-_Electric_Power_Sector.xlsx"\ ,index_col="Month")#指定Month列作为索引列 df 

3）数据预处理
这里介绍几种python中数据简单预处理的函数

• dropna()函数
  该函数能找到DataFrame类型数据的空值，将空值所在行/列删除后将新的DataFrame作为返回值返回
  

#数据预处理 #dropna() #该函数能找到DataFrame类型数据的空值，将空值所在行/列删除后将新的DataFrame作为返回值返回 #参数axis=0表示按行删除，1表示按列删除 #参数how='any'表示该行/列只要有一个以上的空值，就删除该行/列；'all'表示该行/列全部都为空值，就删除该行/列 df=df.dropna(axis=0, how='any') 

• fillna()函数
  该函数能够使用指定的方法填充NA/NaN值
  

#fillna() #能够使用指定的方法填充NA/NaN值 #参数value，用于填充空值的值 #参数method='fill'表示用前面行/列的值，填充当前行/列的空值；'bfill'表示用后面行/列的值，填充当前行/列的空值 #参数axis=0表示按行，1表示按列 #参数limit表示如果method被指定，对于连续的空值，这段连续区域，最多填充前limit个空值 df=df.fillna(method='ffill',axis=0) df 

• drop()函数
  该函数能够去除空值
  

#drop() #去除空值 #利用key_list将df分为df1和df2 key_list=["Geothermal Energy Electric Power Sector CO2 Emissions",\ "Non-Biomass Waste Electric Power Sector CO2 Emissions"] #df1=df[df.keys().drop(key_list)]是ppt上的代码，通过去除指定列的keys值来达到筛选的目的 df1=df.drop(key_list,axis=1) df2=df[key_list] NV=df2.values=="Not Available"#判断df2的数据是否为Not Available，返回值是bool类型 print(NV) index=df2[NV].index#获得"Not Available"的所在行名 print(index) df2=df2.drop(labels=index)#参数labels:待删除的行名/列名 df2 

4）数据可视化
对天然气的CO2数据进行可视化分析

#可视化处理 #天然气CO2排放量 NGE=df1["Natural Gas Electric Power Sector CO2 Emissions"] NGE.head() #折线图 fig, ax=plt.subplots(figsize=(15,15)) NGE.plot(ax=ax,fontsize=15) ax.set_title("天然气碳排放",fontsize=25) ax.set_xlabel("时间（月）",fontsize=25) ax.set_ylabel("碳排放量(百万总吨)",fontsize=25) ax.legend(loc="best",fontsize=15) ax.grid() 

5）分解时序
STL（Seasonal and Trend decomposition using Loess）是一个非常通用和稳健强硬的分解时间序列的方法

#分解时序 #STL（Seasonal and Trend decomposition using Loess）是一个非常通用和稳健强硬的分解时间序列的方法 import statsmodels.api as sm #decompostion=tsa.STL(NGE).fit()报错，这里前面加上索引sm decompostion=sm.tsa.STL(NGE).fit()#statsmodels.tsa.api:时间序列模型和方法 decompostion.plot() #趋势效益 trend=decompostion.trend #季节效应 seasonal=decompostion.seasonal #随机效应 residual=decompostion.resid 

2、平稳性检验

使用ADF平稳性检验方法检验天然气CO2排放量数据的平稳性

#平稳性检验 #自定义函数用于ADF检查平稳性 from statsmodels.tsa.stattools import adfuller as ADF def test_stationarity(timeseries,alpha):#alpha为检验选取的显著性水平 adf=ADF(timeseries) p=adf[1]#p值 critical_value=adf[4]["5%"]#在95%置信区间下的临界的ADF检验值 test_statistic=adf[0]#ADF统计量 if p<alpha and test_statistic<critical_value: print("ADF平稳性检验结果：在显著性水平%s下，数据经检验平稳"%alpha) return True else: print("ADF平稳性检验结果：在显著性水平%s下，数据经检验不平稳"%alpha) return False 

得到的结果为：“ADF平稳性检验结果：在显著性水平0.001下，数据经检验不平稳”

#原始数据平稳性检验 test_stationarity(NGE,1e-3) 

由于原始数据经检验不平稳，需要将不平稳数据化为平稳数据。因此需要对天然气CO2排放量数据作一阶常差分去除趋势（在STL时序分解的第二个子图中可以看出原始数据有明显的上升趋势）和十二步季节性差分去除季节性（由时序分解图和原始图像可以看出原始数据是周期为12的季节性时序数据）。最后将处理后的数据进行平稳性检验，得到的结果为：“ADF平稳性检验结果：在显著性水平0.001下，数据经检验平稳”。

#将数据化为平稳数据 #一阶差分 NGE_diff1=NGE.diff(1) #十二步差分 NGE_seasonal=NGE_diff1.diff(12)#非平稳序列经过d阶常差分和D阶季节差分变为平稳时间序列 print(NGE_seasonal) #十二步季节差分平稳性检验结果 test_stationarity(NGE_seasonal.dropna(),1e-3)#使用dropna()去除NaN值 

3、白噪声检验

在得到平稳数据后，需要对平稳数据进行白噪声检验，验证序列中有用信息是否已经被提取完毕。若序列是白噪声序列，说明序列中有用信息已经被提取完，只剩随机扰动。白噪声序列没有分析价值，应该被舍弃。
本文使用Ljung-Box检验白噪声，即LB白噪声检验方法。通过statsmodel库调用acorr_ljungbox()函数实现LB白噪声检验。

#LB白噪声检验 from statsmodels.stats.diagnostic import acorr_ljungbox def test_white_noise(data,alpha): [[lb],[p]]=acorr_ljungbox(data,lags=1) if p<alpha: print('LB白噪声检验结果：在显著性水平%s下，数据经检验为非白噪声序列'%alpha) else: print('LB白噪声检验结果：在显著性水平%s下，数据经检验为白噪声序列'%alpha) 

白噪声检验的结果为：“LB白噪声检验结果：在显著性水平0.01下，数据经检验为非白噪声序列”。因此，该数据为平稳非白噪声序列，可以继续分析。

#白噪声检验 test_white_noise(NGE_seasonal.dropna(),0.01) 

4、模型定阶

采用网格搜索法定阶（虽然图解法也能定阶，但此方法适用性不强）。

#搜索法定阶 def SARIMA_search(data): p=q=range(0,3) s=[12]#周期为12 d=[1]#做了一次季节性差分 PDQs=list(itertools.product(p,d,q,s))#itertools.product()得到的是可迭代对象的笛卡儿积 pdq=list(itertools.product(p,d,q))#list是python中是序列数据结构，序列中的每个元素都分配一个数字定位位置 params=[] seasonal_params=[] results=[] grid=pd.DataFrame() for param in pdq: for seasonal_param in PDQs: #建立模型 mod= sm.tsa.SARIMAX(data,order=param,seasonal_order=seasonal_param,\ enforce_stationarity=False, enforce_invertibility=False) #实现数据在模型中训练 result=mod.fit() print("ARIMA{}x{}-AIC:{}".format(param,seasonal_param,result.aic)) #format表示python格式化输出，使用{}代替% params.append(param) seasonal_params.append(seasonal_param) results.append(result.aic) grid["pdq"]=params grid["PDQs"]=seasonal_params grid["aic"]=results print(grid[grid["aic"]==grid["aic"].min()]) SARIMA_search(NGE_seasonal.dropna()) 

采用AIC最小原则方法定阶的结果如图所示，定阶模型为$SARIMA(1,1,2)(0,1,2{)}_{12}$

5、模型的建立与检验

建立并训练$SARIMA(1,1,2)(0,1,2{)}_{12}$模型

#建立模型 model=sm.tsa.SARIMAX(NGE,order=(1,1,2),seasonal_order=(0,1,2,12)) SARIMA_m=model.fit() print(SARIMA_m.summary()) 

对建立的模型进行白噪声检验并使用plot_diagnostics()函数画出检验图像，得到的结果为：“LB白噪声检验结果：在显著性水平0.05下，数据经检验为白噪声序列”。这说明该模型已经有效提取了数据信息，且观察图像成正态分布，模型通过检验。

#模型检验 test_white_noise(SARIMA_m.resid,0.05)#SARIMA_m.resid提取模型残差，并检验是否为白噪声 fig=SARIMA_m.plot_diagnostics(figsize=(15,12))#plot_diagnostics对象允许我们快速生成模型诊断并调查任何异常行为 

6、模型预测

#模型预测 from sklearn.metrics import mean_squared_error as mse from sklearn.metrics import mean_absolute_error as mae #获取预测结果，自定义预测误差 def PredictionAnalysis(data,model,start,dynamic=False): pred=model.get_prediction(start=start,dynamic=dynamic,full_results=True) pci=pred.conf_int()#置信区间 pm=pred.predicted_mean#预测值 truth=data[start:]#真实值 pc=pd.concat([truth,pm,pci],axis=1)#按列拼接 pc.columns=['true','pred','up','low']#定义列索引 print("1、MSE:{}".format(mse(truth,pm))) print("2、RMSE:{}".format(np.sqrt(mse(truth,pm)))) print("3、MAE:{}".format(mae(truth,pm))) return pc #绘制预测结果 def PredictonPlot(pc): plt.figure(figsize=(10,8)) plt.fill_between(pc.index,pc['up'],pc['low'],color='grey',\ alpha=0.15,label='confidence interval')#画出置信区间 plt.plot(pc['true'],label='base data') plt.plot(pc['pred'],label='prediction curve') plt.legend() plt.show return True 

• 静态预测：进行一系列的一步预测，即它必须用真实值来进行预测

#静态预测：进行一系列的一步预测，即它必须用真实值来进行预测 pred=PredictionAnalysis(NGE,SARIMA_m,'2015-01-01') PredictonPlot(pred) 

• 动态预测：进行多步预测，除了第一个预测值是用实际值预测外，其后各预测值都是采用递推预测

#动态预测：进行多步预测，除了第一个预测值是用实际值预测外，其后各预测值都是采用递推预测 pred=PredictionAnalysis(NGE,SARIMA_m,'2015-01-01',dynamic=True) PredictonPlot(pred) 

预测未来六十天的变化

#预测未来 forecast=SARIMA_m.get_forecast(steps=60) #预测整体可视化 fig,ax=plt.subplots(figsize=(20,16)) NGE.plot(ax=ax,label="base data") forecast.predicted_mean.plot(ax=ax,label="forecast data") #ax.fill_between(forecast.conf_int().index(),forecast.conf_int().iloc[:,0],\ # forecast.conf_int().iloc[:,1],color='grey',alpha=0.15,label='confidence interval') ax.legend(loc="best",fontsize=20) ax.set_xlabel("时间（年）",fontsize=20) ax.set_ylabel("天然气二氧化碳排放水平",fontsize=18) plt.xticks(fontsize=20) plt.yticks(fontsize=20) plt.show