时间序列预测算法---ARIMA

向云端UP

已于 2025-02-08 14:56:17 修改

阅读量3.7k

点赞数 43

分类专栏：时间序列文章标签： python 时序数据库

于 2024-11-01 21:22:06 首次发布

本文链接：https://blog.csdn.net/u012121721/article/details/143331196

版权

时间序列专栏收录该内容

12 篇文章

订阅专栏

在机器学习和深度学习的世界当中，存在众多经典且有效的（时间）序列模型。这些模型通常通过三种方式来建立样本与样本之间的关联：
ARIMA家族算法群：过去影响未来，因此未来的值由过去的值加权求和而成，以此构建样本与样本之间的关联。
循环网络家族：遍历时间点/样本点，将过去的时间上的信息传递存储在中间变量中，传递给下一个时间点，以此构建样本和样本之间的关联。
卷积网络家族：使用卷积核扫描时间点/样本点，将上下文信息通过卷积计算整合到一起，以此构建样本和样本之间的关联。
时间序列相关参考文章：
时间序列预测算法—ARIMA
基于VARMAX模型的多变量时序数据预测
 基于机器学习时序库pmdarima实现时序预测
 时间序列预测算法—Prophet
时间序列预测算法—LSTM
长时间序列预测算法—Informer
时间序列分类任务—tsfresh
有季节效应的非平稳序列分析
 python时间序列处理
 时间序列异常值检测方法
 时间序列异常值处理方法
现代时间序列分析方法主要有两个不同的方向：一个方向是由外向内的分析视角产生的方法是与确定性因素分解相关的方法；一个方向是由内向外的分析视角产生的方法是时域分析方法。

一、确定性因素分析方法

确定性因素分解方法认为所有的序列波动都可以归纳为受到如下四大类因素的综合影响：
长期趋势(Trend)。序列呈现出明显的长期递增或递减的变化趋势。
循环波动(Circle)。序列呈现出从低到高再由高到低的反复循环波动。循环周期可长可短，不一定是固定的。
季节性变化(Season)。序列呈现出和季节变化相关的稳定周期波动。
随机波动(lmmediate)。除了长期趋势、循环波动和季节性变化之外，其他不能用确定性因素解释的序列波动，都属于随机波动。
常用的模型：
加法模型：x=T+C+S+I
乘法模型：x=TxCxSxI

1.1、确定性因素分析方法

