时间序列算法---ARIMA

  现代时间序列分析方法主要有两个不同的方向：一个方向是由外向内的分析视角产生的方法是与确定性因素分解相关的方法；一个方向是由内向外的分析视角产生的方法是时域分析方法。

一、确定性因素分析方法

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

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模型。

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

  一个观察值序列如果被识别为平稳非白噪声序列，接下来我们将通过考察平稳序列样本自相关系数和偏自相关系数选择适合的模型拟合观察值序列。因此模型拟合的第一步是要根据观察值序列的取值求出该序列的样本自相关系数和样本偏自相关系数的值。
ARMA模型定阶的基本原则：

截尾：落在置信区间内（95%的点都符合该规则）
（1）AR§模型
  描述当前值与历史值之间的关系，用变量自身的历史时间数据对自身进行预测。自回归模型必须满足平稳性的要求。

自回归模型的限制：
1、自回归模型是用自身的数据进行预测
2、必须具有平稳性
3、必须具有自相关性，如果自相关系数小于0.5，则不宜采用
4、自回归只适用预测与自身前期相关的现象
（2）MA(q)模型
  移动平均模型关注的是自回归模型中的误差项的累加，移动平均法能有效地消除预测中的随机波动。

（3）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（差分自回归移动平均模型）—无季节效应的非平稳序列

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

2.2.1、建模流程

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

三、python案例

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)

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}')

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())

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()