时间序列模型笔记（一）

TravelLight92

已于 2024-05-29 13:51:34 修改

阅读量1.4k

点赞数 23

文章标签：学习笔记

于 2024-05-12 19:42:18 首次发布

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

版权

总结一下学习的时间序列相关知识：

一方面把所需要掌握的知识点都说到，同时对一些不太会用到的不作过多讲述，但将深入点的知识链接附上;

二是将完整的流程进行展示，带上数据集案例;

三是作为知识分享，不写屎山，让代码简单直观，因为自己看别人文章，一些简单的功能，非要写得花里胡哨，让人看得脑壳疼。

主要的关注点在ARIMA(SARIMA)：

主要内容：1.重要基础知识点 2.ARIMA及其衍生 3.定p,q,d方法和要点 4.其他应该知道的

参考内容比较多，链接穿插其中：

时间序列，内容其实还真不少，并非一两篇文章能够讲清楚，我也是看了不下二三十篇文章，每个人的路子都有些差异，本人综合一下，帮看官捋清思路。

一、重要基础知识

时间序列，仅含有 连续的时间变量、要预测的变量两列，虽然从现在来看,预测通常不会特别准，但仅依靠一个变量去预测，能取得这样的准确率、并在信息有限的情况下，能实现预测，已经不错了，同时，作为学机器学习、算法的非金融行业者来说，时间序列是要懂的，但没必要掌握得特别深。

同时要知道：时间序列模型本身不太难，只不过网上能简单查到的，很多人实现过程风格迥异，本文经历痛苦地观察各家所长，筛选出比较正规、靠谱的方式，也会提及一些其他库、函数（方法）的实现。

1.1平稳

定义：一组时间序列数据各阶统计特征（如均值、方差、协方差...）不随时间的变化而变化。其数学定义又分为严平稳和宽平稳。

让时间序列满足平稳非常重要，因为只有基于平稳过程的模型才能用来进行预测和估计参数。

严平稳/宽平稳过程-CSDN博客

严平稳，又称为强平稳，是指任意时刻的联合分布都与时间无关，太过理想化，一般到手的数据不会符合这么严格的要求：

宽平稳，又称为弱平稳，如果某个时间序列，均值、方差基本是个确定的值，与时间无关，抽取两个等长的子集，其对应的变量Y值的均值、方差是相同的，协方差只与时间差h有关，则该时间序列是宽平稳的。

不同文章定义可能描述不同，意思差不多；

严平稳比宽平稳的要求更严格，但严平稳并非完全包含宽平稳。

为什么要求时间序列的平稳性？因为非平稳时间序列太过杂乱无章，有的甚至完全无规律可循，而平稳时间序列本身存在某种分布规律，前后具有一定自相关性，进而可以利用这些信息帮助预测未来。

就比如上图数据，原数据周期性上涨，我们可以将原数据，分解成三个部分，逐步上涨的趋势、周期性波动、噪音（难以完全掌握规律的部分），趋势是比较简单的，预测未来的数据，按照趋势往上走一点即可，再叠加上周期性变动+残差（能找多少规律就找多少放进去），即可完成预测。

固其中的周期性变动，得是均值、方差稳定的，得满足弱平稳条件，残差部分，也不能一会很大一会很小太离谱。

1.2白噪声

具有零均值同方差的独立同分布序列，当均值为0，为高斯白噪声或正态白噪声。特点是，虽然平稳，但下一时刻的值跟前面任意时刻没有任何关系，作用是：

①假如要研究的时间序列为白噪声，则数据根本没什么可建模研究的规律，直接放弃，硬要预测则直接预测未来全部是均值，则总体MSE（residual）最小，没有意义。

②假如将原始时间序列的信息，用模型进行提取学习，剩下的内容为白噪声，则认为模型提取了所有可用信息，没有遗漏，则模型起码没有错；

③白噪声是经典的，虽然平稳但是不能预测的数据，在时间序列里面各项检验通常提及它。

简单的白噪声示例：

import numpy as np
from matplotlib import pyplot as plt

white_noise = np.random.standard_normal(size=1000)
plt.figure(figsize=(12, 6))
plt.plot(white_noise)
plt.show()

1.3非平稳时间序列

大多数时间序列都是非平稳的，一般可以通过差分、取对数等方法转化成平稳时间序列，时间序列里常用差分。

1.4随机游走

很经典，是累加和，虽然下一步与上一步息息相关，但无法用时间序列相关模型进行预测，可通过差分后变成白噪声，而白噪声本身也无法进行预测未来，故随机游走也是无解题，其方差逐步增大，不平稳；

概念&术语 | 谈谈时间序列的平稳性

import numpy as np
from matplotlib import pyplot as plt

y1 = np.random.standard_normal(size=1000)
y2 = np.random.standard_normal(size=1000)
y1 = np.cumsum(y1)
y2 = np.cumsum(y2)
plt.figure(figsize=(12, 6))
plt.plot(y1)
plt.plot(y2)
plt.show()

1.5AR模型

时间序列模型(二)：AR模型 - 知乎 (zhihu.com)

AR模型，即自回归模型（Autoregressive Model），一个自回归模型的形式是描述某个变量与其过去值的关系的。对于一个时间序列数据，自回归模型的一阶形式可以表示为：

$Y_{t+1}=c+\phi Y_{t}+\varepsilon_t$

其中c为常数，\phi 为回归系数，Y_t为上一个时刻的值，\varepsilon（类似e的数学符号）是误差项。

当阶数为p个的时候：

$Y_{t+1}=c+\phi_1 Y_{t}+\phi_2 Y_{t-1}+...+\phi_p Y_{t+1-p}+\varepsilon_t$

即下一时刻的值，与之前P个值相关。

1.6MA模型

时间序列模型(三)：MA模型 - 知乎 (zhihu.com)

MA移动平均模型（Moving Average),MA模型的基本思想是：大部分时候时间序列应当是相对稳定的。在稳定的基础上，每个时间点上的标签值受过去一段时间内、不可预料的各种偶然事件影响而波动，MA部分就是寻找残差\噪音的规律。

