#Python 时间序列分析：基于 statsmodels 的 ARMA/ARIMA 自动参数寻优与预测

Python 的科学计算和数据挖掘相关库中，pandas 和 statsmodels 都提供了时间序列相关分析功能。本示例使用的是 statsmodels 做时间序列预测应用。有关时间序列算法的选择，实际场景中最常用的是 ARIMA 或 ARMA，因此本示例将使用 ARIMA/ARMA 来做时间序列分析。

对于这两种时间序列方法而言，应用的难点在于如何根据不同的场景判断参数值（即 p、d、q）。本示例将设置判断阈值，通过自动化的程序方式来完成自动的 ARIMA/ARMA 参数选择以及模型训练，大幅降低时间序列算法的应用难度。

背景说明：示例中模拟的是针对具有时间序列特征的数据集做未来时间序列的预测。数据源文件 time_series.txt 可以从《Python数据分析与数据化运营》的“附件-chapter4”中找到。

---

#一、 核心功能模块封装

为了实现自动化的时间序列分析，我们将整个流程拆解为多个功能函数。

#1. 导入依赖库

statsmodels 库里面的 plot_acf、plot_pacf 用于 ACF 和 PACF 图形检验分析，adfuller、acorr_ljungbox 用于 ADF 检验和随机性检验，ARMA 库做时间序列分析，matplotlib 和 prettytable 做图形和表格格式化输出。

pythonCopy
1# 导入库
2import pandas as pd # pandas库
3import numpy as np # numpy库
4from statsmodels.graphics.tsaplots import plot_acf, plot_pacf # acf和pacf展示库
5from statsmodels.tsa.stattools import adfuller # adf检验库
6from statsmodels.stats.diagnostic import acorr_ljungbox # 随机性检验库
7from statsmodels.tsa.arima_model import ARMA # ARMA库
8import matplotlib.pyplot as plt # matplotlib图形展示库
9import prettytable # 导入表格库

💡 提示：ARIMA 与 ARMA 的选择 ARIMA 相对于 ARMA 多了一步差分的过程。但是由于 statsmodels 中的 ARIMA 对差分的支持不太好，最多只能做 2 阶差分，其实用性会大打折扣（其实 2 阶差分已经可以满足 90% 以上的场景了）。所以在实际应用中，我们通常会手动做数据平稳性处理，然后将平稳数据传给 ARMA 模型做训练，以此完成 ARIMA 的时间序列分析。

#2. 格式化表格输出：pre_table

该函数用来展示格式化表格数据，方便每次展示时调用而无需重复编写代码。

参数说明：

• table_name：表格名称，字符串列表。
• table_rows：表格内容，嵌套列表。考虑到可能会有批量添加的模式，数据记录行需要转换为嵌套列表的形式（例如 [[1,2,3,4], [2,2,3,4]]；即使只有一条记录，也要写成 [[1,2,3,4]]）。

pythonCopy
1# 多次用到的表格
2def pre_table(table_name, table_rows):
3    '''
4    :param table_name: 表格名称，字符串列表
5    :param table_rows: 表格内容，嵌套列表
6    :return: 展示表格对象
7    '''
8    table = prettytable.PrettyTable() # 创建表格实例
9    table.field_names = table_name # 定义表格列名
10    for i in table_rows: # 循环读多条数据
11        table.add_row(i) # 增加数据
12    return table

#3. 数据平稳性处理：get_best_log

该函数用于数据平稳处理。先通过 rule1 和 rule2 判断输入的两个规则是否同时满足，如果满足则说明时间序列是平稳的，无需处理；否则使用对数方法（numpy.log）做平稳性处理。

参数说明：

• ts：时间序列数据，Series 类型。
• max_log：最大 log 处理的次数，默认值为 5，int 型。
• rule1：rule1 规则布尔值，默认值为 True。
• rule2：rule2 规则布尔值，默认值为 True。

