接上文笔记(一),建议先看上文,
上文我们讲述了时间序列基本知识,ARIMA等
本文主要是
(1)SARIMA
(2)自动、手动确定参数
(3)一些其他需要注意的点;
1.接上文,导数据,基本查看
import numpy as np
import pandas as pd
from statsmodels.datasets import co2
from statsmodels.tsa.stattools import adfuller
from matplotlib import pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.stats.stattools import durbin_watson #DW检验
plt.rcParams['font.sans-serif'] = ['SimHei'] # 显示中文
plt.rcParams['axes.unicode_minus'] = False # 显示负号
# 加载数据
data = co2.load().data
# 月度数据
data = data['co2']
data = data.resample('M').mean().ffill()
# 建画布
fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(223)
ax3 = fig.add_subplot(224)
# 可视化整体数据
data.plot(ax=ax1)
ax1.set_title('co2 data')
# 周期性规律
data[-36:].plot(ax=ax2)
ax2.set_title('最后36个数据')
data[:40].plot(ax=ax3)
ax3.set_title('前40个数据')
plt.tight_layout()
plt.show()
我们把数据稍微截取一段,很明显能知道以一年12个月为单位,周期性波动,并且整体往上涨一点的趋势。
上文,我们是先将data手动1阶12滞后项(也可称1阶12步)差分,变成data_diff,然后用data_diff去看图预估最佳参数,再用其他库里面的方法,以AIC、BIC为评价标准找最佳参数(p,q),找到后对两个组合看了一下模型报告summary和四个图(残差折线、直方图、QQ图、自相关图)。
模型报告summary里都不完全符合正态分布,其实是小事,哪能这么精准,主要是p,q两个组合,自相关图都有柱子冒尖,超过了蓝色范围,这是一个不可接受的问题,证明模型中信息还没有完全提取完;
因为单纯的ARIMA处理这种很明显的周期性数据,本身就不太行。
2、用SARIMA
其中的S--seasonal,直译是季节性,其实也包含周期性的意思,数据一个星期一个周期,则周期性为7,数据按月统计每年一个波动周期,则周期为12;
从ARIMA到SARIMA只需要一步,多设置一个seasonal_order参数即可;
ARIMA的order=(p,d,q)
(S)ARIMA的seasonal_order=(P,D,Q,s) 注意前面三个是大写的;
其中 s,季节长度,或者说是周期大小,可以提前定死是12,这是我们人脑判断的,也只能人脑先判断;
其中D,季节差分的次数(最大为1,一般就是0或1);
求中P,季节自回归的阶数(通常不会大于3)
其中Q,季节移动平均的阶数(通常不会大于3)
我们相当于,原本的数据,用ARIMA,把其中的规律,分割出来一部分,交由seasonal_order来处理,当我们用了seasonal_order后,小写的p,q,d基本也会随之变化,一般都会变小一些,毕竟规律要移交一部分给seasonal_order。
2.1 参数的预估
单就本数据集案例,我们之前ARIMA模型自动判断的p,q失效了;
我们目前面临7个参数,
s=12,板上钉钉;
D基本预计为1,原数据我们1阶12滞后差分即变成了非常好的平稳,这12步滞后,交给了s,这1阶我们想要交给大写的D;
d 基本上要么0要么1,因为
其中d+D最大值为2,超过了2绝对有问题!
P,Q不会很大,一般顶多到2或者3;
p,q完全不知道,但感觉不会超过ARIMA的最大值,一般会降低;
2.2 自动确定参数
目前确定这些参数的唯一办法,是比较无脑的网格搜索,因为时间序列的数据维度小,所以搜索挺快的,并且库中自带该方法,其中搜索算法也有相应优化,也不是单纯无脑试;
实现搜索有两个途径,1是自带方法,2是自己写个循环;
用库里自带方法
import pmdarima as pm # 没有这个就pip
model_auto_s = pm.auto_arima(y=data # 用原始数据
,start_p=1 # 网格搜索p起始值
,d=None# None可让程序自己定
,start_q=1 # 初始q值
,max_p=6# 之前搜索的最大值6,
,max_d=1 # 之前我们知道一阶差分即可
,max_q=18 # 前面看到q拖尾很严重
,m=12 # 确定周期性参数,12月一轮回
,seasonal=True # 启动季节性参数
,D=1 # 一次季节性差分即可
,max_P=3 # 一般不会超过3
,max_Q=3 # 一般不会超过3
,information_criterion='aic'# 只能选一个
,n_jobs=-1 # 开启全核心,但没卵用,默认是类似贝叶斯调参
,maxiter=200 # 最多尝试次数
,trace=True# 打印过程信息
,suppress_warnings=True# 默认值,可消除很多不重要警告
,stepwise=True # 默认值,类似贝叶斯调参,说是又快又好
,random=False # 默认值,不搞随机网格搜索,
)
print(model_auto_s.summary())
我们预估D为1,看看程序会怎么搜索最佳值;
注意我们用的是原数据,差分、滞后项,均交给模型去处理,方便很多;
其中stepwise默认为True,参数说明上写着开启则会智能选择最佳参数,有些贝叶斯调参的感觉,具体咱也不清楚。
其中n_jobs=-1只有当stepwise=False时,就变成暴力全部可能的组合搜索;
information_criterion选项,只能选择一个,一般aic、bic可以轮番看看;
其他参数均很容易看懂,亦有注释;
不到半分钟即可找到最佳参数:
以AIC为标准:
可以将标准换成BIC看看
进行12个月的未来预测:
模型检验(是看残差的情况):
①参数折线图,差不多是有个稳定均值
②非常正态
③大部分都在红线范围附近
④残差自相关系数终于没有冒尖的了
可以看到,带季节性的SARIMA感觉上就比较好;
画个图看看
plt.figure(figsize=(12,7))
# 只挑一部分画,全画了看不清细节
data[-36:].plot()
final_df['result'][-12:].plot(label='上文的结果')
y_3.plot(label='原数据SARIMA(1,1,1)(0,1,1)[12]')
plt.legend()
可以看到,上文的ARIMA激进一些,本文的要低一些;如果自己拿到数据,看看波峰,可知不是每个周期顶点都会上涨很多,这第一个波峰和第二个波峰,也就相差0.5不到的样子。
置信区间:
可在predict参数中,让模型返回预测的上、下区间
一般结合fill_between画阴影带:
plt.figure(figsize=(12,7))
data[-36:].plot()
final_df['result'][-12:].plot(label='上文的结果')
y_3.plot(label='原数据SARIMA(1,1,1)(0,1,1)[12]')
plt.fill_between(y_3.index,
bb[:,0],
bb[:,1],color='k',alpha=.2)
plt.legend()
可知,我们不用手动差分滞后,不用还原,对比上文的方法,非常滴方便。
当然,你可以先预留出原数据的最后一年,通过MSE来看ARIMA和SARIMA哪个更好,但是上文文末所说的问题,的确存在且难以完全解决;
但是我们上文的model_easy(1,0,1)和model_hard(6,0,16)的AIC和BIC,都远大于本文中的320,330,本文的SARIMA多项指标均明显不低于上文,我们有理由认为SARIMA更好。
所以,以后碰到明显周期性的数据,可以先用ARIMA试试水,但是要清楚SARIMA很可能更好;
如果嫌麻烦,其实可以通过acf,pacf预估出小p,q ,再直奔无脑搜索,反正周期性人眼定死了,D大概为1,d+D最大为2,大PQ又不会很大,找出最佳的组合,检验一下模型即可。
当然一篇文章不可能写全,可以看看这些参考!
2.3 人工预估+手动写循环找最佳
这又是另外一个思路了:
基本是参考这篇文章:时序预测 | SARIMA(P、D、Q了解一下) - 知乎 (zhihu.com)
原理是,通过原数据(不进行差分,且原数据一般都不平稳),画出ACF,PACF图,像本文数据的周期为12,则P通过PACF的12,24,36阶偏自相关系数,根据蓝色范围来预估,Q则同理看ACF的12,24,36阶,因为P,Q基本不超过3,超过则极其容易过拟合;
同时p,q通过先训练一个ARIMA,搜索出p,q最大值,那么SARIMA的p,q一般不会高于ARIMA里的p,q。
最后再写循环判断P,Q,p,q,d
如果对自动循环判断的觉得不太行,可以手动小调一下p,q的值,如果对时间序列略有经验,人工手动调优目前可能略优于机器搜索。
还有这一篇,也是先心理预估几个参数的最大值,再手动写循环,但是其预估的依据有所不同,这篇是通过差分后的数据,去预估几个参数的最大值。
python 时间序列预测 —— SARIMA_python sarima-CSDN博客
这番操作把我也给整迷糊了,到底是看差分前还是差分后的季节性倍数阶来预估大P,Q,感觉貌似都有点道理,我个人认为让机器去选最省心。
其他参考:
ARIMA Model for Time Series Forecasting (kaggle.com)
(上面有SARIMAX,多一个X实际上就是多了一列数据,就本文案例来讲,就是将月份单独提取出来做一个特征列,也喂给模型,这样模型会更容易知道哪个月份数据是波峰、波谷,可能有更好的效果)
A Guide to Time Series Forecasting with ARIMA in Python 3 | DigitalOcean
时间序列模型(五):时间序列案例_实现销售额预测 - 知乎 (zhihu.com)
(写的不错,建议阅读)