AR（自回归）模型和MA（移动平均）模型虽然在形式上看起来很相似，但是他们的关键区别在于他们对过去信息的处理方式。

对于自回归（AR）模型来说：在这种模型中，当前值是过去值的函数。也就是说，我们是在使用过去的"实际"观察值来预测现在的值。AR模型的基本思想是过去的观测值会对未来的观测值产生影响，即未来的观测值是过去观测值的加权和。

对MA模型来说：在这种模型中，当前值是过去噪声（或称之为误差或冲击）的函数。这里的“白噪声”实际上是模型无法解释的随机部分，是未能被模型捕获的信息。换句话说，MA模型是在试图用过去的"错误"或"冲击"来预测现在的值。

1.7差分与滞后

常用的ARIMA模型由AR+I+MA，其中的I即为差分，可消除数据的波动，使时间序列的值变平稳，符合建模条件。

差分的阶数：一般为0,1,2次，每次差分，原始数据的信息就相当于被稀释，3次差分基本上很少见，1阶差分能使数据平稳，就别用2阶，越少越好。

差分的滞后：在1阶差分中，用后面的数据减前面的数据，则滞后为1，用下周一的数据减上周一的数据，则滞后为7，滞后一般观察数据的周期性，比如时间序列按天统计数值，每周规律类似，则滞后为7，按月统计，周期性一般为年，滞后用12，经常使用几个经验值进行判断。

当然经验判断很可能不准，模型亦有方法判断出最合理的滞后项，当模型判断的最佳==经验、直观主观判断的估计值，则基本稳妥了。

详见时间序列模型(四)：ARIMA模型 - 知乎 (zhihu.com)

1.8ARIMA

经常使用的时间序列模型，在 ARIMA(p, d, q) 模型中：

①p 代表 "自回归部分 (Autoregressive)"：模型中使用观测值的滞后值（即前面 p 个期的值）。自回归模型的出发点是认为观测值是它前面的 p 个值的线性组合，数学形式如下：

②q 代表 "移动平均部分 (Moving Average)"：模型中使用的错误项的滞后值（即前面 q 个期的值）。移动平均模型是将当前值和过去的白噪声之间建立关系，数学形式如下：

③d就是差分的阶数。

1.9ACF自相关系数

自相关函数ACF (Auto-Correlation Function)，

图片摘自上文链接，其公式表达可能不是很直观，首先相关系数，要求两个序列，数量相同，即等长；

$\rho_k=Cov([x_1,x_2,x_3...x_{t-k}],[x_{k-1},x_{k},x_{k+1}...x_t])/\sigma_{1}\sigma_{2})$

其中两段数据等长，simga1，2为两个数据集的标准差，即皮尔逊相关系数，由统计学书本证明过，相关系数的范围一定是【-1,1】之间。

想要实现手算，可参考文章:ACF与PACF

1.10 PACF偏自相关系数

自相关系数比如时刻t的值，衡量其受到t-1,t-2,....时刻的影响，每个数都会受到前面的数的影响，例如想要知道t的值受到t-2的影响:

t-2首先影响t-1,

t-2也会影响t,

同时t-1也会影响t,

此时我们想去掉中间商，了解t受到t-2的单独的影响，就需要用到偏自相关系数（PACF）来衡量。

PACF计算相对复杂一些，可查看偏自相关系数PACF（公式篇） - 知乎 (zhihu.com)

时间序列 | PACF篇

这篇文章从最小二乘法（简单情况下的求解）、尤尔-沃克方程（Yule-Walker equation）和伯格算法（Burg"s method）三种方式给出了求解PACF的方法，实际上可以不需要知道如何求出来的，因为计算方法不是那么简单....

根据作者推荐，第三种burg方法算的结果一般准一些。

那么ACF和PACF用来干什么呢？

作用：

1.判断序列是否平稳，平稳序列，往往带一些自相关性，但不会特别多，即自相关系数一般是个位数，但后面再往上到十位数的时候，会迅速退化到零（滞后期越短相关性越高，滞后期为0时，相关性为1）；反过来，非平稳的数据，退化会发生得更慢，或存在先减后增或者周期性的波动等变动。

2.判断序列是否白噪声，白噪声完全无自相关性，除0阶自相关系数为1外，其他阶基本为0，或在蓝色范围内。

3.判断是否过度差分，1阶的自相关系数一般为正数，如果它为负，此时要格外注意，你可能过度差分了，或者数据不符合时间序列平稳、非白噪声等要求。

4.确定AR MA 中的q,p，即定阶。

可能模型适合用AR，可能适合MA，可能适合ARIMA，一般都要先画图看看，只是先让人肉眼看看心中大致有个数，最终都是用模型自动判断最佳值，并验证它的确是最佳值。

AR中的p往往由PACF确定，MA中的q往往由ACF确定。

①ACF拖尾，PACF在p步后截尾，则用AR（p）模型

②ACF在q步截尾，PACF呈拖尾，则用MA(q)模型

一般情况下，都不会单纯地只用AR、MA，都是用ARIMA，当其中一个为0时，则相当于去掉了另外一个模型的作用。

③ACF,PACF均拖尾，什么意思？两个图都有几步在蓝色范围外，后面数值基本都在蓝色范围内，我们认为其边界即为p,q，认为它是在几步后截尾（实际并没有截尾，我们就只取蓝色范围外的），则用ARMA（p,q）模型，其p,q按照截尾数定，一般要尝试多个组合，按照评估标准选择最佳。

④复杂点的情况，ACF,PACF均让人无法直观判断，statsmodels库里有方法可以自动定最佳的p,q。

最好的情况是，库里面自动判断的方法==看图给出的结果，两个都说他是最佳，则基本稳妥了，如果图和库里面自动判断的不同，则不要相信眼睛，相信数学给出的判断，同时要对比分析下一些metric评估标准，确定下数学真的判断对了。

1.11 拖尾与截尾

