Python 数学建模——ARMA 时间序列分析

前言

  常见的时间序列分析方法有很多，之前介绍了一个稍微新颖的 Prophet 时间序列分析法，这个方法的好处是可以综合考虑节假日的影响（节假日可能导致异常值的出现），并站在不同的时间跨度上考量周际 $(T=7\mathrm{d})$、月际 $(T=30\mathrm{d})$、年际 $(T=365\mathrm{d})$ 的周期性及其影响。Prophet 还能考虑外生变量的影响，这是它的突出特点。详情请看博客：Python 数学建模——Prophet 时间序列预测_多变量prophet-CSDN博客。
  但是 Prophet 的缺点也很多，例如当时间序列不是“日期”的序列时，Prophet 将稍显逊色。本篇博客我们将介绍一个新的时间序列分析方法——ARMA 时间序列分析的使用方法。

ARMA 的具体原理可以参考其他博客，本篇博客主要介绍 ARMA 的使用方法和 Python 实现。

使用前提

  ARMA 模型进行时间序列分析，适合用于平稳非白噪声序列。对于非平稳序列应该使用 ARIMA 模型（即使用差分等方法，后面会稍作介绍），或者取对数等方式，转化为平稳序列进行分析。
  基于这个适用前提，任何一个时间序列在使用 ARMA 之前，都必须进行平稳性检验和白噪声检验。

平稳性检验

  检验一个时间序列的平稳性，有 $3$ 种可行的检验方式。下面将逐一介绍。其中第三种基于假设检验的定量方法很好地摈弃了主观性，因此应该是平稳性检验的最佳选择。前两种方式可以用于数据可视化，带来有关平稳性的直接视觉冲击。

1. 从时间序列上来看，有明显上升、明显下降趋势的为非平稳序列，稳定在一定范围附近的为平稳序列。
2. 从自相关图、偏自相关图上看，具有较强长期相关性（即大部分自相关系数的绝对值都较大，或者说远离0）的是非平稳序列，比如下面最左边的那张图。相对而言自相关系数的绝对值没那么大的，就具有短期相关性，是平稳序列，如下图右边两张图。

图片来自于参考文献。可以看到最左边图，阴影部分达到了 $\pm 0.75$，自相关系数的绝对值已经很大、严重偏离 $0$ 了，说明后期的数据严重依赖于前期的数据，具有比较长期的相关性，不适合使用 ARMA。然而中间这张图，阴影部分在 $\pm 0.4$ 范围左右，则可以认为没有长期相关性。右边这张图也是没有长期相关性。

3. 基于假设检验的平稳性检验方法，或者说是单位根检验。Python 库提供的基于假设检验的方法statsmodels.tsa.stattools.adfuller，对一个时间序列给出 $p$ 值，若 $p<0.05$，则接受“序列是平稳序列”的原假设。

白噪声检验

  白噪声检验也是直接调用 Python 库提供的statsmodels.stats.diagnostic.acorr_ljungbox方法。这个方法基于假设检验给出一个 $p$ 值，当 $p<0.05$ 时接受“序列是非白噪声数据”的原假设。

用法

  ARMA 模型有两个参数 $p,q$，要确定参数，考虑下面的两个指标：

• AIC：赤池信息量准则（Akaike information criterion），定义为$$\mathrm{AIC}=n{\ln \hat{\sigma }}_{\epsilon }^2+2(p+q+1)$$
• BIC：贝叶斯信息量准则（Bayesian information criterion），定义为$$\mathrm{BIC}=n{\ln \hat{\sigma }}_{\epsilon }^2+(p+q)\ln n$$

  这两个指标都是关于 $p,q$ 的函数，它们是研究者针对 ARMA 模型复杂度设计的惩罚项，用于避免过拟合问题。好的参数 $p,q$ 应该要使得上面两个指标尽可能小。
  确定参数是可以只用 AIC/BIC 其中一个，也可以混合起来用，实际操作中通过穷举一些 $(p,q)$ 值，从中挑选最佳参数。即 $(p,q)={\mathrm{argmin}}_{p,q}\mathrm{AIC}(p,q)$，或者 $(p,q)={\mathrm{argmin}}_{p,q}\mathrm{BIC}(p,q)$

