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

Bill.FF

已于 2023-10-13 20:48:15 修改

阅读量2.5w

点赞数 105

分类专栏：时序预测文章标签： python 时序数据库

于 2023-10-13 20:40:35 首次发布

本文链接：https://blog.csdn.net/qq_48521968/article/details/133742331

版权

时序预测专栏收录该内容

3 篇文章 13 订阅

订阅专栏

文章目录

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

该文主要讲解有关季节性时序模型——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等于1
D	季节差分次数（通常不会大于1）	一般就是 0 或 1，做了季节差分就是 1
S	季节/周期长度	该季节性序列的周期大小（注意不一定是12）

$SARIMA(p,d,q)(P,D,Q)_S$ 模型具体公式如下
$\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$ 是不是感觉很复杂（反正博主第一次学的时候是一头雾水），下面我们拆开讲讲公式中各个符号的含义与作用:

延迟算子 $B$
延迟算子类似于一个时间指针，当前序列值乘以一个延迟算子，就相当于把当前序列值的时间向过去拨了一个时刻 $x_{t-p}=B^px_t$ 特别地， $1-B)x_t=x_t-x_{t-1}$ 表示对 $x_t$ 做一阶差分，以此类推 $1-B)^dx_t$ 表示对 $x_t$ 做d阶差分。
季节性延迟算子 $B_s$
该延迟算子相当于把当前序列值的时间向过去回拨了一个周期 $x_{t-s*P}=B_s^Px_t$ 特别地， $1-B_s)x_t=x_t-x_{t-s}$ 可以通俗理解为对 $x_t$ 做了一次周期为s的季节性差分，以此类推 $1-B_s)^Dx_t$ 表示对 $x_t$ 做D次周期为s的季节性差分
。
延迟多项式算子 $\phi(B)$ 与 $\theta(B)$
$\phi(B)=1-\phi_1B-...-\phi_pB^p$ $\theta(B)=1-\theta_1B-...-\theta_qB^q$ 如AR模型可以简化为 $\phi(B)x_t=\epsilon_t$ ，同理ARMA模型可以简化为 $\phi(B)x_t=\theta(B)\epsilon_t$ ，ARIMA模型可以简化为 $\phi(B)(1-B)^dx_t=\theta(B)\epsilon_t$
季节性延迟多项式算子 $\Phi(B_s)$ 与 $\Theta(B_s)$
$\Phi(B_s)=1-\Phi_1B_s-...-\Phi_pB_s^P$ $\Theta(B_s)=1-\Theta_1B_s-...-\Theta_qB^Q_s$

所以，这个看似复杂的SARIMA模型其实就是对一个具有季节性特点的时间序列{ $x_t$ }做D次季节性差分（去周期）和d次常差分（去趋势）得到一个新的序列{ $y_t$ }，这段话用数学公式表示为： $y_t=(1-B)^d(1-B_s)^Dx_t$ 。然后对差分后的{ $y_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$ 这里博主举个栗子给大家加深理解（如果栗子有质量问题，欢迎大家在评论区中留言讨论），假如需要对某品牌羽绒服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