截尾：相关系数在某一步之后为0
拖尾：相关系数不能在某一步之后为0（截尾），而是按指数衰减（或成正弦波形式)

如何判断拖尾和截尾？
1. 如果样本自相关系数（或偏自相关系数）在最初的d阶明显大于2倍标准差范围（在蓝色范围外），而后几乎95%的样本自相关（偏自相关）系数都落在2倍标准差范围以内（蓝圈内），而且其衰减为小值的过程非常突然，这时，通常视为自相关（偏自相关）系数截尾。
2. 如果有超过5%的样本相关系数落在2倍标准差范围以外，或者是由显著非零的相关函数衰减为小值波动的过程比较缓慢或者非常连续，这时，通常视为相关系数不截尾。

平稳的序列的自相关图和偏相关图不是拖尾就是截尾。

1.12 趋势性、季节性等规律

参考：

ACF与PACF

①趋势性，如下图，从上到下依次为原数据、ACF、PACF，可知数据在不停上涨，ACF拖尾且为正，随着滞后阶数的增加而缓慢下降。（乌托邦式理想数据，实际少见）

②季节性

当数据具有季节性时，自相关值在滞后阶数与季节周期相同时（或者在季节周期的倍数）较大。

③组合效应

同时具有趋势和季节性，又周期波动、又不停上涨，如下图所示：

接下来补充几个纯数学的时间序列模型

1.13 简单移动平均（Simple Moving Average）

接下来补充几个简单点的模型：

参考：简单易用 & 炒股必看的时序预测基本方法--移动平均（SMA、EMA、WMA）

较为简单的预测工具，未来时刻的预测值，为之前n个数值的平均，n可以自由选择

$x_{t+1} = \frac{1}{n} \sum_{t-n+1}^{t} x_i$

有个rolling方法可以比较容易实现，因为太无脑+简单，一般不会用，了解即可。

# 利用skshare的某酒厂股价数据（入市有风险，容易血本无归）

# 获取数据
df = ak.stock_zh_a_hist(symbol="600519", start_date="20180308", end_date='20210308')
df = df.set_index('日期')
df.index = pd.to_datetime(df.index)
df = df[['收盘']]
# 进行简单移动平均
df['SMA_5'] = df['收盘'].rolling(window=5).mean()
df['SMA_10'] = df['收盘'].rolling(window=10).mean()
df

1.14加权移动平均（Weighted Moving Average）

与SMA类似，但是在计算时各历史数值权重不同，一般近的权重大、远的权重小，例如最近一个权重为n,最近第二个为n-1,直至到权重为1的；因为权重过大，需要除以n+(n-1)+(n-2)...+1来归一化

$x_{t+1} = \frac{nx_t+(n-1)x_{t-1}+...2x_{t-n+2}+1x_{t-n+1}}{n+(n-1)+...2+1}$

随便写个函数即可,配合rolling窗口函数，会方便不少，代码少几行：

当然，这个方法也因为太简单无脑，基本不会用。

def weighted_ma(data,n):
    weight = np.arange(1,n+1)
    total_weight = weight.sum()
    
    new = data.rolling(window=n).apply(lambda x:(weight*x).sum()/total_weight,raw=False)
    return new

# 参考的作者也分享了一个有意思的计算方式
# 有一种梯形计算面积的感觉
def interesting_wma(data,n):
    return data.rolling(n).apply(lambda x:x[::-1].cumsum().sum()*2/n/(n+1))

df['WMA_5'] = weighted_ma(df['收盘'],5)
df['WMA_5_cumsum'] = interesting_wma(df['收盘'],5)
df

1.15指数移动平均（Exponential Moving Average）

EMA亦称(Exponential Weighted Moving Average)，即EWMA，和WMA类似，但每个值的权重为系数的乘法

$x_{t+1} = \frac{x_t+(1-\alpha)x_{t-1}+(1-\alpha)^2x_{t-2}+...(1-\alpha)^tx_0}{1+(1-\alpha)+(1-\alpha)^2+...(1-\alpha)^t}$

\alpha明显是介于0和1之间，越近的值所占权重明显更大，越远的权重呈指数减少，上述公式中用了前t+1个值。

这个方法，也强不到哪里去，数据不会是这么简单的指数相乘就能拟合的。

没想到pandas功能这么强大，连EWMA的计算方式都自带：

关于ewm的函数具体参数见下链接：

python pandas ewm 一次指数加权移动平均-CSDN博客

------------------------------------------

上面介绍了三种简单的方法，由于太简单，感觉使用场景十分有限，当个基础知识了解即可，不必深究，也没什么深究的点

------------------------------------------

1.16 指数平滑

包含一次指数平滑、二次指数平滑、三次指数平滑（Holt-Winters）等知识内容，内容较多，包含了对趋势、季节性因素、阻尼趋势、加法OR乘法等内容，个人感觉这有点偏金融行业了，可以了解一下。

参考：[Formula&Excel&Python] 一次指数平滑、二次指数平滑、三次指数平滑（Holt-Winters）

接下来，我们按照ARIMA模型顺序，依步骤讲解：

二、检验时间序列的平稳性

时间序列建模：首先要求序列是平稳的，该内容按照写这个文章逻辑划分，应该放在基础知识里面的，但是内容太多，单独作为一个一级标题，可参考：

时间序列的平稳性检验方法汇总

时间序列平稳性检验方法，可分为三类：

图形分析方法
简单统计方法（可忽略）
假设检验方法

2.1图形分析方法

拿到数据，肯定要画个折线图看看，同时把（ACF+PACF）两张图一并画上。

2.1.1可视化原数据

最基本、最简单直接的方法，即绘制图形，肉眼判断，绘制折线图，看曲线是否围绕某一数值上下波动（判断均值是否稳定），看曲线上下波动幅度变化大不大（判断方差是否稳定），看曲线不同时间段波动的频率[~紧凑程度]变化大不大（判断协方差是否稳定），以此来判断时间序列是否是平稳的。

akshare是一个金融方面，获取数据很好的外部库