指数平滑预测方法
简单指数平滑(平稳序列预测)
Holt两参数指数平滑(趋势序列预测)
HoltWinters三参数指数平滑(周期序列预测）
以X11/X12模型为核心的各种季节调整模型
X11模型是第二次世界大战之后，美国人口普查局委托统计学家进行的基于计算机自动进行的时间序列因素分解方法。1954年，X0版本面世，随后十多年陆续推出新的改进版本。1965年，推出成熟版本X11。
1975年，加拿大统计局将ARIMA模型引入X11模型，开发了X11-ARIMA模型。ARIMA模型可以对序列进行向后预测扩充数据，以保证拟合数据的完整性，弥补了中心移动平均方法的缺陷。
1998年，美国人口普查局开发了X12-ARIMA模型。它的改进是将一些特殊因素作为干预变量引入研究。这些干预变量包括:特殊节假日、固定季节因素、工作日因素、交易日因素、闰年因素，以及研究人员自行定义的任意自变量。
2006年美国人口普查局再次推出更新版本X13-ARIMA-Seats，它是在X12-ARIMA的基础上，增加了seats季节调整方法。

二、时域分析方法

时域分析方法主要是从序列自相关的角度揭示时间序列的发展规律。时域分析方法的理论基础：Wold分解定理和Cramer分解定理。
Wold分解定理证明任何平稳序列都可以分解为确定性序列和随机序列之和；Wold分解定理是现代时间序列分析理论的灵魂，是构造ARMA模型拟合平稳序列的理论基础。基于Wold分解定理可以建立：AR模型（自回归模型）、MA模型（移动平均模型）和ARMA模型（自回归移动平均模型）。
Cramer是Wold的指导老师，Cramer分解定理(1962年)是Wold分解定理的理论推广，它是非平稳序列的分解理论，是构造ARIMA模型的理论基础。基于Cramer分解定理可以建立：ARIMA模型。
AR模型只依赖于自身的过去值，MA模型则依赖于过去的误差项。ARMA模型结合了这两者，既考虑过去值也考虑误差项。因此，ARMA模型可以看作是AR和MA模型的结合体，适用于平稳时间序列分析。ARIMA模型（差分自回归移动平均模型）是在ARMA模型的基础上引入了差分运算，以处理非平稳时间序列。基于以上两种定理，下面主要介绍ARMA模型和ARIMA模型。

2.1、ARMA模型（自回归移动平均模型）

同一现象在不同时间上的相继观察值排列而成的序列，这些观察值是一个个随机变量，所以时间序列是随机变量序列。通常可以分为三大类：白噪声序列、平稳非白噪声序列、非平稳序列。

2.1.1、建模流程

在这里插入图片描述

2.1.2、模型的平稳性检验-----单位根检验

对平稳序列建模首先需要确定序列是平稳的。平稳性检验方法有：图检验法、单位根检验法（DF检验和ADF检验）。
图检验法：平稳时间序列具有常数均值和方差。这意味着平稳序列的时序图应该显示出该序列始终在一个常数值附近波动，而且波动的范围有界的特点。图检验方法主要适用于趋势或周期比较明显的序列，对于趋势或周期不太明显的序列，通过图检验方法来判断序列的平稳性具有一定的主观性。图检验法分为：时序图检验和自相关图检验。
自相关图是一个平面二维坐标悬垂线图，横坐标表示延迟时期数，纵坐标表示自相关系数，悬垂线的长度表示自相关系数的大小。平稳序列通常具有短期相关性，这就是我们利用自相关图进行平稳性判别的标准，该性质用自相关系数来描述就是随着延迟阶数k的增加，平稳序列的自相关系数P会很快地衰减向零；而非平稳序列的自相关系数P衰减向零的速度通常比较慢。
单位根检验：如果序列是平稳的，那么该序列的所有特征根都应该在单位圆内。如果序列有特征根在单位圆上或单位圆外，那么该序列就是非平稳序列。

2.1.3、纯随机性（白噪声）检验-----Q统计量和LB统计量

拿到一个观察值序列之后，首先是判断它的平稳性。通过平稳性检验，序列可以分为平稳序列和非平稳序列两大类。但并不是所有的平稳序列都值得建模，只有那些序列值之间具有密切的相关关系、历史数据对未来的发展有一定影响的序列，才值得我们花时间去挖掘历史数据中的有效信息用来预测序列未来的发展。
如果序列值彼此之间没有任何相关性，那就意味着该序列是一个没有记忆的序列，过去的行为对将来的发展没有丝毫影响，这种序列称为纯随机序列。从统计分析的角度来说，纯随机序列是没有任何分析价值的序列。纯随机序列的性质：纯随机性（各序列值之间没有任何相关关系）和方差齐性。
在这里插入图片描述
假设条件：
原假设:延迟期数小于或等于 m 期的序列值之间相互独立
备择假设:延迟期数小于或等于m 期的序列值之间有相关性
检验统计量：Q统计量和LB统计量

2.1.3、模型的参数（p,q）识别-----ACF、PACF

自相关系数ACF衡量当前时间点上的观测值与任意历史时间点的观测值之间的相关性大小，而偏自相关系数PACF衡量当前时间点上的观测值与任意历史时间点的观测值之间的直接相关性的大小。

什么是直接相关性？
时间序列的基本原则之一是“过去影响未来”，假设时间序列 $ts = [y_{t-2}, y_{t-1},y_t]$ ，那 $y_{t-2}$ 的值是如何影响了 $y_t$ 的值的呢？有多种可能：例如， $y_{t-2}$ 的值直接影响了 $y_t$ 的值（比如，双11当天的销量直接冲击了双11两天后的销量），也有可能 $y_{t-2}$ 影响了 $y_{t-1}$ ，再由 $y_{t-1}$ 将这种影响传递给 $y_t$ （比如，前天的降雨量影响了昨天的降雨量，昨天的降雨量又影响今天的降雨量），还有可能两种影响都发生了，又或许 $y_{t-2}$ 与 $y_t$ 之间根本没有互相影响。在这些多种的影响方式中， $y_{t-2}$ 直接对 $y_t$ 产生的影响越大， $y_{t-2}$ 和 $y_t$ 之间的直接相关性越强。
当我们在计算ACF时，我们计算的是多种影响方式带来的“综合效应”，即ACF并不关心两个观测点之间是以什么方式相关，只关心两个观测点之间的相关程度。但PACF只关心“直接相关”的相关程度，并不关心两个观测点之间多种的、其他类型的可能的链接。因此，在求解PACF的时候，我们需要将“不直接”的相关性从ACF的值中“剔除”。

一个观察值序列如果被识别为平稳非白噪声序列，接下来我们将通过考察平稳序列样本自相关系数和偏自相关系数选择适合的模型拟合观察值序列。因此模型拟合的第一步是要根据观察值序列的取值求出该序列的样本自相关系数和样本偏自相关系数的值。
时间序列模型可分为四类:自回归(AR)模型、移动平均(MA)模型自回归移动平均(ARMA)模型和自回归差分移动平均(ARIMA)模型。
ARMA模型定阶的基本原则：
在这里插入图片描述
当ACF/PACF值在蓝色区域之外时，我们就认为当前滞后程度下的ACF/PACF是统计上显著的值，即这个滞后程度下的序列之间的相关性很大程度上是信任的、不是巧合。通常来说，ACF和PACF图有三种常见的形态：拖尾、截尾、既不拖尾也不截尾。其中，拖尾意味着图像呈现按规律衰减、自相关性呈现逐渐减弱的状态（图左）。而截尾特指在某一个滞后程度后ACF/PACF断崖式下跌的状态（图中），例如，滞后0、滞后1对应的ACF都很高，滞后2对应的ACF却断崖式下跌，这种情况被我们称之为“1阶截尾”。（截尾：落在置信区间内（95%的点都符合该规则））
在这里插入图片描述
（1）AR（p）模型
AR（p）称之为p阶的自回归模型(AutoRegressive)，基本思想就是根据自身历史行为预测未来行为的模型。模型非常关注时间序列中的值与它之前和之后的值之间是否存在某种相关性，模型依赖于这种相关性运行。自回归模型必须满足平稳性的要求。

不同时间点的标签值之间强相关(highly-correlated)，位于时间点t的标签值一定强烈地受到t之前的标签值的影响。在数学上，这意味着两个时间点的标签值之间的相关系数会较大。
根据时间的基本属性，两个时间点之间相隔越远，相互之间的影响越弱（例如，昨天是否下雨对今天是否下雨的影响很大，但三个月前的某天是否下雨，对明天是否下雨的影响就相对较小）

自回归模型的限制：
1、自回归模型是用自身的数据进行预测
2、必须具有平稳性
3、必须具有自相关性，如果自相关系数小于0.5，则不宜采用
4、自回归只适用预测与自身前期相关的现象
5、一个自回归模型只能得出一个样本的结果。
（2）MA（q）模型
MA(q)称之为q阶移动平均模型（Moving Average Model），模型在意每个时间点上的数值受到了外界偶然因素多大的影响，关注的是自回归模型中的误差项的累加，移动平均法能有效地消除预测中的随机波动。模型的观点着重于某个时间点的值是被难以想象的复杂因素所主导的，而非只受到历史记录影响。一个MA模型同样只能得出一个时间点上的结果。
时间序列的(长期)趋势与时间序列的(短期)波动受不同因素的影响。从长期来看，时间序列可能呈现缓慢上升、缓慢下降或周期性变动的趋势，这种长时间展现出来的周期受历史标签、及标签本身属性的影响。但在一段较短时间内，时间序列一定是相对稳定的，令时间序列在短期内波动的因子不是历史标签、而是不可预料的各种偶然事件。
不同时间点的标签值之间是关联的，但各种偶然事件在不同时间点上产生的影响之间却是相互独立的。这一点与上面的假设是相辅相成的：历史标签影响着时间序列的长期趋势，那历史标签之间自然是相互关联的。相对的，不可预料的偶然事件之间却没有关联。

公式中的 $y_t$ 表示在时间点 $t$ 时的标签值， $\mu$ 表示当前时间序列标签的均值， $\epsilon_{t}$ 则表示在时间点 $t$ 时、不可预料、不可估计的偶然事件的影响， $\epsilon_{t-1}$ 则代表在时间点 $t - 1$ 时不可预料的、不可估计的偶然事件的影响。不难发现，“偶然事件”这一定义与AR模型中白噪音的定义几乎一模一样，它们都是当前时间点下无法被捕捉到的某些影响，因此 $\epsilon$ 在MA模型中就代表当日的“白噪音”，在理想条件下，MA模型规定 $\epsilon$ 应当服从均值为0、标准差为1的正态分布，但在实际计算时，MA模型规定 $\epsilon$ 等同于模型的预测值 $\hat{y}_{t}$ 与真实标签 $y_{t}$ 之间的差值（Residuals），即：
$\epsilon_t = y_t - \hat{y}_t$
由于偶然事件是无法被预料的、偶然事件带来的影响也是无法被预估的，因此MA模型使用预测标签与真实标签之间的差异就来代表“无法被预料、无法被估计、无法被模型捕捉的偶然事件的影响”。MA模型相信这些影响累加起来共同影响下一个时间点的标签值，因此 $y_t$ 等于所有 $\epsilon$ 的线性组合（加权求和）。
（3）ARMA模型
结合了AR（p）和MA(q)模型的特性，既关心过去对未来的影响，也关心每个时间点上的数值受到外界偶然因素的影响，因此可以应对相对复杂的时间序列数据。具体来说，ARMA模型描述当前值与其自身历史值（自回归部分）以及过去误差项（移动平均部分）之间的关系。优势在于它能够有效捕捉时间序列中的相关性和随机性，适用平稳序列建模与预测。

2.1.4、参数估计

模型识别之后，下一步就是要利用序列的观察值确定该模型的口径，即估计模型中未知参数的值。对于ARMA模型，该模型共有p+q+2个未知参数，参数μ是序列均值，通常采用矩估计方法，用样本均值估计总体均值即可得到它的估计值。原 p+q+2个待估参数减少为 p+q+1个。对于 p+q+1个未知参数的估计方法有三种：矩估计、极大似然估计和最小二乘估计。
在这里插入图片描述

2.1.4.1、矩估计

在这里插入图片描述
矩估计方法，尤其是低阶ARMA 模型场合下的矩估计方法具有计算量小、估计思想简单直观，且不需要假设总体分布的优点。但是在这种估计方法中只用到了p+q个样本自相关系数，即样本二阶矩的信息，观察值序列中的其他信息都被忽略了，这导致矩估计方法是一种比较粗糙的估计方法，它的估计精度一般不高，因此常被用于确定极大似然估计和最小二乘估计迭代计算的初始值。

2.1.4.2、极大似然估计

在极大似然准则下，认为样本来自使该样本出现概率最大的总体，因此未知参数的极大似然估计就是使得似然函数(即联合密度函数)达到最大的参数值。
在这里插入图片描述
使用极大似然估计必须已知总体的分布函数，在时间序列分析中，序列的总体分布通常是未知的。为便于分析和计算，通常假设序列服从多元正态分布。
极大似然估计充分利用了每一个观察值所提供的信息，因而它的估计精度高，同时具有估计的一致性、渐近正态性和渐近有效性等许多优良的统计性质，是一种非常优良的参数估计方法。但它的缺点是需要事先假定序列的分布。

2.1.4.3、最小二乘估计

在这里插入图片描述

2.1.5、模型检验

2.1.5.1、模型显著性检验-----残差序列白噪声检验

模型的显著性检验主要是检验模型的有效性。一个模型是否显著有效主要看它提取的信息是否充分。一个好的拟合模型应该能够提取观察值序列中几乎所有的样本相关信息，换言之，拟合残差项中将不再蕴含任何相关信息，即残差序列应该为白噪声序列，这样的模型称为显著有效模型。反之，如果残差序列为非白噪声序列，那就意味着残差序列中还残留着相关信息未被提取，这就说明拟合模型不够有效，通常需要选择其他模型重新拟合。
因此，模型的显著性检验即残差序列的白噪声检验，原假设和备择假设分别为：
在这里插入图片描述
检验统计量为LB(Ljung-Box)：

如果拒绝原假设，就说明残差序列中还残留着相关信息，拟合模型不显著。如果不能拒绝原假设，就认为拟合模型显著有效。

2.1.5.2、参数的显著性检验----T检验

参数的显著性检验就是要检验每一个未知参数是否显著非零。这个检验的目的是使模型精简。如果某个参数不显著非零，即表示该参数所对应的那个自变量对因变量的影响不明显，该自变量就可以从拟合模型中剔除。最终模型将由一系列参数显著非零的自变量表示。
在这里插入图片描述
检验统计量的P值小于α时，拒绝原假设，认为该参数显著非零，如果参数显著性检验不能拒绝原假设，就应该剔除不显著参数，重新拟合结构更精练的模型。

2.1.6、模型优化

若一个拟合模型通过了检验，说明在一定的置信水平下，该模型能有效拟合观察值序列的波动，但这种有效模型并不一定是唯一的。同一个序列可以构造两个甚至更多个拟合模型，多个模型都显著有效，那么到底该选哪个模型用于统计推断呢?为了解决这个问题，引进AIC和BIC 信息准则的概念，进行模型优化。
（1）AIC：赤池信息准则
该准则的指导思想是一个拟合模型的优劣可以从两方面去考察：一方面是常用来衡量拟合程度的似然函数值；另一方面是模型中未知参数的个数。通常似然函数值越大，说明模型拟合效果越好，模型中未知参数个数越多，说明模型中包含的自变量越多，自变量越多，模型变化越灵活，模型拟合的准确度就会越高。模型拟合程度高是我们所希望的，但是我们又不能单纯地以拟合精度来衡量模型的优劣，因为这样势必会导致未知参数的个数越多越好。
未知参数越多，说明模型中自变量越多，未知的风险越多，而且参数越多，参数估计的难度就越大，估计的精度也越差，因此，一个好的拟合模型应该是一个拟合精度和未知参数个数的综合最优配置。
AIC 准则就是在这种考虑下提出的，它是拟合精度和参数个数的加权函数：
AIC=2(k:模型中未知参数个数)-2ln(L:模型的极大似然函数值)，使AIC达最小的模型被认为是最优模型。
在这里插入图片描述
（2）BIC:贝叶斯信息准则
AIC 准则为选择最优模型提供了有效的规则，但也有不足之处。对于一个观察值序列而言，序列越长，相关信息就越分散，要很充分地提取其中的有用信息，或者说要使拟合精度比较高，通常需要包含多个自变量的复杂模型。在AIC 准则中拟合误差提供的信息会因样本容量而放大，但参数个数的惩罚因子和样本容量没关系，它的权重始终是常数 2。因此，在样本容量趋于无穷大时，由 AIC 准则选择的模型不收敛于真实模型，它通常比真实模型所含的未知参数个数要多。
BIC=(k:模型中未知参数的个数)ln(n:样本数量)-2ln(L:模型的极大似然函数值)
BIC 准则对 AIC准则的改进就是将未知参数个数的惩罚权重由常数2变成了样本容量的对数函数 Inn。

2.2、ARIMA（差分自回归移动平均模型）—无季节效应的非平稳序列

ARIMA模型的基本思想是：一个时间点上的标签值既受过去一段时间内的标签值影响，也受过去一段时间内的偶然事件的影响。这就是说，ARIMA模型假设，标签值是围绕着时间的大趋势而波动的，其中趋势是受历史标签影响构成的，波动是受一段时间内的偶然事件影响构成的，且大趋势本身不一定是稳定的。
$y_t = \beta_0 \ +\ \beta_1y_{t-1} \ + \ \beta_2y_{t-2} +\ ...\ \beta_py_{t-p} \ + \ \epsilon_t\ + \ \theta_1\epsilon_{t-1} \ + \ \theta_2\epsilon_{t-2} \ + \ ... \ \theta_q\epsilon_{t-q}$
通过设置p、q、d的值，我们可以运行一些特殊的ARIMA模型，具体地来说：

参数设置	ARIMA模型的实际情况
ARIMA(p,0,0)	ARIMA等同于自回归模型AR
ARIMA(0,0,q)	ARIMA等同于移动平均模型MA
ARIMA(0,0,0)	不建模，输出等于白噪声 $\epsilon_t$ 或等于白噪声向上或向下平移 $\beta_0$ 个单位
ARIMA(0,1,0) 且不设置常数项 $\beta_0$	在白噪声基础上进行了差分处理输出的结果等同于“随机漫步”（Random Walk）

Cramer分解定理说明任何一个序列的波动都可以视为同时受到确定性影响和随机性影响的作用。平稳序列要求这两方面的影响都是稳定的，而非平稳序列产生的机理就在于它所受到的这两方面的影响至少有一方面是不稳定的。Cramer 分解定理则在理论上保证了适当阶数的差分一定可以充分提取确定性信息。
实践中，我们会根据序列的不同特点选择合适的差分方式，常见情况有以下三种：
（1）序列蕴含显著的线性趋势，1阶差分就可以实现趋势平稳。
（2）序列蕴含曲线趋势，通常低阶(2 阶或 3阶)差分就可以提取出曲线趋势的影响。
（3）蕴含固定周期的序列：对蕴含固定周期的序列进行步长为周期长度（1阶12步）的差分运算，通常可以较好地提取周期信息。

2.2.1、建模流程

在这里插入图片描述
ARIMA主要应用于单变量、同方差场合的线性模型。
ARIMA中AR是自回归，p为自回归项；MA为移动平均，q为移动平均项数，d为时间序列成为平稳时所做的差分次数

三、python案例

通过网盘分享的文件：sunspots_data.xlsx
链接: https://pan.baidu.com/s/1NztswBbfPz4OIGUokG3uHg?pwd=ih1v 提取码: ih1v

3.1、加载数据

import numpy as np  # 导入NumPy库，用于数值计算
import pandas as pd  # 导入Pandas库，用于数据处理和分析
import statsmodels.api as sm  # 导入statsmodels库，用于统计建模
import matplotlib.pyplot as plt  # 导入Matplotlib库，用于绘图
from scipy import stats  # 导入SciPy库，用于科学计算和统计分析
from statsmodels.tsa.arima.model import ARIMA  # 导入ARIMA模型
from statsmodels.graphics.api import qqplot  # 导入QQ图绘制函数
from statsmodels.stats.diagnostic import acorr_ljungbox as lb_test  # 导入Ljung-Box检验函数
from statsmodels.tsa.stattools import adfuller  # 导入ADF单位根检验函数

# 读取 Excel 文件
file_path = r'sunspots_data.xlsx'
df = pd.read_excel(file_path)
# 显示前几行数据
# print(df.head())
# 绘制折线图
plt.plot(df['SUNACTIVITY'])
plt.xlabel('Time')
plt.ylabel('SUNACTIVITY')
plt.title('SUNACTIVITY Over Time')
plt.grid(True)
plt.show()

在这里插入图片描述

3.2、平稳性检验

#单位根检验   原假设是存在单位根，数据不平稳
adfuller(df)

在这里插入图片描述
以上返回的结果分别是检验值、p值、使用的滞后数、使用的样本数、在1%、5%、10%显著性水平上的检验值、以及当前数据上使用ARIMA时可能的AIC值。在这些值当中，我们最为关心前2项数值，因此我们可以将前2项打包到函数中，并让函数直接输出检验的结果：

def adfuller_test(df):
    result=adfuller(df)
    labels = ['ADF Test Statistic','p-value']
    #打包标签和结果
    for value,label in zip(result,labels):
        print(label+' : '+str(value) )
    #如果p值小于等于0.05，则说明我们可以拒绝原假设，认为数据是平稳的。如果p值大于0.05，则说明我们需要接受原假设，数据不平稳。
    if result[1] <= 0.05:
        print("拒绝原假设，数据具有平稳性")
    else:
        print("无法拒绝原假设，数据不平稳")   
    return result[1]
pvalue = adfuller_test(df)

在这里插入图片描述

3.3、差分

# 进行一阶差分
df_diff = df.diff(periods=1).dropna()

# 进行 ADF 检验
result = adfuller(df_diff)
print('ADF Statistic:', result[0])
print('p-value:', result[1])
print('Critical Values:')
for key, value in result[4].items():
    print(f'\t{key}: {value}')

在这里插入图片描述

#通过循环，查看究竟多少步差分能够让这一序列满足单位根检验
#如果许多差分都能够满足单位根检验，我们通常选择ADF值最小、或差分后数据方差最小的步数
columns = ["D1","D2","D3","D4","D5","D6","D7","D8","D9","D10","D11","D12","D24","D36"] #14类
pvalues = []
stds = []

for idx,degree in enumerate([*range(1,13),24,36]):
    #打印差分的步数
    print("{}步差分".format(degree))
    
    #计算差分列
    df[columns[idx]] = df["SUNACTIVITY"] - df["SUNACTIVITY"].shift(idx+1)
    
    #进行ADF检验，提取P值，并计算差分列的标准差
    pvalue = adfuller_test(df[columns[idx]].dropna())
    std_ = df[columns[idx]].std()
    
    #保存P值和标准差
    pvalues.append(pvalue)
    stds.append(std_)
    
    #对结果进行打印
    print("差分后数据的标准差为{}".format(std_))
    print("\n")

在这里插入图片描述

3.4、白噪声检验

lb_test(df, return_df=True,lags=5)#Ljungbox检验发现数据不是白噪声序列

在这里插入图片描述

3.5、模型参数识别

（1）acf、pacf图

import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
# 绘制ACF图
plot_acf(df, lags=154, alpha=0.05)
plt.title('ACF')
plt.show()

# 绘制PACF图
plot_pacf(df, lags=152, alpha=0.05)
plt.title('PACF')
plt.show()

在这里插入图片描述
（2）aic、bic （推荐）

import itertools
import pandas as pd
import statsmodels.api as sm

# 定义参数范围
p_values = range(0, 6)  # p取值范围
q_values = range(0, 6)  # q取值范围

# 初始化最佳模型参数和最小AIC、BIC
best_aic = float('inf')
best_bic = float('inf')
best_order = None

# 遍历参数组合
for p, q in itertools.product(p_values, q_values):
    try:
        # 拟合ARIMA模型
        model = sm.tsa.ARIMA(df, order=(p, 1, q))
        results = model.fit()
        
        # 计算AIC和BIC
        aic = results.aic
        bic = results.bic
        
        # 更新最佳参数和最小AIC、BIC
        if aic < best_aic:
            best_aic = aic
            best_order = (p, 1, q)
        
        if bic < best_bic:
            best_bic = bic
        
        print(f'ARIMA({p}, 1, {q}) - AIC: {aic:.2f}, BIC: {bic:.2f}')
        
    except:
        continue

print(f'Best Model: ARIMA{best_order} - AIC: {best_aic:.2f}, BIC: {best_bic:.2f}')

在这里插入图片描述

3.6、模型拟合

arma_mod021 = ARIMA(df, order=(5, 1, 4)).fit()  # 拟合ARIMA模型到太阳黑子数据集上，指定模型阶数为(2, 0, 0)，即自回归阶数为2，差分阶数为0，移动平均阶数为0，并将结果存储在arma_mod20中
print(arma_mod021.summary())

在这里插入图片描述
summary方法返回的第一部分是当前模型最基本的信息和最基本的评估指标。训练的基本信息包括了当前建立的ARIMA模型、建模时间、样本数量等等信息，这一部分是我们在建模后几乎不会去阅读的内容，在statsmodels中存在6种以上不同的协方差计算方法，默认的方法是opg（即等同于皮尔逊相关系数中的协方差计算方法），值得注意的是模型评估指标部分，这里列出的四个指标分别是对数似然估计、AIC、BIC和HQIC四个评估指标，其中对数似然估计Log Likehood是越大越好，而其他三个指标AIC、BIC和HQIC是越小越好。通常来说我们会直接关注AIC的指标，优先选出AIC最低的模型进行进一步的研究和观察。
summary方法返回的第二部分是整张表单中最为重要的部分，这一部分展示了训练好的ARIMA模型中每一个具体的项的值、以及这些项的显著性水平，结合表单中的内容和我们建立的ARIMA模型时设置的p、d、q值，我们可以写出当前ARIMA模型的具体公式。

coef：该项的系数，比如ar.L1的coef就是 $\beta_1$ 的取值，ma.L1的coef就是 $\theta_1$ 的取值。当该项是常数项时，coef就等于该项本身。
std err：当前项的标准误，衡量当前项对整个ARIMA模型的预测损失的影响。一个项的std err越大，说明这个项对当前模型的损失贡献越大，进一步说明模型对这个项的参数估计可能越不准确。
$\boldsymbol{P>|z|}$ ：该项的系数在统计学上是否显著的p值。在这里我们完成的是针对每一个系数的单变量检验。

单变量检验的原假设（H0）：当前系数coef为0（当前项与标签不相关）
单变量检验的备择假设（H1）：当前系数coef不为0（当前项与标签相关）

很明显，对模型来说理想的状态是拒绝原假设、接受备择假设，因此在我们规定的显著性水平 $\alpha$ （0.05）下，如果有任意项的p值大于 $\alpha$ ，则说明当前项的系数没能通过显著性检验，模型对当前项的系数没有足够的信心。

在一个模型当中，如果很多系数的p值都大于0.05，则可能说明模型处于复杂度过高的状态（即，许多项看起来都与标签无关，所以我们可能创造了许多原本不需要的项）。此时我们可以选择更简单的模型。

Summary表单的第三部分是针对残差完成的各项假设检验，残差是真实标签与当前ARIMA的预测标签之间的差异。当预测值与真实值之间的差异是无法捕捉的、完全无规律的随机数，就说明ARIMA模型已经将能够捕捉的规律全部捕捉到了，因此我们希望残差序列 $r$ 是白噪声序列。在时间序列最严格的定义中，白噪声序列是均值为0、方差为特定 $\sigma^2$ 、服从正态分布、且样本与样本之间相互独立的随机序列，因此针对残差序列 $r$ ，statsmodel中的ARIMA自动完成以下三项检验：

Ljung-Box test

Ljung的发音是"杨"，一般在中文环境中我们称该检验为LBQ检验。LBQ检验是检查样本之间是否存在自相关性的检验。

LBQ检验原假设(H0)：样本之间是相互独立的
LBQ检验备择假设(H1)：样本之间不是相互独立的，存在一些内在联系

如果残差序列是白噪声，那残差序列的样本之间应当是独立的，因此我们想要的是无法拒绝原假设的结果。故而，在查看LBQ测试的p值时，我们希望这是一个大于显著性水平（0.05）、令我们接受原假设的p值。总结一下，LBQ下的Prob(Q)项大于0.05时模型可用，否则可考虑更换模型的超参数组合。

Heteroskedasticity test

这是大名鼎鼎的异方差检验。在一个序列当中，随机抽出不同的子序列并计算这些子序列的方差，如果每个子序列的方差不一致，则称该序列存在“异方差现象”，否则则称该序列为“同方差序列”。在时间序列中，白噪音应当是同方差序列，因此我们对残差做异方差检验，可以帮助我们进一步确认残差是否高度类似于白噪音数据。在统计学中，我们有多种方法来完成异方差检验，Statsmodels执行的是the White’s test。

异方差检验（White’s test）原假设(H0)：该序列同方差
异方差检验（White’s test）备择假设(H1)：该序列异方差

很明显我们想要的是接受原假设的结果，因此同LBQ测试，我们希望异方差检验的p值是大于显著性水平的（大于0.05）。总结一下，异方差检验的Prob(H)大于0.05时模型可用，否则可考虑更换超参数组合。

Jarque-Bera test

雅克贝拉检验（根据葡萄牙语发音音译）是正态性检验中的一种，用于检验当前序列是否符合正态分布。

雅克-贝拉检验原假设(H0)：序列是正态分布的
雅克-贝拉检验备择假设(H1)：序列不是正态分布的

时间序列中的白噪声是需要满足正态性条件的，因此我们希望检验结果无法拒绝原假设，因此我们希望该检验的p值是大于显著性水平的。总结一下，雅克-贝拉检验的Prob(JB)大于0.05时模型可用，否则可考虑更换超参数组合。

对ARIMA模型而言，最佳的情况是三个检验都无法拒绝原假设，但现实往往不会这么完美。在实际选择模型过程中，如果实在找不到能让三个检验都接受原假设的参数组合，也可以选择没通过的检验中、p值最大的参数组合。

3.7、模型检验

（1）残差检验
残差自相关性检验
目的：检验模型的残差是否呈现白噪声，确保模型已经捕捉到数据中的主要模式。
方法：常用 Ljung-Box Q 检验来判断残差的自相关性。
指标：Ljung-Box Q 统计量和相应的 p 值。若 p 值较大（通常大于 0.05），说明残差序列无显著自相关性，可以认为模型有效。

arma_mod021.resid.plot(figsize=(10,3)) #残差时序图

在这里插入图片描述
残差正态性检验
目的：检查残差是否符合正态分布假设。正态性是很多统计检验的前提，残差符合正态分布有助于模型的稳定性。
方法：常用 Jarque-Bera 检验。
指标：Jarque-Bera (JB) 统计量及其 p 值。若 p 值较大，残差分布近似正态；若 p 值小（通常小于 0.05），则残差不符合正态分布假设。

qqplot(arma_mod021.resid, line="q", fit=True)#QQ图检验是否服从正态分布
print(sm.stats.durbin_watson(arma_mod021.resid.values))#DW检验
lb= lb_test(arma_mod021.resid, return_df=True,lags=5)#LB检验
print(lb)
#上述检验都是为了确保残差数据是随机的，不存在趋势或周期性，残差是白噪声

在这里插入图片描述
残差异方差性检验
目的：检查残差的方差是否随时间变化。若存在异方差性，模型在不同时间段的拟合效果可能不一致。
方法：使用 Heteroskedasticity Test（异方差检验）。
指标：H 统计量和相应的 p 值。若 p 值较大（通常大于 0.05），表明残差方差一致。

残差偏度和峰度检验
目的：进一步描述残差的分布特性。
方法：通过偏度 (Skewness) 和峰度 (Kurtosis) 指标描述残差分布的形状。
指标：
偏度 (Skew)：反映残差分布的对称性。接近 0 表示对称，正值表示右偏，负值表示左偏。
峰度 (Kurtosis)：反映残差分布的尖峰程度。接近 3 表示正态分布，较大值表示重尾分布。

（2）参数的显著性检验
ar.L1 至 ar.L5: 自回归（AR）部分的5个滞后项。可以看到其中 ar.L1、ar.L2、ar.L3 和 ar.L4 的系数在统计上显著。
ma.L1 至 ma.L4: 移动平均（MA）部分的4个滞后项。同样，ma.L1、ma.L2 和 ma.L3 系数显著。
sigma2: 残差的方差，反映模型未解释部分的波动性。

3.8、模型拟合预测

# 进行预测
# 拟合ARIMA模型
model = sm.tsa.ARIMA(df, order=(5, 1, 4))
results = model.fit()
forecast_10_periods = results.predict(start=len(df)+1, end=len(df) + 9)
forecast_10_periods

# 绘制原始数据和预测结果的折线图
plt.figure(figsize=(12, 6))
plt.plot(df, label='origin')
plt.plot(forecast, label='predict_data')
plt.plot(range(len(df)+1, len(df) + 10), forecast_10_periods, label='predict_data_10')
plt.xlabel('time')
plt.ylabel('SUNACTIVITY')
plt.title('origin_and_predict_and_10')
plt.legend()
plt.show()