代码实例

  有一份太阳黑子数据，展示了从 $1700$ 到 $1988$ 年每年太阳黑子的观测数量。部分如下所示：

点击下载路径，Ctrl+S下载上表所示sunspots.csv表单。

第一步——平稳性分析

  这是一个“年份”的时间序列，不适合使用 Prophet 时间预测模型。考虑使用 ARMA，先进行平稳性分析。

方法一

基于方法一，可以画出时间序列，代码如下：

import pandas as pd, numpy as np
import statsmodels.api as sm
import matplotlib.pylab as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
plt.rc('axes',unicode_minus=False)
plt.rc('font',family='SimHei'); plt.rc('font',size=16)
d=pd.read_csv('sunspots.csv'); dd=d['counts']
years=d['year'].values.astype(int)
plt.plot(years,dd.values,'-*');
plt.savefig("原始数据折线图.pdf")

  画出时间序列图如下图所示。肉眼观察图像，认为太阳黑子数量稳定在一定范围附近，没有明显上升或者明显下降的趋势，原序列基本平稳。

方法二

  基于方法二，我们还可以画自相关图和偏自相关图。代码如下所示：

# 接着上面的代码
ax1=plt.subplot(121); plot_acf(dd,ax=ax1,title='自相关',lags = len(dd) - 1)
ax2=plt.subplot(122); plot_pacf(dd,ax=ax2,title='偏自相关',lags = len(dd) - 1)
plt.savefig("原序列的自相关与偏自相关图.pdf")
# lags 参数是我自己加的，原代码里面没有，其含义应该是相关图横坐标的范围
# 如果不加 lags 参数，自相关、偏自相关图的横坐标就只到 25 左右，很不能说明问题

  画出结果如下图所示。自相关图蓝色阴影在 $\pm 0.4$ 之内，偏自相关图几乎没有阴影。从而认为原序列没有很长的长期相关性，从而是平稳序列。

方法三

  基于方法三，我们可以调库，使用ADF方法得到置信因子 $p$ 的值：

# 接着上面的代码
from statsmodels.tsa.stattools import adfuller as ADF
print(ADF(dd))
"""
结果如下，返回值依次为adf、pvalue、usedlag、nobs、critical values、icbest、regresults、resstore
(-2.3842262328920083, 0.1462380194095098, 8, 280, {'1%': -3.453922368485787, '5%': -2.871918329081633, '10%': -2.5723001147959184}, 2267.166855421211)

看出 pvalue = 0.1462380194095098
"""

  遗憾的是，这个最严谨的方法认为太阳黑子是非平稳序列，理论上原数据不适合使用 ARMA 预测。对于非平稳序列我们的处理手段是：差分，或者取对数。取对数就是将原数据 $x_i$ 转换为 ${\hat{x}}_i=\log x_i$，而一阶差分则是将原数据 $x_i$ 转换为 ${\hat{x}}_i=x_i-x_{i-1}$。
  ARIMA 相比于 ARMA 的区别，就是引入了 $n$ 阶差分，使得原来的非平稳序列 { x i } \{x_i\} {xi​} 转变为平稳序列 { x ^ i } \{\hat x_i\} {x^i​}，从而继续使用 ARMA 进行时间序列分析。我们看一下太阳黑子一阶差分数据的平稳性：

# 接着上面的代码
print(ADF(dd.diff().dropna()))
"""
(-14.076125927559811, 2.8827295545409215e-26, 7, 280, {'1%': -3.453922368485787, '5%': -2.871918329081633, '10%': -2.5723001147959184}, 2263.2365299490502)

可以看出 pvalue = 2.8827295545409215e-26，相当地平稳
"""

  结果相当平稳，因此可以使用太阳黑子的一阶差分数据进行时间序列分析，分析结果求前缀和得到原来的太阳黑子数据。下面我将分别使用原数据dd（虽然不适合用 ARMA，但还是使用着看看）以及它的一阶差分数据dd.diff().dropna()进行时间序列分析，然后对比它们的结果。