你不能仅凭肉眼判断一个序列是平稳的，但一个序列不平稳，你可以一眼看出来！

import akshare as ak
# -------------- 准备数据 --------------
# 白噪声
white_noise = np.random.standard_normal(size=1000)

# 随机游走
x = np.random.standard_normal(size=1000)
random_walk = np.cumsum(x)

# GDP
df = ak.macro_china_gdp()
df = df.set_index('季度')
new_index = [x[:5] for x in df.index] # 索引太长改短点
df.index = new_index
gdp = df['国内生产总值-绝对值'][::-1].astype('float')

# GDP DIFF
gdp_diff = gdp.diff(4).dropna()
# -------------- 绘制图形 --------------
fig, ax = plt.subplots(2, 2,figsize=(10,8))

ax[0][0].plot(white_noise)
ax[0][0].set_title('white_noise')
ax[0][1].plot(random_walk)
ax[0][1].set_title('random_walk')

ax[1][0].plot(gdp)
ax[1][0].set_title('gdp')
ax[1][0].set_xticklabels([i[:-1] for i in gdp.index],rotation=270)

ax[1][1].plot(gdp_diff)
ax[1][1].set_title('gdp_diff')
ax[1][1].set_xticklabels([i for i in gdp_diff.index],rotation=270)
plt.tight_layout()
plt.show()

2.1.2可视化统计特征

指绘制时间序列的自相关图和偏自相关图，根据自相关图的表现来判断序列是否平稳。

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

fig, ax = plt.subplots(4, 2,figsize=(10,8))
fig.subplots_adjust(hspace=0.5)

plot_acf(white_noise, ax=ax[0][0])
ax[0][0].set_title('ACF(white_noise)')
plot_pacf(white_noise, ax=ax[0][1])
ax[0][1].set_title('PACF(white_noise)')

plot_acf(random_walk, ax=ax[1][0])
ax[1][0].set_title('ACF(random_walk)')
plot_pacf(random_walk, ax=ax[1][1])
ax[1][1].set_title('PACF(random_walk)')

plot_acf(gdp, ax=ax[2][0])
ax[2][0].set_title('ACF(gdp)')
plot_pacf(gdp, ax=ax[2][1])
ax[2][1].set_title('PACF(gdp)')

plot_acf(gdp_diff, ax=ax[3][0])
ax[3][0].set_title('ACF(gdp_diff)')
plot_pacf(gdp_diff, ax=ax[3][1])
ax[3][1].set_title('PACF(gdp_diff)')

plt.tight_layout()
plt.show()

①白噪声，虽然是平稳序列，但后面的数据跟前面的没有关系，故白噪声无建模价值；

②随机游走，虽然自相关系数拖尾，由于均值不停变化、方差不停增大，非平稳；

③GDP数据，原数据图明显看到均值在往上涨的趋势，且有周期性，不平稳，ACF和PACF图都不用看了；

④GDP差分后的数据，看原数据图，可以勉强认为是平稳的，算通过肉眼考验，ACF和PACF图不太明朗，就不用肉眼看了，用模型自带的方法判断。

肉眼看，带有较强主观性，判断不平稳比较准，判断平稳，则不是那么敢确定。

就像判断一个人漂亮否，各有所爱，但判断一个人丑不丑，大家都能一致地确定他/她很丑。

2.2简单统计方法（不重要）

只是作为一个补充，了解即可，简单将数据，切成前后两半，看均值、方差是否差异很大。

2.3假设检验方法（重要）

假设 $x_{t+1}=\phi x_t+\varepsilon_t$ ,未来的数据是之前数据乘以一个系数+一个随机影响数.

如果 $\phi$ 的绝对值>1,则该序列爆炸式增长，且非平稳；

如果 $\phi$ 的绝对值=1,则该序列无规律，非平稳；

如果 $\phi$ 的绝对值<1,则该序列随时间t变化最终会收敛，长期来看是平稳的；

如果 $\phi$ =1，即是存在单位根，我们希望没有单位根，不然没法进行时间序列建模。

单位根（unit root）检验就是检验该差分方程的特征方程（characteristic equation）的各个特征根（characteristic root）是均小于1，还是存在等于1的情况。没有检验均大于1的情况，是因为当根均大于1时为爆炸型发散序列，日常数据中基本不存在。

2.3.1 DF检验

即ADF检验的前一版本，DF检验只能用于一阶自回归过程，即下一时刻只受到上一时刻的影响，不能考虑更早之前的影响，较少用。

2.3.2 ADF检验

能适用于高阶自回归过程的平稳性检验

ADF假设条件同DF：

原假设H0：存在单位根( $\phi=1$ )，时间序列不平稳

备择假设H1：不存在单位根( $\phi <1$ ),时间序列是平稳的

我们希望要研究的序列，是平稳的，或者通过差分等方差，变成平稳的，希望能拒绝原假设。

一般可用两个外部库，还是用上面的GDP数据：

① statsmodels库中的adfuller

先不差分，第一个值为检验统计量，第二个为p值，第三个为滞后数，第四个为样本数，括号中为不同的置信水平下检验统计量临界值。

简单版：先看第二个p值，原假设数据是不平稳的，所以我们希望p<0.05，从而拒绝原假设，此处p快到1了，证明原数据肯定不平稳；其次看第一个检验统计量Z值，都知道正态分布中，均值 $\mu$ - $2\sigma$ , $\mu-3\sigma$ 的范围，包含了95%和99.7%的数据，其对应了x轴的一个点位，如果发生了极小概率时间，那么Z值对应的X轴上的点肯定会在非常靠左或者靠右的位置（这里是单尾检验，一般我们希望Z值在1%的更左边才能拒绝原假设），通过对比第一个值和括号内三个点位的值，看Z值处于范围内还是范围外，此处0.79处于接受域内，证明在原假设为真的情况下，发生了一个很经常出现、概率很大的事，就是原数据极大概率是不平稳的。

上面说了看两个东西，一般p值很大，就不需要看检验统计量Z值，p值很小，一般还要看看Z值，如同生活中一件事明显错误，一眼可看出错得离谱，而想要确定正确的事，得多看两眼；

