计量经济学(五)——时间序列分析之ARIMA模型预测
时间序列分析（ARIMA）模型是一种广泛用于预测和分析随时间变化的数据模型。ARIMA模型由自回归（AutoRegressive，AR）、差分（Integrated，I）和移动平均（Moving Average，MA）三部分构成。它通过对过去数据的自回归和移动平均来预测未来数据点，广泛应用于经济学、金融、气象学等领域中的时间序列预测。
| 时间序列 | ARIMA流程 |
|---|---|
 |
 |
一、ARIMA模型的基本结构
RIMA（AutoRegressive Integrated Moving Average）模型通过结合自回归（AR）、差分（I）和移动平均（MA）三种思想，能够有效处理非平稳的时间序列数据。它对时间序列的历史数据进行建模，以预测未来趋势，尤其在金融、经济、气象和其他领域的预测任务中应用广泛。ARIMA模型的关键在于通过差分处理非平稳数据，并结合AR和MA模型捕捉数据之间的线性依赖关系。为了更深入理解其核心结构，我们分别扩展讨论自回归（AR）、移动平均（MA）、差分（I）以及整体的ARIMA模型。
1.1 自回归模型（AR）
自回归模型（AR）假设一个时间序列数据点可以通过之前几个数据点的线性组合来预测。换句话说，序列的当前值可以用过去的值来回归。例如，AR(1)模型（即自回归阶数为1）表示当前值 \(Y_t\) 与前一时刻值 \(Y_{t-1}\) 之间存在线性关系：
其中，\(\phi_1\) 是自回归系数，\(\epsilon_t\) 是误差项。对于更高阶的AR模型，如AR(p)，则当前值 \(Y_t\) 取决于之前 \(p\) 个时刻的值：
自回归模型反映了当前值与历史值之间的自相关性，其核心思想是利用序列中前几期的数据来预测后续的值。AR模型在描述那些自相关性明显的序列时尤其有效，例如经济数据中的GDP、股票价格等。通过对AR模型的参数估计，我们可以确定序列中各滞后期的影响程度。当AR系数$ \phi$ 较大时，表示过去的数据对当前时刻有很大的影响；当 \(\phi\) 较小时，历史数据的影响逐渐减弱。
1.2 移动平均模型（MA）
移动平均模型（MA）与自回归模型不同，它假设序列的当前值是过去误差项的线性组合。误差项即是当前值与真实值之间的差异。MA模型的一个典型应用是在金融时间序列中，来平滑数据中的波动。对于MA(q)模型，其数学表达式为：
其中，\(\epsilon_t\) 是当前时刻的随机误差，\(\theta_1, \theta_2, ..., \theta_q\) 是移动平均系数。
相比AR模型，MA模型主要捕捉的是误差项之间的依赖关系，而不是数据值本身的滞后影响。例如，股票市场的价格波动通常受到不可预见的随机事件影响，MA模型可以通过这些误差项来建模价格的变化规律。在实际应用中，MA模型的优点在于它能够捕捉到序列中的短期波动特征，尤其是噪音较大的数据序列。MA模型的阶数 \(q\) 越高，模型对过去误差的依赖就越强。通过选择适当的 \(q\)，可以在平滑噪声和保留数据中的趋势特征之间取得平衡。
1.3 差分过程（I）
差分过程是将非平稳的时间序列转化为平稳序列的关键步骤。在时间序列分析中，平稳序列是指统计特性（如均值、方差等）随时间不发生变化的序列。许多实际应用中的时间序列（如GDP、股票价格等）往往是非平稳的，表现为随时间呈现上升或下降趋势。这时我们通过差分运算将其转化为平稳序列。
差分的基本形式是计算相邻时刻之间的差值：
这称为一阶差分。如果一阶差分后序列仍然不平稳，可以继续进行二阶差分，公式为：
差分次数 \(d\) 表示进行了多少次差分操作以使序列变得平稳。通常，对于大多数的经济和金融数据，一阶差分已经足够，但某些更复杂的数据可能需要更高阶的差分。
差分过程在ARIMA模型中是至关重要的步骤，它消除了趋势或周期性影响，使得模型能更好地捕捉序列中的本质变化。差分次数\(d\) 越大，模型对趋势的消除能力越强，但也增加了模型的复杂度。
1.4 ARIMA模型
综合AR、I和MA的思想，ARIMA模型将自回归、差分和移动平均部分结合在一起，用于处理具有趋势或非平稳性的数据。ARIMA模型的一般形式为：
在这个表达式中：
- \(\Delta^d Y_t\) 表示进行了 \(d\) 次差分后的平稳序列；
- $ \phi_1, \phi_2, ..., \phi_p $ 是自回归项的系数；
- $ \theta_1, \theta_2, ..., \theta_q$ 是移动平均项的系数；
- $ \epsilon_t$ 是随机误差项。
ARIMA模型通过差分消除非平稳性，用AR部分捕捉数据中的长期依赖关系，用MA部分处理短期随机波动。这种混合模型能够有效应对经济、金融等领域中常见的复杂时间序列数据。通过调节模型的参数 \(p\) 、 $d $ 、 \(q\) ，可以灵活地适应不同特征的时间序列数据。具体来说，\(p\) 和 \(q\) 的值可以通过对序列的自相关函数（ACF）和偏自相关函数（PACF）图进行分析来确定，而 \(d\) 则通过平稳性检验来确定。在实际应用中，模型的选择和参数的设定通常需要不断的实验和检验，以确保模型的准确性。
二、ARIMA模型的检验步骤
建立一个ARIMA模型的过程通常包括以下步骤：
- 数据预处理
在使用很多时间序列模型的时候，如ARMA、ARIMA，都会要求时间序列是平稳的，所以一般在研究一段时间序列的时候，第一步都需要进行平稳性检验，除了用肉眼检测的方法，另外比较常用的严格的统计检验方法就是ADF检验，也叫做单位根检验。时间序列如果存在明显的趋势或季节性，就需要通过差分操作来去除这些非平稳因素。
平稳性检验：最常用的平稳性检验方法是ADF检验（Augmented Dickey-Fuller Test），该检验用于判断时间序列是否存在单位根，即是否非平稳。若序列通过检验，说明它是平稳的；若没有通过，则需要通过差分处理。 - 模型识别与阶数确定
确定ARIMA模型的p、d、q值是关键的一步，可以通过以下步骤进行：
确定d值（差分次数）：通过平稳性检验（ADF检验）或观察时间序列的自相关图（ACF图）和偏自相关图（PACF图）来判断。若序列非平稳，则进行差分处理，直至序列平稳。
确定p值（自回归阶数）：通过观察偏自相关函数（PACF）图来确定。若PACF图在某个滞后阶数上截尾（显著），则这个阶数就是p。
确定q值（移动平均阶数）：通过观察自相关函数（ACF）图来确定。若ACF图在某个滞后阶数上截尾，则这个阶数就是q。 - 模型估计
一旦确定了ARIMA模型的阶数\((p, d, q)\)，可以使用最小二乘法或最大似然法来估计模型中的参数。常用的软件如Python中的statsmodels库可以方便地进行ARIMA模型的参数估计。 - 模型检验
模型的检验包括以下几个方面：
残差分析：检查模型的残差是否满足白噪声的假设，即残差应当是独立同分布的随机噪声。如果残差表现出自相关性，说明模型不够充分，需要进一步修正。
AIC/BIC准则：使用Akaike信息准则（AIC）或贝叶斯信息准则（BIC）来评估模型的优劣，选择AIC或BIC值较小的模型。
Ljung-Box检验：用于检验残差的自相关性。如果Ljung-Box检验结果表明残差具有显著的自相关性，则表明模型需要进一步改进。 - 模型预测
在模型通过了各项检验后，可以使用ARIMA模型对未来的时间序列数据进行预测。预测过程基于估计出的ARIMA模型参数，对未来的数据点进行推算。
三、ARIMA模型的阶数确定方法
在ARIMA模型中，三个关键参数 \(p\)（自回归项阶数）、\(d\)（差分阶数）和 \(q\)（移动平均项阶数）的选择对模型的预测效果至关重要。合理的阶数选择可以保证模型的预测精度，同时避免模型过拟合或欠拟合。为了确定最优的ARIMA模型，常用的方法包括自相关函数（ACF）和偏自相关函数（PACF）图法、AIC/BIC信息准则和模型逐步检验。在扩展的讨论中，我们将详细介绍每种方法的原理及其应用。
3.1 ACF和PACF图法
自相关函数（ACF）图和偏自相关函数（PACF）图是时间序列分析中最常用的工具之一，它们用于描述一个时间序列与其不同滞后期的相关性，从而帮助确定模型中自回归和移动平均项的阶数。
- 自相关函数（ACF）图确定 \(q\) 值
ACF图显示了时间序列与其不同滞后期之间的自相关程度。对于确定 \(q\) 值（移动平均项阶数），我们可以通过观察ACF图的截尾特性来推断 \(q\) 值。通常情况下，当ACF图在滞后 \(q\) 阶后迅速截尾，即自相关系数从某一滞后期之后迅速衰减至接近零，则可以推测 \(q\) 为该截尾点的滞后阶数。 - 偏自相关函数（PACF）图确定 \(p\) 值
偏自相关函数（PACF）用于表示序列在排除所有其他滞后期的影响后，仅与特定滞后期之间的相关性。通过PACF图可以帮助确定自回归项的阶数 \(p\)。通常，PACF图如果在某一滞后期之后突然截尾，即在特定滞后阶数之后的偏自相关系数迅速衰减至零，意味着模型中的自回归项只需要包含这几个滞后项。 - 使用ACF和PACF图的注意事项
尽管ACF和PACF图在初步判断模型的阶数时非常有用，但它们也存在一定的局限性，特别是在面对复杂或高噪声的序列时，截尾特征可能不明显或出现混淆。因此，在实际应用中，ACF和PACF图常与其他模型选择标准（如AIC/BIC准则）结合使用，以确定最优模型。
3.2 AIC/BIC信息准则
为了精确量化模型的优劣，AIC（Akaike信息准则）和BIC（贝叶斯信息准则）是常用的两个评估指标。这两个准则综合考虑了模型的拟合优度和复杂性，帮助选择最优的 \(p\)、\(d\)、\(q\) 组合。
- AIC和BIC的定义
- AIC（Akaike Information Criterion）：由日本统计学家赤池弘次提出，用于权衡模型的拟合效果和复杂性。其计算公式为：\[\text{AIC} = -2\log(L) + 2k \]其中，\(L\) 是最大似然估计值，\(k\) 是模型中的参数数量。AIC值越小，模型越优。
- BIC（Bayesian Information Criterion）：与AIC类似，但更注重惩罚复杂模型。其计算公式为：\[\text{BIC} = -2\log(L) + k\log(n) \]其中，\(n\) 为样本大小，\(k\) 为模型参数数量。与AIC相比，BIC对参数数量的惩罚更加严格，因此在样本量较大时，BIC更倾向于选择简单模型。
- AIC（Akaike Information Criterion）：由日本统计学家赤池弘次提出，用于权衡模型的拟合效果和复杂性。其计算公式为：
- AIC和BIC的使用
通过比较不同ARIMA模型的AIC和BIC值，可以有效确定模型的最优阶数 \(p\)、\(d\)、\(q\)。通常，较低的AIC和BIC值表示模型具有较好的拟合能力和较小的复杂性。值得注意的是，AIC和BIC虽都能用于模型选择，但它们有时会给出不同的结果。在样本量较小时，AIC更受欢迎，而在样本量较大时，BIC则因其更严格的复杂性惩罚更为常用。
3.3 模型逐步检验法
模型逐步检验是一种实用且广泛使用的模型选择策略。在该方法中，研究人员通常从低阶模型开始，如ARIMA(0,1,0)或ARIMA(1,0,0)，并通过逐步增加 \(p\) 和 \(q\) 值来比较不同模型的拟合效果。
- 逐步检验的过程
- 初步模型选择：首先进行平稳性检验（如ADF检验）来确定差分次数 \(d\)，然后通过ACF和PACF图初步确定自回归和移动平均项的阶数。
- 模型调整：在确定了初步模型后，通过逐步调整 \(p\) 和 \(q\) 的值构建一系列ARIMA模型。对于每个模型，计算其AIC和BIC值，并选择使得这些值最小的模型。
- 逐步比较：在逐步增加 \(p\) 和 \(q\) 的过程中，需要不断对比各个模型的AIC/BIC值。如果模型复杂度的增加（即增加 \(p\) 或 \(q\) 的值）导致AIC/BIC值上升，说明进一步增加模型复杂度没有提升模型性能，应停止调整。
- 残差检验：除了AIC/BIC外，还需对模型的残差进行诊断。理想情况下，模型的残差应当是白噪声（无自相关性）。通过Ljung-Box检验或ACF图可以验证残差的独立性和随机性。如果残差存在显著的自相关性，则需要进一步调整模型。
- 优缺点分析
逐步检验法的优点在于其操作简单，适用于经验不足的研究者。然而，该方法也有局限性，如模型选择过程中可能陷入局部最优解。此外，该方法的效率较低，尤其在处理大规模数据或复杂模型时，需要进行大量的计算。为此，在实践中，逐步检验法通常与自动化算法（如网格搜索或BIC/AIC自动选择）结合使用，以提高模型选择的效率。
确定ARIMA模型的 \(p\)、\(d\)、\(q\) 值是模型建模中的核心步骤，直接影响模型的拟合效果和预测精度。ACF和PACF图提供了直观的滞后期相关性分析，AIC和BIC准则为模型的复杂性和拟合优度提供了定量评估标准，而逐步检验法则是将理论与实践相结合的有效方法。在实际应用中，研究人员通常会结合这些方法进行模型选择，以获得最佳的ARIMA模型。
四、案例分析
基于1985-2021年某杂志的销售量，预测某商品的未来五年的销售量。
| 年份 | 年度销量 | 年份 | 年度销量 | 年份 | 年度销量 |
|---|---|---|---|---|---|
| 1985 | 100 | 1997 | 111.9 | 2009 | 159.1 |
| 1986 | 101.6 | 1998 | 122.9 | 2010 | 155.1 |
| 1987 | 103.3 | 1999 | 131.9 | 2011 | 161.2 |
| 1988 | 111.5 | 2000 | 134.2 | 2012 | 171.5 |
| 1989 | 116.5 | 2001 | 131.6 | 2013 | 168.4 |
| 1990 | 120.1 | 2002 | 132.2 | 2014 | 180.4 |
| 1991 | 120.3 | 2003 | 139.8 | 2015 | 201.6 |
| 1992 | 100.6 | 2004 | 142 | 2016 | 218.7 |
| 1993 | 83.6 | 2005 | 140.5 | 2017 | 247 |
| 1994 | 84.7 | 2006 | 153.1 | 2018 | 253.7 |
| 1995 | 88.7 | 2007 | 159.2 | 2019 | 261.4 |
| 1996 | 98.9 | 2008 | 162.3 | 2020 | 273.2 |
| 2021 | 279.4 |
4.1 序列平稳化检验，确定\(d\)值
假定某个时间序列是由一系列随机过程生成的，即假定时间序列的每一个数值都是从一个概率分布中随机得到，如果满足下列条件：
- 均值\(\mu\)是与时间\(t\)无关的常数；
- 方差是与时间$t无关的常数；
- 协方差是只与时间间隔\(K\)有关，与时间\(t\)无关的常数；
则称该随机时间序列是平稳的，而该随机过程是平稳随机过程。
ADF 大致的思想就是基于随机游走的。对于随机游走，可以看到比白噪声平滑很多，并且呈现出一些“趋势性”的感觉。它的均值为0，方差与时间\(t\)有关，它不满足平稳性要求。如果一个时间序列是非平稳的，它常常可以通过取差分的方法而形成平稳序列。
import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller
import matplotlib.pyplot as plt
# 创建年度销量数据
data = {
'Year': [
1985, 1986, 1987, 1988, 1989, 1990, 1991, 1992, 1993, 1994,
1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004,
2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014,
2015, 2016, 2017, 2018, 2019, 2020, 2021
],
'Sales': [
100, 101.6, 103.3, 111.5, 116.5, 120.1, 120.3, 100.6, 83.6, 84.7,
88.7, 98.9, 111.9, 122.9, 131.9, 134.2, 131.6, 132.2, 139.8, 142,
140.5, 153.1, 159.2, 162.3, 159.1, 155.1, 161.2, 171.5, 168.4, 180.4,
201.6, 218.7, 247, 253.7, 261.4, 273.2, 279.4
]
}
# 将数据转化为Pandas DataFrame
df = pd.DataFrame(data)
df.set_index('Year', inplace=True)
# 绘制原始时间序列图
plt.figure(figsize=(10, 6))
plt.plot(df.index, df['Sales'], marker='o', color='b', label='Original Sales')
plt.title('Original Time Series')
plt.xlabel('Year')
plt.ylabel('Sales')
plt.grid(True)
plt.legend()
plt.show()
# 执行ADF单位根检验 (原始序列)
result = adfuller(df['Sales'])
output = pd.Series(result[:4], index=['Test Statistic', 'p-value', 'Lags Used', 'Number of Observations'])
output['Critical Values'] = result[4]
print("ADF检验结果 (原始序列):")
print(output)
if result[1] < 0.05:
print("\n该序列没有单位根，可以认为是平稳的。")
else:
print("\n该序列存在单位根，非平稳。")
# 对数据进行差分操作
df_diff = df['Sales'].diff().dropna()
# 绘制差分后时间序列图
plt.figure(figsize=(10, 6))
plt.plot(df_diff.index, df_diff, marker='o', color='r', label='Differenced Sales')
plt.title('Differenced Time Series')
plt.xlabel('Year')
plt.ylabel('Differenced Sales')
plt.grid(True)
plt.legend()
plt.show()
# 对差分后的数据进行ADF单位根检验
result_diff = adfuller(df_diff)
output_diff = pd.Series(result_diff[:4], index=['Test Statistic', 'p-value', 'Lags Used', 'Number of Observations'])
output_diff['Critical Values'] = result_diff[4]
print("\nADF检验结果 (差分序列):")
print(output_diff)
if result_diff[1] < 0.05:
print("\n该差分序列没有单位根，可以认为是平稳的。")
else:
print("\n该差分序列存在单位根，非平稳。")
ADF检验结果 (原始序列):
Test Statistic 1.813771
p-value 0.998376
Lags Used 10.0
Number of Observations 26.0
Critical Values {'1%': -3.7112123008648155, '5%': -2.981246804...
dtype: object
该序列存在单位根，非平稳。
ADF检验结果 (差分序列):
Test Statistic -3.156056
p-value 0.022673
Lags Used 0.0
Number of Observations 35.0
Critical Values {'1%': -3.6327426647230316, '5%': -2.948510204...
dtype: object
该差分序列没有单位根，可以认为是平稳的。
所以差分阶数\(d=1\)。
| 原时间序列 | 一阶差分后的时间序列 |
|---|---|
 |
 |
4.2 确定p值和q值
绘制ACF 、PACF图
拖尾和截尾：拖尾就是序列缓慢衰减，“尾巴”慢慢拖着滑下来，或者震荡衰减；截尾则是突然截断了，像个悬崖，指序列从某个时点变得非常小。专业点来说呢，如果样本自相关系数和样本偏自相关系数在最初的阶明显大于2倍标准差（下图虚线），而后几乎95%的系数都落在2倍标准差的范围内，且非零系数衰减为小值波动的过程非常突然，通常视为k阶截尾。如果有超过5%的样本相关系数大于2倍标准差，或者非零系数衰减为小值波动的过程比较缓慢或连续，通常视为拖尾。
自相关函数 (ACF) 是序列与其滞后自身之间的相关性。假设序列为\(\{X_t\}\)，其滞后\(k\)的自相关系数\(\rho_k\)，定义为：
\[\rho_k = \frac{\sum_{t=1}^{n-k} (X_t - \bar{X})(X_{t+k} - \bar{X})}{\sum_{t=1}^{n} (X_t - \bar{X})^2} \]其中，\(\bar{X}\)为序列的均值。
偏自相关函数 (PACF) 是在控制了所有较小滞后影响后，序列与其滞后项的相关性。其表达较为复杂，一般通过递归算法计算。对于滞后\(k\)的偏自相关系数\(\phi_k\)，它代表序列与滞后\(k\)项的条件自相关，剔除了1至\(k-1\)项的影响。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
# 创建年度销量数据
data = {
'Year': [
1985, 1986, 1987, 1988, 1989, 1990, 1991, 1992, 1993, 1994,
1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004,
2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014,
2015, 2016, 2017, 2018, 2019, 2020, 2021
],
'Sales': [
100, 101.6, 103.3, 111.5, 116.5, 120.1, 120.3, 100.6, 83.6, 84.7,
88.7, 98.9, 111.9, 122.9, 131.9, 134.2, 131.6, 132.2, 139.8, 142,
140.5, 153.1, 159.2, 162.3, 159.1, 155.1, 161.2, 171.5, 168.4, 180.4,
201.6, 218.7, 247, 253.7, 261.4, 273.2, 279.4
]
}
# 将数据转化为Pandas DataFrame
df = pd.DataFrame(data)
df.set_index('Year', inplace=True)
# 对数据进行一阶差分操作
df_diff = df['Sales'].diff().dropna()
# 绘制差分后时间序列的自相关图（ACF）和偏自相关图（PACF）
plt.figure(figsize=(12, 6))
# 自相关图 (ACF)
plt.subplot(121)
plot_acf(df_diff, ax=plt.gca(), title='Autocorrelation (ACF) of Differenced Series')
# 偏自相关图 (PACF)
plt.subplot(122)
plot_pacf(df_diff, ax=plt.gca(), title='Partial Autocorrelation (PACF) of Differenced Series')
plt.tight_layout()
plt.show()
从上图可以看到：趋势序列 ACF 有 1 阶截尾，PACF 有 1 阶截尾尾。因此可以选 p=1， q=1。通过拖尾和截尾对模型定阶，具有很强的主观性。
AIC、BIC准则
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
# 创建年度销量数据
data = {
'Year': [
1985, 1986, 1987, 1988, 1989, 1990, 1991, 1992, 1993, 1994,
1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004,
2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014,
2015, 2016, 2017, 2018, 2019, 2020, 2021
],
'Sales': [
100, 101.6, 103.3, 111.5, 116.5, 120.1, 120.3, 100.6, 83.6, 84.7,
88.7, 98.9, 111.9, 122.9, 131.9, 134.2, 131.6, 132.2, 139.8, 142,
140.5, 153.1, 159.2, 162.3, 159.1, 155.1, 161.2, 171.5, 168.4, 180.4,
201.6, 218.7, 247, 253.7, 261.4, 273.2, 279.4
]
}
# 将数据转化为Pandas DataFrame
df = pd.DataFrame(data)
df.set_index('Year', inplace=True)
# 对数据进行一阶差分操作
df_diff = df['Sales'].diff().dropna()
# 准备存储AIC和BIC值
aic_values = []
bic_values = []
p_values = range(0, 6) # 选择不同的AR阶数
q_values = range(0, 6) # 选择不同的MA阶数
# 遍历不同的AR和MA阶数，拟合ARIMA模型并记录AIC和BIC值
for p in p_values:
for q in q_values:
try:
model = ARIMA(df_diff, order=(p, 1, q)) # 这里d=1代表差分
results = model.fit()
aic_values.append((p, q, results.aic))
bic_values.append((p, q, results.bic))
except:
continue
# 将AIC和BIC值转化为数组
aic_values = np.array(aic_values)
bic_values = np.array(bic_values)
# 找到AIC最小值对应的(p, q)
aic_min_index = np.argmin(aic_values[:, 2])
p_aic_min, q_aic_min, aic_min_value = aic_values[aic_min_index]
# 找到BIC最小值对应的(p, q)
bic_min_index = np.argmin(bic_values[:, 2])
p_bic_min, q_bic_min, bic_min_value = bic_values[bic_min_index]
# 输出最优的(p, q)值和AIC/BIC最小值
print(f"AIC最小值对应的p, q为: p={int(p_aic_min)}, q={int(q_aic_min)}, AIC={aic_min_value}")
print(f"BIC最小值对应的p, q为: p={int(p_bic_min)}, q={int(q_bic_min)}, BIC={bic_min_value}")
# 绘制AIC/BIC图
plt.figure(figsize=(12, 6))
# 绘制AIC图
plt.subplot(121)
plt.scatter(aic_values[:, 0], aic_values[:, 1], c=aic_values[:, 2], cmap='viridis')
plt.colorbar(label='AIC Value')
plt.title('AIC Criterion for ARIMA(p, d=1, q)')
plt.xlabel('AR order (p)')
plt.ylabel('MA order (q)')
# 绘制BIC图
plt.subplot(122)
plt.scatter(bic_values[:, 0], bic_values[:, 1], c=bic_values[:, 2], cmap='plasma')
plt.colorbar(label='BIC Value')
plt.title('BIC Criterion for ARIMA(p, d=1, q)')
plt.xlabel('AR order (p)')
plt.ylabel('MA order (q)')
plt.tight_layout()
plt.show()
从评价准则BIC的结果看：p = 0， q = 2时，两者值最小，BIC最小值为252.58。
ARIMA模型 (p,d,q)
由上述步骤，我们已知d=1，p=0，q=2，故拟合模型为ARIMA（0,1,2），你和后ARIMA模型的表达式(见下面程序)为：
ΔY_t = 0.74 + 0.08 * ε_(t-1) + 61.68 * ε_(t-2) + ε_t。
4.3 拟合预测
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
import numpy as np
# 创建年度销量数据
data = {
'Year': [
1985, 1986, 1987, 1988, 1989, 1990, 1991, 1992, 1993, 1994,
1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004,
2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014,
2015, 2016, 2017, 2018, 2019, 2020, 2021
],
'Sales': [
100, 101.6, 103.3, 111.5, 116.5, 120.1, 120.3, 100.6, 83.6, 84.7,
88.7, 98.9, 111.9, 122.9, 131.9, 134.2, 131.6, 132.2, 139.8, 142,
140.5, 153.1, 159.2, 162.3, 159.1, 155.1, 161.2, 171.5, 168.4, 180.4,
201.6, 218.7, 247, 253.7, 261.4, 273.2, 279.4
]
}
# 将数据转化为Pandas DataFrame
df = pd.DataFrame(data)
df.set_index('Year', inplace=True)
# 创建ARIMA(0,1,2)模型并进行拟合
model = ARIMA(df['Sales'], order=(0, 1, 2))
model_fit = model.fit()
# 输出模型摘要信息
print(model_fit.summary())
# 获取ARIMA模型的参数
params = model_fit.params
print("\nARIMA模型的表达式为：")
print(f"ΔY_t = {params[0]:.2f} + {params[1]:.2f} * ε_(t-1) + {params[2]:.2f} * ε_(t-2) + ε_t")
# 预测未来5年，即从2022年到2026年的销量
forecast_years = 5
forecast = model_fit.forecast(steps=forecast_years).round(2) # 保留两位小数
# 生成预测的年份
future_years = np.arange(2022, 2022 + forecast_years)
# 创建包含实际值和预测值的新DataFrame
future_sales = pd.DataFrame({'Year': future_years, 'Forecasted_Sales': forecast})
future_sales.set_index('Year', inplace=True)
# 将预测值添加到原始数据中，使用pd.concat替代append
df_forecasted = pd.concat([df, future_sales])
# 绘制实际值和预测值的对比图
plt.figure(figsize=(10, 6))
plt.plot(df_forecasted.index, df_forecasted['Sales'], label='Actual Sales', marker='o')
plt.plot(future_sales.index, future_sales['Forecasted_Sales'], label='Forecasted Sales', linestyle='--', marker='x')
plt.title('Sales Forecast for 2022-2026 (ARIMA(0,1,2))')
plt.xlabel('Year')
plt.ylabel('Sales')
plt.legend()
plt.grid(True)
plt.show()
# 输出未来5年的预测值，保留两位小数
print("预测的2022-2026年销量（保留两位小数）:")
print(future_sales.round(2))
SARIMAX Results
==============================================================================
Dep. Variable: Sales No. Observations: 37
Model: ARIMA(0, 1, 2) Log Likelihood -125.601
Date: Tue, 15 Oct 2024 AIC 257.202
Time: 18:21:44 BIC 261.952
Sample: 0 HQIC 258.860
- 37
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ma.L1 0.7433 0.161 4.625 0.000 0.428 1.058
ma.L2 0.0831 0.202 0.411 0.681 -0.313 0.479
sigma2 61.6759 13.336 4.625 0.000 35.537 87.815
===================================================================================
Ljung-Box (L1) (Q): 0.29 Jarque-Bera (JB): 1.74
Prob(Q): 0.59 Prob(JB): 0.42
Heteroskedasticity (H): 2.30 Skew: -0.14
Prob(H) (two-sided): 0.16 Kurtosis: 4.04
===================================================================================
ARIMA模型的表达式为：
ΔY_t = 0.74 + 0.08 * ε_(t-1) + 61.68 * ε_(t-2) + ε_t
预测的2022-2026年销量（保留两位小数）:
Forecasted_Sales
Year
2022 281.85
2023 282.10
2024 282.10
2025 282.10
2026 282.10
总结
ARIMA模型能够处理具有不同趋势和波动特征的时间序列数据，适用于很多实际应用场景。相比之下，简单的线性回归模型无法处理时间序列中的自相关性，而更加复杂的结构化模型如状态空间模型则往往需要更高的计算复杂度。此外，ARIMA模型在处理非平稳数据时具有独特的优势，通过差分操作可以有效应对数据中的趋势和季节性变化。然而，ARIMA模型的局限性在于其假设时间序列的线性关系，因此对于高度非线性的数据，可能需要扩展模型，如SARIMA（季节性ARIMA）或引入非线性建模工具（如GARCH模型）来应对复杂的时间序列特征。ARIMA模型通过自回归、差分和移动平均的有机结合，能够很好地对非平稳时间序列进行建模和预测，尤其在经济和金融领域具有广泛的应用。