第二步——白噪声分析

  使用下面的代码，分别对原数据、一阶差分数据进行白噪声分析。结果显示，它们都不是白噪声数据。光看着一点，它们都符合 ARMA 使用条件。

# 接着上面的代码
from statsmodels.stats.diagnostic import acorr_ljungbox
print(acorr_ljungbox(dd,lags=1)) # 原数据
# 结果 (array([193.5490947]), array([5.34166945e-44]))
# 其中 pvalue = 5.34166945e-44 < 0.05，确定为非白噪声数据
print(acorr_ljungbox(dd.diff().dropna(),lags=1)) # 一阶差分数据
# 结果 (array([80.51235627]), array([2.88892777e-19]))
# 其中 pvalue = 2.88892777e-19 < 0.05，确定为非白噪声数据

第三步——确定参数

  需要确定参数 $p,q$。按照参考文献中的说法，我们需要把满足 $0\le p,q\le \lfloor n/10\rfloor$ 的参数全部尝试一遍，在这之中找出最好的参数对 $(p,q)$ 。但是我们的数据横跨 $289$ 年，$\lfloor n/10\rfloor =28$，这样时间开销太大了，程序会跑得很慢。
  实际上，$(p,q)$ 太大就没必要尝试下去了。为了更快地跑出结果，我这里尝试 $0\le p,q<6$ 的所有参数（也跑了很久），并且只使用 BIC 指标。
  先确定原数据的参数：

# 接着上面的代码
bic_matrix = []  # BIC矩阵
for p in range(6):
    tmp = []
    for q in range(6):
        try:  # 存在部分报错，所以用try来跳过报错。
            tmp.append(sm.tsa.ARMA(dd, (p, q)).fit().bic)
        except:
            tmp.append(float("inf"))
    bic_matrix.append(tmp)
print(np.argmin(np.array(bic_matrix)))
# 结果为 26，即 p = 4,q = 2

  再试一试一阶差分数据的参数：

# 接着上面的代码
bic_matrix = []  # BIC矩阵
for p in range(6):
    tmp = []
    for q in range(6):
        try:  # 存在部分报错，所以用try来跳过报错。
            tmp.append(sm.tsa.ARMA(dd, (p, q)).fit().bic)
        except:
            tmp.append(float("inf"))
    bic_matrix.append(tmp)
print(np.argmin(np.array(bic_matrix)))
# 结果为 14，即 p = 2,q = 2

  至此，确定原数据 ARMA 的参数 $p=4,q=2$，一阶差分数据的参数 $p=2,q=2$。

第四步——模型构建与检验

检验模型效果

  按照下面的代码训练模型：

# 接着上面的代码

################ 获取原数据的预测值 ################
zmd=sm.tsa.ARMA(dd,(4,2)).fit()
dhat=list(zmd.predict())
#################################################

############## 获取一阶差分数据的预测值 ##############
zmd2=sm.tsa.ARMA(dd.diff.dropna(),(2,2)).fit()
dhat2=list(zmd2.predict())
# 进行前缀和操作，恢复为原数据
dhat2.insert(0,dhat[0])
from itertools import accumulate
dhat2 = list(accumulate(dhat2))
##################################################

  上面获取的数据，是对已知的 $1700-1988$ 年数据的拟合，用于检验模型的正确性的。实际上，在获取zmd对象后，可以调用zmd.summary()对此次 ARMA 模型训练结果进行总结。比如print(zmd.summary())的结果就是：

                              ARMA Model Results                              