各项数值含义，可参考本文：时间序列数据的ADF检验—adfuller()函数的模型拟合系数介绍-CSDN博客

from statsmodels.tsa.stattools import adfuller
result = adfuller(gdp,regression='c',maxlag=4)
result

② arch库中的ADF

from arch.unitroot import ADF
adf = ADF(gdp,lags=4,trend='c')
print(adf.summary())

一般地，statsmodels库中的adfuller较为常用，其中参数，看看英文解释即可秒懂！

接下来，我们的演示数据中，gdp进行了1阶4滞后的差分，来看看差分后是否平稳。

两个库中方法，都是一样的，只要参数都设成相同的，结果必然相同。

可以看到，差分后平稳了。

2.3.3 PP检验、KPSS检验等其他方法：

参考

时间序列的平稳性检验方法汇总

这篇文章介绍了很多方法，该系列作者文章水平较高，写得很深入，好处是了解比较全面深入，缺点是很多公式要仔细了解原理，会有一定难度和门槛。

三、建模完整流程

偷一个图过来,整体流程如下，建议看完全文，再回头看看这个图，本文也是按正规流程来的，接下来每一步用代码+数据+讲解演示：

是时候换一拨数据集了，用statsmodels中自带的美国夏威夷莫纳罗亚天文台大气二氧化碳数据，并转换为月度数据。

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检验
# 加载数据
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')
# 相关图
plot_acf(data,ax=ax2)
ax2.set_title('自相关系数')
plot_pacf(data, ax=ax3)
ax3.set_title('偏自相关系数')
plt.tight_layout()
plt.show()

3.1 看看数据

原数据明显的周期性波动且趋势上涨！明显均值肯定不平稳！

一般都一次性把三个图都画了，省事！

3.1.1 分解

其实没必要做这一步，只是觉得这个方法很重要，写上去。

from statsmodels.tsa.seasonal import seasonal_decompose
decomposition = seasonal_decompose(data)
trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid
# 画图
plt.figure(figsize=(16, 12))
plt.subplot(411)
plt.plot(data)
ax = plt.gca()
ax.set_title('原始数据',fontsize=18)

plt.subplot(412)
plt.plot(trend)
ax = plt.gca()
ax.set_title('趋势',fontsize=18)

plt.subplot(413)
plt.plot(seasonal)
ax = plt.gca()
ax.set_title('周期性变动',fontsize=18)

plt.subplot(414)
plt.plot(residual)
ax = plt.gca()
ax.set_title('残差\噪音',fontsize=18)

plt.tight_layout()
plt.show()

3.2 平稳性检验

3.2.1 原数据是否平稳

数据稳步增长，肯定不平稳，但是我们还是检验一下，同时进行1阶差分，差分的滞后项，我们希望是12，因为是月度数据。

result = adfuller(data)
result

如果你觉得这些数据看着不明所以，那么就这样！

result = adfuller(data)
result_beauty = pd.Series(result[0:4], index=['Test Statistic','p-value','Lags Used','Number of Observations Used'])
for k, v in result[4].items():
    result_beauty['Critical Value (%s)' % k] = v
result_beauty['maximized information criterion if autolag is not None.'] = result[5]
print(result_beauty)

当然，还有这个库！

# 另外一个库的实现方法
from arch.unitroot import ADF
adf = ADF(data)
print(adf)

3.2.2 变平稳

如果时间序列不平稳，使用一些方法使其变平稳；

一般有差分、取对数变化、分解（分解成长期趋势、季节变动、不规则波动）

常用差分，差分几阶，滞后为多少？怎么确定呢？

一般差分多为1阶、2阶

滞后项，如果原数据画图，有很明显的规律，我们希望滞后项就是那个规律，同时要进行adf检验，检验结果跟我们预想的一样，则稳妥了。

预期结果---1阶12滞后差分：

# 原数据是明显的以一年为单位周期性变化，选12滞后差分
# 季节差分
data_diff = data.diff(12).dropna()

data_diff.plot(figsize=(12,4))
plt.title('co2 - 1阶12滞后差分后数据')
plt.show()
# ADF检验
res = adfuller(data_diff)
res
# 差分后序列平稳。

循环检验最佳滞后项：

# 一个简单的循环
from statsmodels.tsa.stattools import adfuller
result_df = pd.DataFrame() # 存储结果
row_index=0 # 行索引
for i in range(1,25):
    data_diff = data.diff(i).dropna()
    result = adfuller(data_diff,regression='ctt')
    if result[1]<0.05: # p值小于0.05,就保存之
        result_df.loc[row_index,'滞后项']=i
        result_df.loc[row_index,'ADF统计量']=result[0]
        result_df.loc[row_index,'P值']=result[1]      
        result_df.loc[row_index,'标准差']=data_diff.std(ddof=0)
        result_df.loc[row_index,'方差']=data_diff.var(ddof=0)
        row_index+=1
        
result_df

虽然滞后项20从ADF统计量、P值上最佳，但12滞后项在经验预期、标准差上获胜，并且12<20,一般能小则不取大，所以还是选12，反正数值上都是不错的。

需要注意的是参数regression:其中默认c，表示数据基本是一个稳定的值，ct则数据是稳定并有上升或下降趋势，ctt则是ct基础上还有周期性变化，n则表示数据没有常数项和趋势，如果不选择正确的参数，结果可能会大大不同。

有的文章可能直接看原数据折线图，将滞后项定为12，不会进行循环检验，这个循环检验，是我个人添加的认为比较稳妥的一种方式，经验判断（心理预期）它不一定对，只有数据才有说服力，反正代码就那么几行也不费事。

3.3 检查白噪声

如果变平稳了，下一步要判断是否为白噪声，因为白噪声序列没有可预测规律，不能用于研究，

非白噪声，才能进行下一步，故我们要进行白噪声检验：

有以下3个方法，实际上我们必然要经过画图肉眼看+数据检验两个流程：

①画相关图：

