文章目录

一、基本概念1、季节性序列2、SARIMA模型

二、SARIMA模型求解1、数据准备2、平稳性检验3、白噪声检验4、模型定阶5、模型的建立与检验6、模型预测

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

一、基本概念

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

1、季节性序列

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

2、SARIMA模型

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

参数概念特点p非季节自回归的阶数对于AR模型，可以借助PACF确定阶数p（ACF：拖尾，PACF：p阶截尾）q非季节移动平均的阶数对于MA模型，可以借助PACF确定阶数q（PACF：拖尾，ACF：q阶截尾）d一步差分次数做了多少次一步差分就是多少P季节自回归的阶数（通常不会大于3）从PACF上推断，如果季节长度为12，看PACF图上滞后12阶、24阶、48阶时的偏自相关系数，如果滞后到24阶时表现显著，那么P等于2。Q季节移动平均的阶数（通常不会大于3）从ACF上推断，如果季节长度为12，看ACF图上滞后12阶、24阶、48阶时的自相关系数，如果滞后到12阶时表现显著，那么Q等于1D季节差分次数（通常不会大于1）一般就是 0 或 1，做了季节差分就是 1S季节/周期长度该季节性序列的周期大小（注意不一定是12）

S

A

R

I

M

A

(

p

,

d

,

q

)

(

P

,

D

,

Q

)

S

SARIMA(p,d,q)(P,D,Q)_S

SARIMA(p,d,q)(P,D,Q)S​模型具体公式如下

ϕ

p

(

B

)

Φ

P

(

B

s

)

(

1

−

B

s

)

d

(

1

−

B

s

)

D

x

t

=

θ

q

(

B

)

Θ

Q

(

B

s

)

ϵ

t

\phi_p(B)\Phi_P(B_s)(1-B_s)^d(1-B_s)^Dx_t=\theta_q(B)\Theta_Q(B_s)\epsilon_t

ϕp​(B)ΦP​(Bs​)(1−Bs​)d(1−Bs​)Dxt​=θq​(B)ΘQ​(Bs​)ϵt​是不是感觉很复杂（反正博主第一次学的时候是一头雾水），下面我们拆开讲讲公式中各个符号的含义与作用:

延迟算子

B

B

B 延迟算子类似于一个时间指针，当前序列值乘以一个延迟算子，就相当于把当前序列值的时间向过去拨了一个时刻

x

t

−

p

=

B

p

x

t

x_{t-p}=B^px_t

xt−p​=Bpxt​特别地，

(

1

−

B

)

x

t

=

x

t

−

x

t

−

1

(1-B)x_t=x_t-x_{t-1}

(1−B)xt​=xt​−xt−1​表示对

x

t

x_t

xt​做一阶差分，以此类推

(

1

−

B

)

d

x

t

(1-B)^dx_t

(1−B)dxt​表示对

x

t

x_t

xt​做d阶差分。季节性延迟算子

B

s

B_s

Bs​ 该延迟算子相当于把当前序列值的时间向过去回拨了一个周期

x

t

−

s

∗

P

=

B

s

P

x

t

x_{t-s*P}=B_s^Px_t

xt−s∗P​=BsP​xt​特别地，

(

1

−

B

s

)

x

t

=

x

t

−

x

t

−

s

(1-B_s)x_t=x_t-x_{t-s}

(1−Bs​)xt​=xt​−xt−s​可以通俗理解为对

x

t

x_t

xt​做了一次周期为s的季节性差分，以此类推

(

1

−

B

s

)

D

x

t

(1-B_s)^Dx_t

(1−Bs​)Dxt​表示对

x

t

x_t

xt​做D次周期为s的季节性差分 。延迟多项式算子

ϕ

(

B

)

\phi(B)

ϕ(B)与

θ

(

B

)

\theta(B)

θ(B)

ϕ

(

B

)

=

1

−

ϕ

1

B

−

.

.

.

−

ϕ

p

B

p

\phi(B)=1-\phi_1B-...-\phi_pB^p

ϕ(B)=1−ϕ1​B−...−ϕp​Bp

θ

(

B

)

=

1

−

θ

1

B

−

.

.

.

−

θ

q

B

q

\theta(B)=1-\theta_1B-...-\theta_qB^q

θ(B)=1−θ1​B−...−θq​Bq如AR模型可以简化为

ϕ

(

B

)