pythonCopy
1# 数据平稳处理
2def get_best_log(ts, max_log=5, rule1=True, rule2=True):
3    '''
4    :param ts: 时间序列数据，Series类型
5    :param max_log: 最大log处理的次数，int型
6    :param rule1: rule1规则布尔值，布尔型
7    :param rule2: rule2规则布尔值，布尔型
8    :return: 达到平稳处理的最佳次数值和处理后的时间序列
9    '''
10    if rule1 and rule2: # 如果两个规则同时满足
11        return 0, ts # 直接返回0和原始时间序列数据
12    else: # 只要有一个规则不满足
13        for i in range(1, max_log): # 循环做log处理
14            ts = np.log(ts) # log处理
15            lbvalue, pvalue1 = acorr_ljungbox(ts, lags=1) # 白噪声检验结果
16            adf, pvalue2, usedlag, nobs, critical_values, icbest = adfuller(ts) # ADF检验
17            rule_1 = (adf < critical_values['1%'] and adf < critical_values['5%'] and adf < critical_values['10%'] and pvalue1 < 0.01) # 稳定性检验
18            rule_2 = (pvalue2 < 0.05) # 白噪声检验
19            rule_3 = (i < 5)
20            if rule_1 and rule_2 and rule_3: # 如果同时满足条件
21                print('The best log n is: {0}'.format(i)) # 打印输出最佳次数
22                return i, ts # 返回最佳次数和处理后的时间序列

💡 知识点拓展：时间序列数据的平稳性 平稳性是做时间序列分析的前提条件。通俗理解就是数据没有随着时间呈现明显的趋势和规律（如剧烈波动、递增、递减等），而是相对均匀且随机地分布在均值附近。在 ARIMA 模型中的 I 就是对数据做差分以实现平稳，关键参数 p、d、q 中的 d 即时间序列成为平稳时所做的差分次数。

如何判断是否需要平稳性处理？

1. 观察法：通过输出时间序列图发现数据是否平稳。
2. 自相关和偏相关法：通过观察自相关和偏相关的系数分析数据是否平稳。
3. ADF 检验：通过 ADF 检验得到的显著性水平分析数据是否平稳。

实现数据平稳的常用方法：

• 对数法：减小数据的波动，使其线性规律更加明显（仅限值大于 0 的数据集）。
• 差分法：非纯随机的时间序列经一阶或二阶差分后通常会变得平稳。
• 平滑法：分为移动平均法和指数平均法。
• 分解法：将时序数据分离成成长期趋势、季节趋势和随机成分等。

#4. 平稳化数据还原：recover_log

在经过平稳性处理之后，原有的时间序列训练集已经被更改，通过 ARMA 得到的预测数据也是经过平稳处理后的数据，必须经过还原才能投入实际使用。

pythonCopy
1# 还原经过平稳处理的数据
2def recover_log(ts, log_n):
3    '''
4    :param ts: 经过log方法平稳处理的时间序列，Series类型
5    :param log_n: log方法处理的次数，int型
6    :return: 还原后的时间序列
7    '''
8    for i in range(1, log_n + 1): # 循环多次
9        ts = np.exp(ts) # log方法还原
10    return ts # 返回时间序列

#5. 平稳性检验：adf_val

该函数实现了时间序列图、ACF 图、PACF 图、ADF 检验数据等多种检验方法的数据和图形输出。返回的 ADF 值、P 值等将用来构成 rule1 的判断规则。

pythonCopy
1# 平稳性检验
2def adf_val(ts, ts_title, acf_title, pacf_title):
3    '''
4    :param ts: 时间序列数据，Series类型
5    :param ts_title: 时间序列图的标题名称，字符串
6    :param acf_title: acf图的标题名称，字符串
7    :param pacf_title: pacf图的标题名称，字符串
8    :return: adf值、adf的p值、三种状态的检验值
9    '''
10    plt.figure()
11    plt.plot(ts) # 时间序列图
12    plt.title(ts_title) # 时间序列标题
13    plt.show()
14    
15    plot_acf(ts, lags=20, title=acf_title).show() # 自相关检测
16    plot_pacf(ts, lags=20, title=pacf_title).show() # 偏相关检测
17    
18    adf, pvalue, usedlag, nobs, critical_values, icbest = adfuller(ts) # 平稳性检验
19    table_name = ['adf', 'pvalue', 'usedlag', 'nobs', 'critical_values', 'icbest'] # 表格列名列表
20    table_rows = [[adf, pvalue, usedlag, nobs, critical_values, icbest]] # 表格行数据，嵌套列表
21    adf_table = pre_table(table_name, table_rows) # 获得平稳性展示表格对象
22    
23    print('stochastic score') # 打印标题
24    print(adf_table) # 打印展示表格
25    return adf, pvalue, critical_values # 返回adf值、adf的p值、三种状态的检验值

#6. 白噪声（随机性）检验：acorr_val

该函数用于白噪声检验，返回的 P 值用来做 rule2 判断规则的参数值。