差分后，我们定然还要再画一波图。

由定义知，白噪声完全无自相关性，除0阶自相关系数为1外，理想情况下，延迟k阶的样本自相关系数均为0。实际上由于样本序列的有限性，延迟k阶自相关系数并不完全为0，只要在0值附近即认为无自相关性。

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
fig = plt.figure(figsize=(20, 10))
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(223)
ax3 = fig.add_subplot(224)

ax1.plot(data_diff)
ax1.set_title('差分后的数据')
plot_acf(data_diff, ax=ax2)
plot_pacf(data_diff, ax=ax3)
plt.show()

自相关图和偏自相关图中浅蓝色条形的边界为 $\pm2 / \sqrt T$ ，T为序列的长度。由于随机扰动的存在，自相关系数并不严格等于0，如果一个序列中有较多自相关系数的值在边界之外，那么该序列很可能不是白噪声序列。

很多在蓝色范围以外，不像白噪声的样子。

②Box-Pierce检验

在Ljung-Box中指定boxpierce=True即可

详见时间序列预处理 | 白噪声检验

③Ljung-Box检验

是对Box-Pierce的更新增强版，看名字可知都有Box这个人的参与，意在解决小样本下精确度不高的缺陷：

原假设：序列滞后m项的值之间相互独立，即相关性系数为0，即序列为白噪声

备择假设：m项的值之间有相关性，即系数非0，非白噪声

我们对GDP差分数据试一下：

from statsmodels.stats.diagnostic import acorr_ljungbox
res = acorr_ljungbox(data_diff
                     , lags=24
                     , return_df=True,boxpierce=True,auto_lag=False)
print(res)

重点看两个卡方检验的p值，如果p值都很大，不能拒绝原假设（系数为0，很可能是白噪声），而差分后的数据，p值很小，系数不为0，该方法比画相关图看蓝色范围靠谱许多，这样我们肉眼看+数据检验答案相同，这步又稳妥了。

④还有一种实现Ljung-Box检验的方式为，调用statsmodels包中的acf函数，随便用哪个都行。

import statsmodels as sm
import pandas as pd
# 把白噪声拉出来对比
r, q, p = sm.tsa.stattools.acf(white_noise, nlags=12, qstat=True)
df = pd.DataFrame(zip(range(1,41), q, p), columns=['lag', 'lb_stat', 'lb_pvalue'])
print(df)
# 差分后的数据
print('*'*49)
r, q, p = sm.tsa.stattools.acf(data_diff, nlags=12, qstat=True)
df = pd.DataFrame(zip(range(1,41), q, p), columns=['lag', 'lb_stat', 'lb_pvalue'])
print(df)

可以看到，白噪声的相关系数，存在极大可能性为0，而我们差分后的数据，统计学上为0的可能性极小，这证明了其中包含规律，可用于建模研究。

四、定阶

看到这里，会发现，我们的ARIMA模型，其中的I部分，实际上我们已经手动进行了1阶12滞后差分，剩下的就是确定AR和MA部分，由于ARMA模型，很多21、22年的文章还能用，但在我写这篇文章的时候，已经被更新抛弃了，但是机智的我想到了用ARIMA模型，将其中的I部分设为0，来实现用ARMA模型的目的；

但是单纯地直接用ARIMA模型又有个问题，我们的I部分，实际上有两个参数，差分阶数1和滞后项12，而ARIMA模型只能传阶数这一个参数，所以我们只能将手动进行了1阶12滞后的数据，用ARMA模型训练。

此时你肯定会想，这12滞后项无法传入模型当参数，要ARIMA有何用，其实要想解决这个问题，要使用SARIMA模型，就是在ARIMA模型中再传入一个参数，即带了S这个季节性波动（季节性即稳定的周期性波动），可以将我们想要的滞后项传入当参数，但是我们必须先掌握ARIMA，才能迈向SARIMA。

回到正题，用ARMA建模时，确定p、q的问题称为定阶，一般有以下方法：

4.1画ACF和PACF图

在我们差分后，画了那三张图，这里再贴出来，在看下面两张图的时候，我们就应该对p,q有个大致预期,ACF图确定q,PACF图确定p。

先看ACF，第一根柱子是自相关，自己跟自己绝对是1，从第二根开始为1阶自相关系数，可以看到蓝色范围外的还不少，因为咱们这原数据太整齐了，甚至闭着眼按照之前的趋势，往上画一下也差不太远...可以看到，第16阶刚好卡在蓝色范围之上，第17阶差不多正好，第18阶就在范围内，我们预计q是16左右。

再看PACF，就稍微简单点，3阶基本卡在边界了，4阶又冒了点头，7阶基本很接近0，所以我们考虑p为2或4或6。

实际中数据可能更加复杂，反复横跳让人懵逼，所以先有个心理预期，提前知道大概的答案，

再交给模型自动判断。

4.2 模型自动判定

时间序列模型这点好，自带网格搜索方法自动判断最佳的p和q,如果和我们观察图像上的一致，则稳妥；不一致则要对比一下哪个好。（当然很多文章，手动写循环进行自动判断最佳的p,q，既然有现成的方法，为什么要自己写）

这个方法可以自动判断p和q,我们手动传入,max_ar和 max_ma，即p和q搜索的最大值，前面我们看图，心中预估出了最大值，就设得比预估值再稍微大一点，设太长搜索的时间会久一些，不过对比机器学习训练模型的时间，那是快得多。

from statsmodels.tsa.stattools import arma_order_select_ic
bic_min_order = arma_order_select_ic(data_diff, max_ar=7, max_ma=18,
                                     ic=["aic", "bic"],trend='c')
print(bic_min_order)

'aic_min_order': (6, 16), 'bic_min_order': (1, 1)}

根据两个评价标准，给了两个最佳结果。

4.3模型判断的评价标准

①信息准则 AIC

AIC（Akaike information criterion，赤池信息量准则）：

$AIC=2K-2 \ln(L)$

其中，k表示模型参数的个数(p+q+d),

