Python实现SARIMA-LSTM季节性差分自回归移动平均模型结合长短期记忆神经网络时间序列预测（完整源码和数据)

#!/usr/bin/env python # coding: utf-8 # In[1]: import itertools import math import numpy as np import pandas as pd import matplotlib.pyplot as plt from numpy import concatenate from pandas import concat, DataFrame # from statsmodels.tsa.statespace.sarimax import SARIMAX from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score from sklearn.preprocessing import MinMaxScaler from tensorflow.python.keras import Sequential from tensorflow.python.keras.layers import LSTM, Dropout, Dense from statsmodels.tsa.stattools import adfuller from statsmodels.stats.diagnostic import acorr_ljungbox from statsmodels.graphics.tsaplots import plot_acf, plot_pacf import seaborn as sns from statsmodels.graphics.api import qqplot import statsmodels.api as sm # import matplotlib import warnings import statsmodels from scipy import stats # In[2]: warnings.filterwarnings('ignore')#忽略警告 plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签 # In[3]: ''' 首先，通过pandas的read_csv()方法读取名为'焦作.csv'的CSV文件，并选取文件中的第一列和第二列。 接着，使用features列表将数据集中的数据限制在'AQI'这一列。 最后，将限制后的数据集重新赋值给data_raw变量。 该代码的作用是从文件中读取数据并选取特定的列用于分析。 ''' data_raw = pd.read_csv('焦作.csv', usecols=[0, 1]) features=['AQI'] data_raw=data_raw[features] # In[4]: ''' 创建一个大小为(10, 3)的图形对象plt.figure(figsize=(10, 3))。 设置图形标题为'数据AQI'。设置横坐标标签为'time'，纵坐标标签为'AQI'。 使用matplotlib库的plot()方法，绘制data_raw数据集的折线图，将线条颜色设置为蓝色，标签设置为'AQI'。 添加图例legend()。 最后使用show()方法显示图形。 这段代码的目的是可视化数据集中的'AQI'列。 '''# plt.figure(figsize=(10, 3)) plt.title('数据AQI') plt.xlabel('time') plt.ylabel('AQI') plt.plot(data_raw, 'blue', label='AQI') plt.legend() plt.show() # In[5]: #定义稳定性检验函数 def adf_val(ts, ts_title):#定义名为“adf_val”的函数，它有两个参数：时间序列ts和时间序列名称ts_title。 # 稳定性（ADF）检验 adf, pvalue, usedlag, nobs, critical_values, icbest = adfuller(ts)#使用adfuller函数对时间序列ts进行ADF检验，并将返回值分配给adf、pvalue、usedlag、nobs、critical_values和icbest这些变量。 #将ADF检验结果存储在两个列表name和values中，并使用zip函数将它们合并为一个元组的列表，最后打印出来。 name = ['adf', 'pvalue', 'usedlag', 'nobs', 'critical_values', 'icbest'] values = [adf, pvalue, usedlag, nobs,critical_values, icbest] print(list(zip(name, values))) # return adf, pvalue, critical_values, # 返回adf值、adf的p值、三种状态的检验值 # In[6]: #白噪声检验也称为纯随机性检验，当数据是纯随机数据时，再对数据进行分析就没有任何意义了，所以拿到数据后最好对数据进行一个纯随机性检验。 def acorr_val(ts):#定义名为“acorr_val”的函数，它有一个参数：时间序列ts。 lbvalue, pvalue = acorr_ljungbox(ts, lags=1) # 使用acorr_ljungbox函数对时间序列ts进行Ljung-Box检验，并将返回值分配给lbvalue和pvalue这两个变量。 # return lbvalue, pvalue#返回lbvalue和pvalue这两个变量的值。 # In[7]: def tsplot(y, lags=None, figsize=(14, 8)):#定义名为“tsplot”的函数，它有三个参数：时间序列y、滞后值lags和图形大小figsize。 #创建一个图形对象，并将其大小设置为figsize。然后，定义一个2x2的布局，并将时间序列y、直方图、自相关图和偏自相关图分别放置在四个子图中。 fig = plt.figure(figsize=figsize) layout = (2, 2) ts_ax = plt.subplot2grid(layout, (0, 0)) hist_ax = plt.subplot2grid(layout, (0, 1)) acf_ax = plt.subplot2grid(layout, (1, 0)) pacf_ax = plt.subplot2grid(layout, (1, 1)) #在时间序列子图中绘制时间序列y，并将子图的标题设置为“price”。 y.plot(ax=ts_ax) ts_ax.set_title('price') #在直方图子图中绘制时间序列y的直方图，并将子图的标题设置为“Histogram”。 y.plot(ax=hist_ax, kind='hist', bins=25) hist_ax.set_title('Histogram') #自相关： 对一个时间序列，现在值与其过去值的相关性。如果相关性为正，则说明现有趋势将继续保持。 plot_acf(y, lags=lags, ax=acf_ax)#在自相关图子图中绘制时间序列y的自相关图，如果有指定滞后值，则使用lags。 #可以度量现在值与过去值更纯正的相关性 plot_pacf(y, lags=lags, ax=pacf_ax)#在偏自相关图子图中绘制时间序列y的偏自相关图，如果有指定滞后值，则使用lags #将自相关图和偏自相关图的x轴范围限制为非负值，并去掉图形的上、右边框。然后，调整子图的布局，使它们不重叠。最后，显示图形。 [ax.set_xlim(0) for ax in [acf_ax, pacf_ax]] sns.despine() fig.tight_layout() fig.show() # return ts_ax, acf_ax, pacf_ax#返回时间序列子图、自相关图子图和偏自相关图子图。 # In[8]: data = data_raw.astype('float32')#将data_raw转换为float32数据类型，赋值给变量data。 # In[9]: #adf结果如果小于三个level的统计值。pvalue也是接近于0 的，那么是平稳的 adf, pvalue1, critical_values = adf_val(data, 'raw time series')#对数据框df进行ADF单位根检验，检验结果赋值给变量ADF、pvalue1和critical_values print('adf',adf) print('pvalue1',pvalue1) print('critical_values',critical_values) #若p值远小于0.01，认为该时间序列是平稳的 aco=acorr_val(data)#对data进行自相关性检验，检验结果赋值给变量aco print('aco',aco) # In[10]: diff1 = data.diff(1)#对data数据框进行一阶差分操作，即每个时间点的值减去前一个时间点的值，得到的差分结果赋值给变量diff1。 #画出一阶差分后的数据 plt.figure(figsize=(12,3)) plt.plot(diff1, 'blue', label='diff1') plt.title('diff1') plt.xlabel('time') plt.ylabel('price') plt.legend() plt.show() # In[11]: ''' 将diff1中的缺失值（NaN）用0进行填充，得到填充后的差分结果。 然后，将差分结果中的正无穷和负无穷值（Inf）用0进行替换，得到最终处理后的差分结果。 这两行代码的目的是处理差分结果中可能出现的异常值。 '''# diff1[np.isnan(diff1)] = 0 diff1[np.isinf(diff1)] = 0 # In[12]: #对差分后的数据再次进行检验 adf, pvalue1, critical_values = adf_val(diff1, 'raw time series')# print('adf',adf) print('pvalue1',pvalue1) print('critical_values',critical_values) aco=acorr_val(diff1)# print('aco',aco) # In[13]: ##画出自相关、偏自相关、直方图 tsplot(data, lags=20)# # In[14]: ''' 定义了3个range对象，分别表示ARIMA模型中p、d和q的取值范围。 使用itertools.product()方法生成pdq元组列表，其中包含了从p、d和q的所有可能组合。 seasonal_pdq是一个元组列表，每个元组中包含了p、d、q和季节周期的值，这些值用于季节性ARIMA模型中。 该代码的作用是生成ARIMA模型中各参数的取值范围。 ''' # p = range(0, 3) d = range(0, 1) q = range(0, 3) pdq = list(itertools.product(p, d, q)) seasonal_pdq = [(x[0], x[1], x[2], 6) for x in list(itertools.product(p, d, q))] # In[15]: ''' 使用两个嵌套循环遍历pdq和seasonal_pdq中所有可能的参数组合，并构建SARIMA模型。 在每次循环中，使用SARIMAX()方法构建SARIMA模型，并使用模型拟合数据，获得模型拟合结果。 使用format()方法将模型的AIC值输出到控制台上，方便比较每个模型的优劣。 使用if语句更新AIC最小的模型。 如果模型构