pythonCopy
1# 白噪声（随机性）检验
2def acorr_val(ts):
3    '''
4    :param ts: 时间序列数据，Series类型
5    :return: 白噪声检验的P值和展示数据表格对象
6    '''
7    lbvalue, pvalue = acorr_ljungbox(ts, lags=1) # 白噪声检验
8    table_name = ['lbvalue', 'pvalue'] # 表格列名列表
9    table_rows = [[lbvalue, pvalue]] # 表格行数据，嵌套列表
10    acorr_ljungbox_table = pre_table(table_name, table_rows)
11    
12    print('stationarity score') # 打印标题
13    print(acorr_ljungbox_table) # 打印展示表格
14    return pvalue # 返回白噪声检验的P值和展示数据表格对象

💡 知识点拓展：白噪声检验 白噪声（White Noise）检验也称为随机性检验，用于检验时间序列的各项数值之间是否具有任何相关关系。检测方法包括：

• 图形法：时间序列应该是围绕均值随机性上下分布的状态。
• Ljung-Box 法：对时间序列是否存在滞后相关的一种统计检验。

白噪声检验通常和平稳性检验协同进行，如果平稳性检验通过，白噪声检验一般也会通过。

#7. ARMA 最优模型训练：arma_fit

该模型实现的是基于给定的时间序列数据，自动寻找 BIC 最小时的 p 和 q 的阶数，并形成训练好的 ARMA 模型。

pythonCopy
1# arma最优模型训练
2def arma_fit(ts):
3    '''
4    :param ts: 时间序列数据，Series类型
5    :return: 最优状态下的arma模型对象
6    '''
7    max_count = int(len(ts) / 10) # 最大循环次数最大定义为记录数的10%
8    bic = float('inf') # 初始值为正无穷
9    tmp_score = [] # 临时p、q、aic、bic和hqic的值的列表
10    
11    for tmp_p in range(max_count + 1): # p循环max_count+1次
12        for tmp_q in range(max_count + 1): # q循环max_count+1次
13            model = ARMA(ts, order=(tmp_p, tmp_q)) # 创建ARMA模型对象
14            try:
15                results_ARMA = model.fit(disp=-1, method='css') # ARMA模型训练
16            except:
17                continue # 遇到报错继续
18            finally:
19                try:
20                    tmp_aic = results_ARMA.aic # 模型的获得aic
21                    tmp_bic = results_ARMA.bic # 模型的获得bic
22                    tmp_hqic = results_ARMA.hqic # 模型的获得hqic
23                    tmp_score.append([tmp_p, tmp_q, tmp_aic, tmp_bic, tmp_hqic]) # 追加每个模型的训练参数和结果
24                    if tmp_bic < bic: # 如果模型bic小于最小值，那么获得最优模型ARMA的下列参数：
25                        p = tmp_p # 最优模型ARMA的p值
26                        q = tmp_q # 最优模型ARMA的q值
27                        model_arma = results_ARMA # 最优模型ARMA的模型对象
28                        aic = tmp_aic # 最优模型ARMA的aic
29                        bic = tmp_bic # 最优模型ARMA的bic
30                        hqic = tmp_hqic # 最优模型ARMA的hqic
31                except:
32                    continue
33                    
34    pdq_metrix = np.array(tmp_score) # 将嵌套列表转换为矩阵
35    pdq_pd = pd.DataFrame(pdq_metrix, columns=['p', 'q', 'aic', 'bic', 'hqic']) # 基于矩阵创建数据框
36    table_name = ['p', 'q', 'aic', 'bic', 'hqic'] # 表格列名列表
37    table_rows = [[p, q, aic, bic, hqic]] # 表格行数据，嵌套列表
38    parameter_table = pre_table(table_name, table_rows) # 获得最佳ARMA模型结果展示表格对象
39    
40    print('each p/q traning record') # 打印标题
41    print(pdq_pd) # 打印输出每次ARMA拟合结果，包含p、d、q以及对应的AIC、BIC、HQIC
42    print('best p and q') # 打印标题
43    print(parameter_table) # 输出最佳ARMA模型结果展示表格对象
44    return model_arma # 最优状态下的arma模型对象

⚠️ 注意与知识点：statsmodels 中的 ARMA statsmodels 中的 ARMA 通过应用 fit 方法后，返回的对象是一个 ARMAResults 类，而不再是原有的 model 本身（这一点跟 sklearn 不同）。因此无法直接通过 model 本身获取参数和评估指标，必须通过 ARMAResults 来获取。

