时间序列预测-ARIMA

参考资料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.html​pyflux.readthedocs.io

 

文章中完整的jupyter notebook程序可在下面链接查看。

Jupyter Notebook Viewer​nbviewer.jupyter.org图标

数据集来自:

https://www.itl.nist.gov/div898/handbook/datasets/GAS_FURNACE.DAT​www.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: float64

p-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()

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值