Python自动确定SARIMA模型参数

简介

ARIMA (AutoRegressive Integrated Moving Average) 模型是一种广泛用于时间序列数据分析和预测的统计模型。

ARIMA 模型结合了自回归 (AR) 和移动平均 (MA) 模型，并且可以处理非平稳时间序列数据，通过差分操作来使其平稳。ARIMA 模型的三个主要参数是 ( p ), ( d ), 和 ( q )。

• ( p ): 自回归 (AR) 项数，即模型使用的过去观测值的数量。
• ( d ): 差分次数，即将时间序列数据差分多少次以使其平稳。
• ( q ): 移动平均 (MA) 项数，即模型使用的过去预测误差的数量。

ARIMA 模型的构建步骤

1. 确定时间序列的平稳性:

• 绘制时间序列图，观察数据的趋势和季节性。
• 使用单位根检验（如ADF检验）来检测时间序列的平稳性。
• 如果数据非平稳，可以进行差分处理，使其平稳。

2. 模型识别与选择:

• 使用自相关函数 (ACF) 和偏自相关函数 (PACF) 图来识别合适的 ( p ) 和 ( q ) 值。
• ACF 图用于识别 MA 项数 ( q )，PACF 图用于识别 AR 项数 ( p )。

3. 参数估计:

• 使用最小二乘法或最大似然估计法来估计模型参数。

4. 模型诊断:

• 检查模型的残差，确保残差是白噪声（即没有显著的自相关）。
• 使用Ljung-Box检验或其他统计检验来验证残差的独立性。

5. 模型预测:

• 使用拟合好的 ARIMA 模型对未来的时间点进行预测。

ARIMA 模型的应用示例 (Python)

下面是一个简单的 Python 示例，使用 statsmodels 库构建和使用 ARIMA 模型：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# 读取时间序列数据
data = pd.read_csv('time_series_data.csv', index_col='Date', parse_dates=True)
ts = data['Value']

# 绘制时间序列图
plt.plot(ts)
plt.title('Time Series Data')
plt.show()

# 检验平稳性
result = adfuller(ts)
print('ADF Statistic:', result[0])
print('p-value:', result[1])

# 如果数据非平稳，进行差分处理
if result[1] > 0.05:
    ts_diff = ts.diff().dropna()

# 绘制 ACF 和 PACF 图
plot_acf(ts_diff)
plot_pacf(ts_diff)
plt.show()

# 构建 ARIMA 模型
p = 1  # 根据 PACF 图选择
d = 1  # 差分一次
q = 1  # 根据 ACF 图选择

model = ARIMA(ts, order=(p, d, q))
model_fit = model.fit()

# 打印模型摘要
print(model_fit.summary())

# 进行预测
forecast = model_fit.forecast(steps=10)
print(forecast)

ARIMA时间序列预测模型Python代码

import pandas 
# 读取数据，指定日期为索引列
data = pandas.read_csv(
    'D:\\DATA\\pycase\\number2\\9.3\\Data.csv' ,
    index_col='日期'
)
# 绘图过程中
import  matplotlib.pyplot as plt
# 用来正常显示中文标签
plt.rcParams['font.sans-serif']=['SimHei']
# 用来正常显示负号
plt.rcParams['axes.unicode_minus'] = False 
# 查看趋势图
data.plot() #有增长趋势，不平稳
# 附加：查看自相关系数合片自相关系数（查分之后），可以用于平稳性的检测，也可用于定阶系数预估
#自相关图（）
from statsmodels.graphics.tsaplots import plot_acf
plot_acf(data).show() #自相关图既不是拖尾也不是截尾。以上的图的自相关是一个三角对称的形式，这种趋势是单调趋势的典型图形，说明这个序列不是平稳序列

# 1 平稳性检测
from statsmodels.tsa.stattools import adfuller as ADF
def tagADF(t):
    result = pandas.DataFrame(index=[
            "Test Statistic Value", "p-value", "Lags Used", 
            "Number of Observations Used", 
            "Critical Value(1%)", "Critical Value(5%)", "Critical Value(10%)"
        ], columns=['销量']
    );
    result['销量']['Test Statistic Value'] = t[0]
    result['销量']['p-value'] = t[1]
    result['销量']['Lags Used'] = t[2]
    result['销量']['Number of Observations Used'] = t[3]
    result['销量']['Critical Value(1%)'] = t[4]['1%']
    result['销量']['Critical Value(5%)'] = t[4]['5%']
    result['销量']['Critical Value(10%)'] = t[4]['10%']
    return result;
print('原始序列的ADF检验结果为:',tagADF(ADF(data[u'销量'])))  # 添加标签后展现
 
# 平稳判断：得到统计量大于三个置信度(1%,5%,10%)临界统计值，p值显著大于0.05，该序列为非平稳序列。
# 备注：得到的统计量显著小于3个置信度（1%，5%，10%）的临界统计值时，为平稳 此时p值接近于0 此处不为0，尝试增加数据量，原数据太少

# 2 进行数据差分，一般一阶差分就可以
D_data = data.diff(1).dropna()
D_data.columns = [u'销量差分']
#差分图趋势查看
D_data.plot() 
plt.show()
# 附加：查看自相关系数合片自相关系数（查分之后），可以用于平稳性的检测，也可用于定阶系数预估
#自相关图
plot_acf(D_data).show()
plt.show()
#偏自相关图
from statsmodels.graphics.tsaplots import plot_pacf
plot_pacf(D_data).show()
 
# 3 平稳性检测
print(u'差分序列的ADF检验结果为：', tagADF(ADF(D_data[u'销量差分']))) 
# 解释：Test Statistic Value值小于两个水平值，p值显著小于0.05，一阶差分后序列为平稳序列。
 
# 4 白噪声检验
from statsmodels.stats.diagnostic import acorr_ljungbox
#返回统计量和p值
print(u'差分序列的白噪声检验结果为：', acorr_ljungbox(D_data, lags=1))  # 分别为stat值（统计量）和P值
# P值小于0.05，所以一阶差分后的序列为平稳非白噪声序列。
 
# 5 p，q定阶
from statsmodels.tsa.arima_model import ARIMA
#一般阶数不超过length/10
pmax = int(len(D_data)/10) 
#一般阶数不超过length/10
qmax = int(len(D_data)/10) 
#bic矩阵
bic_matrix = [] 
for p in range(pmax+1):
  tmp = []
  for q in range(qmax+1):
#存在部分报错，所以用try来跳过报错。
    try: 
      tmp.append(ARIMA(data, (p,1,q)).fit().bic)
    except:
      tmp.append(None)
  bic_matrix.append(tmp)
#从中可以找出最小值
bic_matrix = pandas.DataFrame(bic_matrix) 
#先用stack展平，然后用idxmin找出最小值位置。
p,q = bic_matrix.stack().idxmin() 
print(u'BIC最小的p值和q值为：%s、%s' %(p,q))
# 取BIC信息量达到最小的模型阶数，结果p为0，q为1，定阶完成。
 
# 6 建立模型和预测
model = ARIMA(data, (p,1,q)).fit() 
#给出一份模型报告
model.summary2() 
#作为期5天的预测，返回预测结果、标准误差、置信区间。
model.forecast(5)