关于 p 和 q 的选择：

1. 人工观察法：观察 ACF 和 PACF 图，通过截尾、拖尾等信息分析。
2. 程序遍历法（本例采用）：通过循环遍历所有条件，并依据 AIC、BIC 最小原则自动确定阶数。

#8. 模型训练与效果评估：train_test

该函数可将训练的时间序列内预测值与实际值做比较，并输出 RMSE（均方根误差） 和时间序列对比图。

pythonCopy
1# 模型训练和效果评估
2def train_test(model_arma, ts, log_n, rule1=True, rule2=True):
3    '''
4    :param model_arma: 最优ARMA模型对象
5    :param ts: 时间序列数据，Series类型
6    :param log_n: 平稳性处理的log的次数，int型
7    :param rule1: rule1规则布尔值，布尔型
8    :param rule2: rule2规则布尔值，布尔型
9    :return: 还原后的时间序列
10    '''
11    train_predict = model_arma.predict() # 得到训练集的预测时间序列
12    if not (rule1 and rule2): # 如果两个条件有任意一个不满足
13        train_predict = recover_log(train_predict, log_n) # 恢复平稳性处理前的真实时间序列值
14        ts = recover_log(ts, log_n) # 时间序列还原处理
15        
16    ts_data_new = ts[train_predict.index] # 将原始时间序列数据的长度与预测的周期对齐
17    RMSE = np.sqrt(np.sum((train_predict - ts_data_new) ** 2) / ts_data_new.size) # 求RMSE
18    
19    # 对比训练集的预测和真实数据
20    plt.figure() # 创建画布
21    train_predict.plot(label='predicted data', style='--') # 以虚线展示预测数据
22    ts_data_new.plot(label='raw data') # 以实线展示原始数据
23    plt.legend(loc='best') # 设置图例位置
24    plt.title('raw data and predicted data with RMSE of %.2f' % RMSE) # 设置标题
25    plt.show() # 展示图像
26    return ts # 返回还原后的时间序列

#9. 未来数据预测：predict_data

该函数用来预测未来指定时间项的数据，并将原始时间序列（实线）和预测的时间序列（虚线）输出到一个图像中。

pythonCopy
1# 预测未来指定时间项的数据
2def predict_data(model_arma, ts, log_n, start, end, rule1=True, rule2=True):
3    '''
4    :param model_arma: 最优ARMA模型对象
5    :param ts: 时间序列数据，Series类型
6    :param log_n: 平稳性处理的log的次数，int型
7    :param start: 要预测数据的开始时间索引
8    :param end: 要预测数据的结束时间索引
9    :param rule1: rule1规则布尔值，布尔型
10    :param rule2: rule2规则布尔值，布尔型
11    :return: 无
12    '''
13    predict_ts = model_arma.predict(start=start, end=end) # 预测未来指定时间项的数据
14    print('-----------predict data----------') # 打印标题
15    if not (rule1 and rule2): # 如果两个条件有任意一个不满足
16        predict_ts = recover_log(predict_ts, log_n) # 还原数据
17    print(predict_ts) # 展示预测数据
18    
19    # 展示预测趋势
20    plt.figure() # 创建画布
21    ts.plot(label='raw time series') # 设置推向标签
22    predict_ts.plot(label='predicted data', style='--') # 以虚线展示预测数据
23    plt.legend(loc='best') # 设置图例位置
24    plt.title('predicted time series') # 设置标题
25    plt.show() # 展示图像

---

#二、 实战应用：时间序列预测完整流程

在完成上述核心模块的封装后，我们将按照以下 8 个步骤执行完整的预测流程。

#步骤 1：读取与预处理数据

先通过 lambda 函数定义一个日期解析器，配合 pandas.datetime.strptime 将字符串解析为日期格式。

pythonCopy
1# 读取数据
2date_parse = lambda dates: pd.datetime.strptime(dates, '%m-%d-%Y') # 创建解析列的功能对象
3df = pd.read_table('time_series.txt', delimiter='\t', index_col='date', date_parser=date_parse)
4ts_data = df['number'].astype('float32') # 将列转换为float32类型
5print('data summary') # 打印标题
6print(ts_data.describe()) # 打印输出时间序列数据概况

#步骤 2：原始数据检验

直接调用定义好的稳定性检验和白噪声检验函数，传入原始时间序列数据。