L表示似然函数,具体数值，在稍微下面一点的模型报告中有写，具体公式不太需要深究，主要是包含了一个负的SSE（sum of squared error）,可知-SSE越负，则整体的AIC越大，模型越差，如果模型效果比较好，则-SSE应该负得不多，AIC越小，模型效果越好，有点信息熵的感觉。

②BIC

AIC的改进版本，考虑到了样本数量影响，因为AIC在样本量大的场景下，通常愿意选择更复杂的模型，即便实际模型应该是一个复杂度很低的，因为AIC惩罚项为参数个数和样本容量没关系，但是似然函数却会受到样本容量的放大。

$BIC=k\ln(n)-2\ln(L)$

其中n为样本容量，其他同AIC，可以消除样本量大时，模型太过复杂的问题。

4.4 建模+检查

#  'aic_min_order': (6, 16), 'bic_min_order': (1, 1)}
# from statsmodels.tsa.arima_model import ARMA 这个方法已经被放弃了，用不了
from statsmodels.tsa.arima.model import ARIMA
# (1, 1)
model_easy = ARIMA(data_diff,order=(1,0,1)).fit()
print(model_easy.summary())
# (6, 16)
model_hard = ARIMA(data_diff,order=(6,0,16)).fit()
print(model_hard.summary())

ARMA已经被放弃了，用不了，我们用ARIMA，将I设为0，也是一样。

此时我们将两组系数都放模型里，查看模型报告！

看这个Log Likelihood,AIC=2L-2In(L)=2*2+2*299.647=607.294;

BIC就是把观察数加了进去；

想知道具体每个含义每个公式，请自己找找相关的文章，其实也没有太大必要每个的公式、计算原理都搞清楚。

重点看圈起来的部分，其中p>z的绝对值部分，说明这是用的Z检验，原假设是该系数为0，p如果小于0.05，则证明该系数应该不为0，则没问题。如果大于0.05，则证明该系数很可能为0，或者接近0，这个时候就要怀疑，这个参数是不是多余了，模型是不是过拟合了，此时我们闷不做声，知道有这个潜在隐患，其对应的coef项也可知，系数很接近0；