==============================================================================
Dep. Variable:                 counts   No. Observations:                  289
Model:                     ARMA(4, 2)   Log Likelihood               -1197.676
Method:                       css-mle   S.D. of innovations             15.159
Date:                Sat, 14 Sep 2024   AIC                           2411.353
Time:                        15:20:12   BIC                           2440.684
Sample:                             0   HQIC                          2423.106
                                                                              
================================================================================
                   coef    std err          z      P>|z|      [0.025      0.975]
--------------------------------------------------------------------------------
const           49.7380      6.211      8.009      0.000      37.566      61.910
ar.L1.counts     2.8101      0.086     32.694      0.000       2.642       2.979
ar.L2.counts    -3.1179      0.218    -14.294      0.000      -3.545      -2.690
ar.L3.counts     1.5248      0.213      7.165      0.000       1.108       1.942
ar.L4.counts    -0.2366      0.080     -2.954      0.003      -0.394      -0.080
ma.L1.counts    -1.6480      0.057    -29.016      0.000      -1.759      -1.537
ma.L2.counts     0.7885      0.055     14.396      0.000       0.681       0.896
                                    Roots                                    
=============================================================================
                  Real          Imaginary           Modulus         Frequency
-----------------------------------------------------------------------------
AR.1            0.8639           -0.5664j            1.0330           -0.0924
AR.2            0.8639           +0.5664j            1.0330            0.0924
AR.3            1.0928           -0.0000j            1.0928           -0.0000
AR.4            3.6248           -0.0000j            3.6248           -0.0000
MA.1            1.0450           -0.4197j            1.1262           -0.0608
MA.2            1.0450           +0.4197j            1.1262            0.0608
-----------------------------------------------------------------------------

  这里面可能有一些比较重要的信息。接下来我们就原数据，原数据的预测数据、一阶差分数据的预测数据进行作图：

# 接着上面的代码
plt.figure(figsize=(10, 6))
plt.plot(years,dd,'-')
plt.plot(years,dhat,'--')
plt.plot(years,dhat2,'--')
plt.rcParams['font.family'] = 'KaiTi'
plt.legend(('原始观测值','原数据预测值','一阶差分数据预测值'))
plt.grid()
plt.show()

  作图结果如下所示。本不适合 ARMA 预测的原数据，预测后与原始观测值非常符合；而适合 ARMA 预测的一阶差分数据，进行前缀和操作后，预测结果反而非常糟糕。这并不是因为一阶差分数据拟合效果太差，而是它在通过前缀和操作变回原数据的过程中，误差不断积累变大，最终导致结果产生严重的偏离。
  手动差分数据不可取，要想对非平稳数据进行时间序列分析，还应用 ARIMA 模型合适。

经过尝试，对于原数据 { x i } \{x_i\} {xi​}，使用 ${\hat{x}}_i=\ln (x+0.02)$ （因为原数据中含有 $0$ 值，不可直接取对数，故加 $0.02$）会通过平稳性分析（$p=0.03817<0.05$）和白噪声分析，获取最佳参数 $p=2,q=4$ 且获取比较良好的拟合效果，有兴趣的读者可以尝试一下。使用这样的 { x ^ i } \{\hat x_i\} {x^i​} 进行预测的结果如下图所示：

预测未来数据

  若要预测接下来的数据，可以用下面的代码（二选一，都行）：

# 从已知数据的后一天开始预测
begin = d.shape[0]
# 预测 5 天
end = begin + 5
dnext = zmd.predict(begin, end)
print(dnext)  #显示预测值
"""结果如下：
289    139.440487
290    153.467594
291    143.350740
292    114.224154
293     76.024077
294     40.746472
"""

  或者可以这样：

# 预测后面 5 天的数据
dnext = zmd.forecast(5)
print(dnext[0])  #显示预测值
"""结果如下：
[139.44048728 153.46759385 143.35074018 114.22415414  76.02407734]
"""

  这里的zmd.forecast得到的数据分别是预测结果、标准误差、置信区间。dnext[0]只打印预测结果。
原理和代码参考文献：python之时间序列算法(ARMA)_python arma-CSDN博客