pythonCopy
1# 原始数据检验 
2adf, pvalue1, critical_values = adf_val(ts_data, 'raw time series', 'raw acf', 'raw pacf') # 稳定性检验 
3pvalue2 = acorr_val(ts_data) # 白噪声检验 

原始时间序列图：4 月份之前波动较为明显，4 月份之后基本处于平稳波动状态。

原始 ACF 检验图：前 20 个滞后点基本已经看出明显趋势。

原始 PACF 检验图

ADF 稳定性检验与白噪声检验结果

指标分析：如果时间序列具有稳定性，ADF 的值需要小于 critical_values 中的 1%、5% 和 10% 的三个值，且 P 值一般小于 0.01。在白噪声检验中，P 值需要小于 0.05 即说明数据是随机分布的。从结果看，示例数据恰好满足条件。

#步骤 3：定义平稳性处理规则

基于步骤 2 的返回结果，定义 ADF 的 P 值和白噪声检验的 P 值的阈值规则。

pythonCopy
1# 创建用于区分是否进行平稳性处理的规则 
2rule1 = (adf < critical_values['1%'] and adf < critical_values['5%'] and adf < critical_values['10%'] and pvalue1 < 0.01) # 稳定性检验 
3rule2 = (pvalue2[0,] < 0.05) # 白噪声检验 

#步骤 4：执行平稳性处理

如果数据符合平稳性，则不需要处理；否则循环处理直到满足阈值要求。由于示例数据已满足条件，此处直接返回 0 和原始数据集。

pythonCopy
1# 对时间序列做稳定性处理 
2log_n, ts_data = get_best_log(ts_data, max_log=5, rule1=rule1, rule2=rule2)

#步骤 5：复检数据平稳性

再次做检验。由于没有做任何处理，返回的结果仍然跟步骤 2 相同。

pythonCopy
1adf, pvalue1, critical_values = adf_val(ts_data, 'final time series', 'final acf', 'final pacf') # 稳定性检验 
2pvalue2 = acorr_val(ts_data) # 白噪声检验 

#步骤 6：训练最佳 ARMA 模型

调用 arma_fit 自动寻找最优参数并训练模型。（交互窗口可能会有警告信息，忽略即可）。

pythonCopy
1# 训练最佳ARMA模型并输出相关参数和对象 
2model_arma = arma_fit(ts_data)

输出每次训练不同 p 和 q 对应的 AIC、BIC 和 HQIC 值，最终选出满足 BIC 最小原则的最优参数。

#步骤 7：模型内测与效果评估

通过最优模型对训练集做时间序列内的预测，并与原始数据做对比，输出 RMSE。

pythonCopy
1ts_data = train_test(model_arma, ts_data, log_n, rule1=rule1, rule2=rule2)

ARMA 预测结果与实际结果对比：预测值与实际值差异不大，较好地还原了实际数据分布规律。

#步骤 8：未来时间序列预测

设置预测起止时间（1991-07-28 至 1991-08-02），调用 predict_data 做预测分析。

pythonCopy
1start = '1991-07-28' # 设置预测开始的时间索引 
2end = '1991-08-02' # 设置预测结束的时间索引 
3predict_data(model_arma, ts_data, log_n, start, end, rule1=rule1, rule2=rule2) # 预测时间序列数据 

预测结果与实际结果合成图

---

#三、 关键点与实操小结

在上述完整的预测流程中，核心需要关注以下几个业务逻辑节点：

1. 数据平稳性处理：判断与执行。
2. p 和 q 阶数评估：除了人工观察 ACF/PACF 图，采用 AIC/BIC 遍历是一种高效的自动化方法。
3. 数据还原：在做数据平稳性处理（如 log）之后，务必注意逆向还原数据，否则预测结果将失去业务意义。
4. 模型复用：ARMA 模型在做最优训练后，可直接返回模型对象给其他程序调用，无需重复训练。

代码实操技巧总结：

• 采用 def 函数式编程封装独立功能模块，并编写规范的 Doc String。
• 使用 prettytable 库（add_row、field_names）实现终端的优雅表格输出。
• 熟练运用 statsmodels 库中的 plot_acf、plot_pacf、adfuller、acorr_ljungbox 等检验工具。
• 掌握基于 AIC/BIC 最小化原则 自动选择 ARMA 模型 p、q 阶数的代码逻辑。
• 灵活使用 lambda 表达式与 datetime.strptime 配合 pandas.read_table 解决非标准时间格式的解析问题。