参考资料1
https://pyflux.readthedocs.io/en/latest/arima.html#example
参考资料2
多元序列分析ARIMAX(p,I,q)
没事瞎折腾
4 人赞同了该文章
这里借助Python的statsmodels库和pyflux库进行多元时间序列分析，建立ARIMAX(p,I,q)模型用来预测二氧化碳浓度数据。其中pyflux库是一个专门用来建立时间序列模型的python库。该库的文档地址为：
https://pyflux.readthedocs.io/en/latest/getting_started.htmlpyflux.readthedocs.io
文章中完整的jupyter notebook程序可在下面链接查看。
Jupyter Notebook Viewernbviewer.jupyter.org
数据集来自：
https://www.itl.nist.gov/div898/handbook/datasets/GAS_FURNACE.DATwww.itl.nist.gov
1:数据准备和可视化
## 加载包
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
%config InlineBackend.figure_format = "retina"
from matplotlib.font_manager import FontProperties
fonts = FontProperties(fname = "/Library/Fonts/华文细黑.ttf")
import statsmodels.api as sm
import pyflux as pf
from sklearn.metrics import mean_absolute_error,mean_squared_error
## 读取数据,数据来自https://www.itl.nist.gov/div898/handbook/datasets/GAS_FURNACE.DAT
datadf = pd.read_csv(".../data/gas furnace data.txt",sep="\s+")
datadf.columns = ["GasRate","C02"]
## GasRate:输入天然气速率，C02：输出二氧化碳浓度
datadf.head()
接下来可视化两列数据
## 可视化两列数据
plt.figure(figsize=(15,5))
plt.subplot(1,2,1)
datadf.GasRate.plot(c="r")
plt.xlabel("Observation")
plt.ylabel("Gas Rate")
plt.subplot(1,2,2)
datadf.C02.plot(c="r")
plt.xlabel("Observation")
plt.ylabel("C02")
plt.show()
切分数据集，前面百分之75做训练集，后面百分之25做测试集
## 前面百分之75做训练集，后面百分之25做测试集
trainnum = np.int(datadf.shape[0]*0.75)
traidata = datadf.iloc[0:trainnum,:]
testdata = datadf.iloc[trainnum:datadf.shape[0],:]
print(traidata.shape)
print(testdata.shape)
(222, 2)
(74, 2)2：平稳时间序列建模ARIMAX
2.1:单位根检验检验序列的平稳性
因为ARIMAX(p,i,q)要求所有的序列的是平稳的，所以要对序列进行单位根检验，判断序列的平稳性。
## 1:单位根检验检验序列的平稳性,ADF 检验
dftest = sm.tsa.adfuller(datadf.GasRate,autolag='BIC')
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','Lags Used','Number of Observations Used'])
print("GasRate 检验结果：")
print(dfoutput)
dftest = sm.tsa.adfuller(datadf.C02,autolag='BIC')
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','Lags Used','Number of Observations Used'])
print("C02 检验结果：")
print(dfoutput)
GasRate 检验结果：
Test Statistic -4.878952
p-value 0.000038
Lags Used 2.000000
Number of Observations Used 293.000000
dtype: float64
C02 检验结果：
Test Statistic -2.947057
p-value 0.040143
Lags Used 3.000000
Number of Observations Used 292.000000
dtype: float64p-value均小于0.05，说明在置信度为95％水平下，两个序列均是平稳序列。
2.2:接下来可以可视化C02数据的自相关系数和偏相关系数。
## 可视化序列的自相关和偏相关图
fig = plt.figure(figsize=(10,5))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(traidata.C02, lags=30, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(traidata.C02, lags=30, ax=ax2)
plt.subplots_adjust(hspace = 0.3)
plt.show()
2.3:接下来建立一个简单的ARIMAX模型，探索库中函数的相关使用方法。
## 建立ARIMAX（1，0，2）模型
model = pf.ARIMAX(data=traidata,formula="C02~GasRate",ar=1,ma=2,integ=0)
model_1 = model.fit("MLE")
model_1.summary()
Normal ARIMAX(1,0,2)
======================================================= ==================================================
Dependent Variable: C02 Method: MLE
Start Date: 2 Log Likelihood: -71.9362
End Date: 221 AIC: 155.8725
Number of observations: 220 BIC: 176.2343
==========================================================================================================
Latent Variable Estimate Std Error z P>|z| 95% C.I.
======================================== ========== ========== ======== ======== =========================
AR(1) 0.9086 0.0191 47.5425 0.0 (0.8712 | 0.9461)
MA(1) 1.0231 0.0552 18.5272 0.0 (0.9149 | 1.1314)
MA(2) 0.6231 0.0442 14.1127 0.0 (0.5365 | 0.7096)
Beta 1 4.8793 1.0166 4.7996 0.0 (2.8868 | 6.8719)
Beta GasRate -0.4057 0.0533 -7.613 0.0 (-0.5102 | -0.3013)
Normal Scale 0.3356
==========================================================================================================
model_1.summary()方法可以输出模型的拟合的相关信息。
## 可视化模型在训练集上的拟合情况
model.plot_fit(figsize=(10,5))
model.plot_fit()方法可以将在训练数据中模型的拟合情况进行可视化。可以发现在训练集上模型的拟合效果很好。
## 可视化模型的在测试集上的预测结果
model.plot_predict(h=testdata.shape[0], ## 往后预测多少步
oos_data=testdata, ## 测试数据集
past_values=traidata.shape[0], ## 图像显示训练集的多少数据
figsize=(15,5))
## 可视化原始数据
datadf.C02.plot(figsize=(15,5))
plt.xlabel("Time")
plt.ylabel("C02")
plt.show()model.plot_predict()方法可以将预测结果可视化
使用model.predict()方法可以用来预测新的数据集
## 预测新的数据
C02pre = model.predict(h=testdata.shape[0], ## 往后预测多少步
oos_data=testdata, ## 测试数据集
intervals=True, ## 同时预测置信区间
)
print("在测试集上mean absolute error:",mean_absolute_error(testdata.C02,C02pre.C02))
print("在测试集上mean squared error:",mean_squared_error(testdata.C02,C02pre.C02))
C02pre.head()
在测试集上mean absolute error: 1.5731456243696424
在测试集上mean squared error: 3.8376299478820215
## 可视化原始数据和预测数据进行对比
datadf.C02.plot(figsize=(15,5),c="b",label="C02")
C02pre.C02.plot(c = "r",label="Prediction")
plt.xlabel("Time")
plt.ylabel("C02")
plt.legend(loc=0)
plt.show()
从图像上可以看出预测的效果和实际数据的差异。
3.通过遍历寻找合适的P，Q
p = np.arange(6)
q = np.arange(6)
pp,qq = np.meshgrid(p,q)
resultdf = pd.DataFrame(data = {"arp":pp.flatten(),"mrq":qq.flatten()})
resultdf["bic"] = np.double(pp.flatten())
resultdf["mae"] = np.double(qq.flatten())
## 迭代循环建立多个模型
for ii in resultdf.index:
model_i = pf.ARIMAX(data=traidata,formula="C02~GasRate",ar=resultdf.arp[ii],ma=resultdf.mrq[ii],integ=0)
try:
modeli_fit = model_i.fit("MLE")
bic = modeli_fit.bic
C02_pre = model.predict(h=testdata.shape[0],oos_data=testdata)
mae = mean_absolute_error(testdata.C02,C02_pre.C02)
except:
bic = np.nan
resultdf.bic[ii] = bic
resultdf.mae[ii] = mae
print("模型迭代结束")
模型迭代结束
## 按照BIC寻找合适的模型
resultdf.sort_values(by="bic").head()
## 重新建立效果较好的模型
model = pf.ARIMAX(data=traidata,formula="C02~GasRate",ar=4,ma=1,integ=0)
model_1 = model.fit("MLE")
model_1.summary()
Normal ARIMAX(4,0,1)
======================================================= ==================================================
Dependent Variable: C02 Method: MLE
Start Date: 4 Log Likelihood: 21.1205
End Date: 221 AIC: -26.2409
Number of observations: 218 BIC: 0.835
==========================================================================================================
Latent Variable Estimate Std Error z P>|z| 95% C.I.
======================================== ========== ========== ======== ======== =========================
AR(1) 2.4472 0.0647 37.8416 0.0 (2.3205 | 2.574)
AR(2) -2.1984 0.1625 -13.5296 0.0 (-2.5169 | -1.8799)
AR(3) 0.7821 0.149 5.2483 0.0 (0.49 | 1.0741)
AR(4) -0.0596 0.0496 -1.2008 0.2298 (-0.1569 | 0.0377)
MA(1) -0.9278 0.0268 -34.6168 0.0 (-0.9803 | -0.8752)
Beta 1 1.5274 0.1027 14.869 0.0 (1.326 | 1.7287)
Beta GasRate -0.1015 0.0071 -14.2897 0.0 (-0.1154 | -0.0876)
Normal Scale 0.2189
==========================================================================================================
## 可视化潜在变量
model.plot_z()model.plot_z()可以对潜在变量进行可视化
预测新的数据并可视化
## 预测新的数据
C02pre = model.predict(h=testdata.shape[0], ## 往后预测多少步
oos_data=testdata, ## 测试数据集
)
print("在测试集上mean absolute error:",mean_absolute_error(testdata.C02,C02pre.C02))
## 可视化原始数据和预测数据进行对比
datadf.C02.plot(figsize=(15,5),c="b",label="C02")
C02pre.C02.plot(c = "r",label="Prediction")
plt.xlabel("Time")
plt.ylabel("C02")
plt.legend(loc=0)
plt.show()
扫一扫
1695
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