双横线下面的是对残差部分的检验，我们希望左上的P >0.05(原假设是残差为白噪声，不能拒绝原假设则证明残差为白噪声，有用的信息都被获取)；我们希望右上角的p>0.05(原假设残差是符合正态分布，如果p>0.05,则残差是符合正态的，此处不符合要求）；我们希望左下角的p>0.05(原假设残差的方差不因自变量时间的选取范围改变，如果p>0.05,则没有异方差性，此处不符合要求）；

但时间序列模型，并不是总是那么完美如意的，很多数据本来就是瞎乱波动的，很难碰到完全符合要求的数据，我们的检验报告，至少能充分证明残差是白噪声...

4.5检查残差

理论上，如果模型拟合充分，则残差应该是服从零均值正态分布且完全随机无相关性，所以还需要分析残差的正态性和无关性，以此检验模型拟合是否良好。如果检验未通过，说明模型不合适或者参数没调好。

上面的model.summary()方法给出了冰冷的数字，一般还会进行下面这一步：

①画残差数据图，看残差是否剧烈变化；

②画残差直方图图，看是不是大概正态分布的；

③画残差QQ图，即Quantile-Quantile Plot，即分位数-分位数图，QQ图中若残差基本完全落在45°线上即为符合正态性假设；

④画残差ACF图，看自相关系数是否有超出蓝色范围的。

有两种方法：

4.5.1挨个画

最重要的，是看ACF图有没有出头的！此处有一根柱子冒出很远，则证明模型参数要调整。

import statsmodels.api as sm
from statsmodels.graphics.tsaplots import plot_acf

fig, axs = plt.subplots(2, 2,figsize=(12,10))
fig.subplots_adjust(hspace=0.3)

model_easy.resid.plot(ax=axs[0][0])
axs[0][0].set_title('residual')

model_easy.resid.plot(kind='hist', ax=axs[0][1])
axs[0][1].set_title('histogram')

sm.qqplot(model_easy.resid, line='45', fit=True, ax=axs[1][0])
axs[1][0].set_title('Normal Q-Q')

plot_acf(model_easy.resid, ax=axs[1][1])

plt.show()

4.5.2自带方法全部画

model_hard.plot_diagnostics(figsize=(14,12))
plt.show()

比较推荐下面这个方法，简单，注意两个图，一个是model_easy,一个是model_hard;看起来model_hard的ACF，虽然后面有两个冒头，但没有冒得那么过分，model_hard的表现是要优于model_easy，但我们的模型还有优化空间，那就是带季节性的SARIMA了。

4.6 SARIMA

其实本文所用的数据，有很明显的趋势和周期性波动，用SARIMA更合适！

SARIMA是在ARIMA的类中，添加一个参数seasonal_order=(P,D,Q,S),可参考：

我写的第二部分，不过还是要先把ARIMA的部分搞明白，再去看SARIMA。

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

五、建模后实现预测

先掌握基础，再往上发展，咱们继续ARIMA的部分：

5.1 预测须知

ARIMA模型适合短期预测（时间序列模型基本都不适合长期预测），对于长期预测，预测结果偏向于均值。

对于AR模型来说，预测得越远，均方误差越大，预测能力越弱；

对于MA模型来说，到q步之后就没了预测能力，比如我们的model_hard=（6,0,16），则最多预测后面16步的数据，并且这16个预测值中，越靠后可信程度越低；

但话又说回来，仅依靠一个时间自变量，去预测未来，本身就是不靠谱的，机器学习拼的就是数据集的多少，时间序列只是能实现一个预测，有总比没有强。

5.2 模型预测

在使用ARIMA模型的情况下，我们基本接近最佳答案了，接下来是预测部分：

直接调用模型自带的方法：

方法1：model.forecast(n)

preds = model_hard.forecast(12)
fig,axs = plt.subplots(1,2,figsize=(12,6))

data.plot(ax=axs[0],color='blue')
axs[0].set_title('原数据')

preds.plot(ax=axs[1],color='green')
axs[1].set_title('纯预测部分')

介于q=16,我们预测12个月份的未来数据，看Y轴可知，预测部分数值很小，而原数据都是300+，这是因为我们模型学习的是差分部分的规律，还需要把减去的那部分加回来，进行还原：

fig,axs = plt.subplots(2,1,figsize=(10,6))

data_diff.plot(ax=axs[0],color='blue')
axs[0].set_title('一阶12滞后差分数据')

preds.plot(ax=axs[1],color='green')
axs[1].set_title('纯预测部分')

如上图所示，差分部分的未来规律！

差分还原：

df1 = pd.DataFrame(data)
df2 = pd.DataFrame(preds.values, columns=['predict'],index=preds.index)
final_df = pd.concat([df1, df2], axis=1)
final_df['result'] = final_df['predict'] + final_df['co2'].shift(12)
final_df.iloc[-24:,:]

我们将预测的未来12个月的数据，和现有的最近12个月的数据一相加，即为未来12个月的数据。

就是上图中红色框框两个部分一加，就是未来的绿色框框的数据，使用shift滑动窗口轻松实现。

plt.figure(figsize=(12, 4))
final_df['result'].plot(color='r', label='arima result')
final_df['co2'].plot(color='g', label='co2 data')
plt.legend()
plt.show()

红色部分，即为未来预测数据，看起来有个断点，实际上是因为plot本身是把一个个点连成线，因为有两个部分，所以plot方法就没有连两个边界点。

方法2：model.predict(start,end)

temp_preds = model_hard.predict(len(data_diff),len(data_diff)+12)
plt.figure(figsize=(12, 4))
data_diff.plot(color='g', label='data_diff',alpha=.5)
temp_preds.plot(color='r', label='total_predict',alpha=.5)
plt.legend()
plt.show()

基本同forecast，不过predict可以将初始值设为过去的，可以"展示"学习到的过去的规律，类似机器学习里其他模型拟合的数据。

predict方法中有个参数dynamic(bool)略微要说下,经过四处查找，终于明白其含义；

我们的模型不是线性的嘛：

预测值=常数（如果有）+AR的几个系数*过去的值+MA的几个系数*过去的残差扰动+随机扰动

当我们调用predict方法的时候，可以让模型展示所拟合到的过去的数据。

如果dynamic = False，则可称为静态预测，其中过去的值，会用到数据中真实的值；

如果dynamic = True，则可称为动态预测，其中过去的值，一开始会用到真实的值，再往后，会逐步用模型预测的值当做真实值。

当然在预测未来之时，没有区别，但在展示过去（画模拟拟合真实值情况），会有较大不同；一般静态预测肯定准一些，毕竟都是真实值，动态预测如果一开始有点偏差，后面会稍稍增大。

history_false = model_hard.predict(len(data_diff)-120,len(data_diff)+12,dynamic=False)
history_true = model_hard.predict(len(data_diff)-120,len(data_diff)+12,dynamic=True)
plt.figure(figsize=(12, 4))
data_diff[-120:].plot(color='g', label='data_diff',alpha=.5)
history_false.plot(color='r', label='静态预测',alpha=.5)
history_true.plot(color='purple', label='动态预测',alpha=.5)
plt.legend()
plt.show()

如上图所示，我们的原数据是到21年12月截止，我们用差分后的数据，进行模型训练，然后预测的，我们对差分的数据训练情况画图时，静态预测一般离真实值不会差太远，说白了就是像线性回归一样，拿着系数乘以真实值；但动态预测，后面就会一直预测为均值，因为后面都是用模型预测的数据去乘以系数，会震荡接近均值。这下你知道为什么说，ARIMA模型预测，预测太远了就没有预测能力了吧，这条紫色线直接就均值干到底了。

为什么还要还原回来这么麻烦？

还不是因为我们训练用的数据，是我们手动进行了1阶12滞后的数据，如果直接用ARIMA的话，直接预测即可，但是这个12滞后参数又不能传进去，所以我们只能进行曲线救国。

虽然要先预测再加回来，略显麻烦，但用SARIMA模型，直接用原数据训练，就可省掉这一步骤，所以看完了笔记（一），就可以看笔记（二）

5.3 模型分析

可知我们是通过AIC、BIC，用自带的方法网格搜索到的最佳p,q值，同时结合我们plot_acf,plot_pacf看图前期预估的p,q值定的阶，而一般机器学习都会有train,test数据。

比如我们预先留出已有数据的最近12个月作为测试数据，通过回归的MSE、MAE等来对模型进行预测，判断初次的p,q值效果如何，对比选出了最佳的p,q值，然后就用这组p,q值重新训练全部数据，来预测未来12个月的值，比较符合机器学习的思路。

优点：过去表现得比较好的p,q值，感觉应该不会太差

缺点：但在时间序列里会存在一个问题，由于我们自变量仅有时间一列，假如最近12个月的波动较大，或者最近一年，跟最近的第二年，差异比较大，随便看看股价图、销量走势图，哪里像本数据这样单纯简单无脑，以MSE为标准，定下认为最好的p,q,数据再加上一段后，那个p,q可能又不是最好的，所以用MSE先验证定p、q,再加数据重新训练再预测，这一招虽然是不错的，但它不一定完全对。

不过也可以看看MSE的效果怎么样，一些文章也把MSE作为一个次要的评价标准，但关键依据还是AIC、BIC，归根结底还是自变量太少，稍微的风吹草动带来的影响就比较大。

TravelLight92

关注

23
点赞
踩
30

收藏

觉得还不错? 一键收藏
0
评论
时间序列模型笔记（一）

把所需要掌握的知识点都说到，同时对一些不太会用到的不作过多讲述，但将深入点的知识链接附上;将完整的流程进行展示，带上数据集案例;作为知识分享，不写屎山，让代码简单直观，因为自己看别人文章，一些简单的功能，非要写得花里胡哨秀技术，让人看得脑壳疼。
复制链接

扫一扫