时间序列模型笔记（二）

TravelLight92

已于 2024-05-29 13:50:13 修改

阅读量843

点赞数 7

文章标签：笔记

于 2024-05-15 14:47:05 首次发布

本文链接：https://blog.csdn.net/TravelLight92/article/details/138905241

版权

接上文笔记（一），建议先看上文，

时间序列模型笔记（一）-CSDN博客

上文我们讲述了时间序列基本知识，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()