x

t

=

ϵ

t

\phi(B)x_t=\epsilon_t

ϕ(B)xt​=ϵt​，同理ARMA模型可以简化为

ϕ

(

B

)

x

t

=

θ

(

B

)

ϵ

t

\phi(B)x_t=\theta(B)\epsilon_t

ϕ(B)xt​=θ(B)ϵt​，ARIMA模型可以简化为

ϕ

(

B

)

(

1

−

B

)

d

x

t

=

θ

(

B

)

ϵ

t

\phi(B)(1-B)^dx_t=\theta(B)\epsilon_t

ϕ(B)(1−B)dxt​=θ(B)ϵt​季节性延迟多项式算子

Φ

(

B

s

)

\Phi(B_s)

Φ(Bs​)与

Θ

(

B

s

)

\Theta(B_s)

Θ(Bs​)

Φ

(

B

s

)

=

1

−

Φ

1

B

s

−

.

.

.

−

Φ

p

B

s

P

\Phi(B_s)=1-\Phi_1B_s-...-\Phi_pB_s^P

Φ(Bs​)=1−Φ1​Bs​−...−Φp​BsP​

Θ

(

B

s

)

=

1

−

Θ

1

B

s

−

.

.

.

−

Θ

q

B

s

Q

\Theta(B_s)=1-\Theta_1B_s-...-\Theta_qB^Q_s

Θ(Bs​)=1−Θ1​Bs​−...−Θq​BsQ​

所以，这个看似复杂的SARIMA模型其实就是对一个具有季节性特点的时间序列{

x

t

x_t

xt​}做D次季节性差分（去周期）和d次常差分（去趋势）得到一个新的序列{

y

t

y_t

yt​}，这段话用数学公式表示为：

y

t

=

(

1

−

B

)

d

(

1

−

B

s

)

D

x

t

y_t=(1-B)^d(1-B_s)^Dx_t

yt​=(1−B)d(1−Bs​)Dxt​。然后对差分后的{

y

t

y_t

yt​}建立如下模型：

y

t

=

ϕ

1

x

t

−

1

+

ϕ

2

x

t

−

2

+

.

.

.

+

ϕ

p

x

t

−

p

+

Φ

1

x

t

−

s

+

Φ

2

x

t

−

2

s

+

.

.

.

Φ

P

x

t

−

P

∗

s

+

θ

1

ϵ

t

−

1

+

θ

2

ϵ

t

−

2

+

.

.

.

+

θ

q

ϵ

t

−

q

+

Θ

1

ϵ

t

−

s

+

Θ

2

ϵ

t

−

2

s

+

.

.

.

+

Θ

Q

ϵ

t

−

Q

∗

s

+

ϵ

t

y_t=\phi_1x_{t-1}+\phi_2x_{t-2}+...+\phi_px_{t-p}\\ +\Phi_1x_{t-s}+\Phi_2x_{t-2s}+...\Phi_Px_{t-P*s}\\ +\theta_1\epsilon_{t-1}+\theta_2\epsilon_{t-2}+...+\theta_q\epsilon_{t-q}\\ +\Theta_1\epsilon_{t-s}+\Theta_2\epsilon_{t-2s}+...+\Theta_Q\epsilon_{t-Q*s}\\ +\epsilon_t

yt​=ϕ1​xt−1​+ϕ2​xt−2​+...+ϕp​xt−p​+Φ1​xt−s​+Φ2​xt−2s​+...ΦP​xt−P∗s​+θ1​ϵt−1​+θ2​ϵt−2​+...+θq​ϵt−q​+Θ1​ϵt−s​+Θ2​ϵt−2s​+...+ΘQ​ϵt−Q∗s​+ϵt​这里博主举个栗子给大家加深理解（如果栗子有质量问题，欢迎大家在评论区中留言讨论），假如需要对某品牌羽绒服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

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

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最小原则方法定阶的结果如图所示，定阶模型为

S

A

R

I

M

A

(

1

,

1

,

2

)

(

0

,

1

,

2

)

12

SARIMA(1,1,2)(0,1,2)_{12}

SARIMA(1,1,2)(0,1,2)12​

5、模型的建立与检验

建立并训练

S

A

R

I

M

A

(

1

,

1

,

2

)

(

0

,

1

,

2

)

12

SARIMA(1,1,2)(0,1,2)_{12}

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