ARIMA PYTHON 计算RMSE python arima预测

ARIMA模型预测

时间序列分析预测就是在已有的和时间有关的数据序列的基础上构建其数据模型并预测其未来的数据，例如航空公司的一年内每日乘客数量、某个地区的人流量，这些数据往往具有周期性的规律。如下图所示，有的数据呈现出简单的周期性循环，有的呈现出周期性循环变化。

   

ARIMA（Autoregressive Integrated Moving Average model，差分整合移动平均自回归模型）不仅可以很好地对已有数据进行拟合，并且还可以在此基础上对未来数据进行预测。其函数原型为：ARIMA(data,order=(p,d,q))，因此除了需要传入原始数据data之外，还需要三个参数p、d、q，这三个参数与ARIMA的原理有关。

ARIMA可以拆分为AR、I、MA，分别对应p、d、q三个参数。

1. AR（AutoRegressive，自回归模型）描述当前值与历史值之间的关系，用变量自身的历史时间数据对自身进行预测，自回归过程的公式定义： yt是当前值 u是常数项 P是阶数 ri是自相关系数 et是误差。因此参数p为自回归模型的阶数。
2. MA（Moving Average，移动平均值）利用一定时间间隔内的平均值作为某一期的估计值，这样可以消除数据的周期性影响。其公式为。因此参数q为移动平均的阶数。
3. I代表差分操作，时间序列最常用来剔除周期性因素的方法，它主要是对等周期间隔的数据进行线性求减。从而使数据变得平稳，ARIMA一般进行一次差分即可稳定，因此d一般取值为0、1、2，这里取1，进行一次差分。

确定参数值

根据BIC

因此要使用ARIMA模型进行预测首先需要确定参数p、q的值。因为一般阶数不超过整体数据的十分之一，因此这里采用穷举法。分别从0~10取p、q，并且得到该模型的BIC值（贝叶斯信息量 bayesian information criterion）。BIC准则综合考虑了残差大小和自变量的个数，残差越小、自变量个数越少BIC值越小，模型越优。因此从所有的p、q取值中找出bic值最小的即为理想的模型参数。经过如下函数，我得到p、q为5、5

def get_order(data):
    pmax = int(len(data) / 10)    #一般阶数不超过 length /10
    qmax = int(len(data) / 10)
    bic_matrix = []
    for p in range(pmax +1):
        temp= []
        for q in range(qmax+1):
            try:
                temp.append(ARIMA(data, order=(p, 1, q)).fit().bic)    # 将bic值存入二维数组
            except:
                temp.append(None)
        bic_matrix.append(temp)    
    bic_matrix = pd.DataFrame(bic_matrix)   #将其转换成Dataframe 数据结构
    p,q = bic_matrix.stack().astype('float64').idxmin()        # 找出bic值最小对应的索引
    return p,q

根据图像

还可以根据自相关系数ACF与偏自相关系数PACF图像来确定阶数p、q的值，若下所示使用库函数输出数据的图像

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

def draw_acf_pacf(data,lags):
    f = plt.figure(facecolor='white')
    ax1 = f.add_subplot(211)
    plot_acf(data,ax=ax1,lags=lags)
    ax2 = f.add_subplot(212)
    plot_pacf(data,ax=ax2,lags=lags)
    plt.subplots_adjust(hspace=0.5)
    plt.show()

绘制图像有如下两种类型，拖尾是指序列以指数率单调递减或震荡衰减，如下左图所示。截尾指序列从某个时点变得非常小，如下右图所示。

 

根据ACF图像与PACF图像的拖尾结尾选取不同的模型：

• 自相关系数拖尾，则q=0，偏自相关系数p阶截尾：AR(p)模型，
• 自相关系数q阶截尾，偏自相关函数拖尾，p=0：MA(q)模型
• 自相关函数和偏自相关函数均拖尾，阶数分别为q、p：ARMA(p,q)模型

 那么图像如何确定阶数呢？即观察图片从第几个数据之后收敛到置信区间，例如上面右图从24之后再无超出置信区间的数据，所以认为其24阶截尾。

进行预测

接下来就可以使用arima模型进行模型拟合与预测了，这里使用的是python第三方包statsmodels.tsa.arima.model中的ARIMA模型。这是Statsmodels自从0.11版本新独立的模块，其原来的模块为statsmodels.tsa.arima_model.ARIMA，二者在功能上都实现了arima模型，并且具有相同的属性和方法名，其返回值均为ARIMAResults对象，通过该对象的predict()、forecast()得到预测值。值得注意的是新旧模块的方法的传入参数和返回值类型不太相同，因此使用时需要注意。其官方文档连接为：statsmodels.tsa.arima.model.ARIMAResults、statsmodels.tsa.arima_model.ARIMAResults

这里使用predict()函数进行预测，其可以传入start、end参数代表预测数据的开始和结束坐标，坐标值不仅可以是int、str，还可以是 datetime时间类型，如果start、end介于原有数据区域内，即为对原数据的预测拟合，如果end超过了原有数据长度，即代表对未来数据进行预测。forecast()函数接收一个steps参数，代表对未来多少个数据进行预测，也可以传入datetime时间，代表从样本数据结束一直预测到该时间。特别注意的是在使用pandas的Datetime时间索引的数据进行训练和预测时，要保证时间间隔是有规律的，否则在预测时索引会出错而报错如下：

KeyError: 'The `end` argument could not be matched to a location related to the index of the data.'

import h5py
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from statsmodels.tsa.arima.model import ARIMA

%matplotlib inline  
plt.rcParams['font.sans-serif'] = ['SimHei']    # 设置字符集显示中文
plt.rcParams['axes.unicode_minus'] = False      # 设置负号正确显示

# 显示数据趋势图
def draw_data(data):
    plt.figure()
    plt.plot(data)
    
with h5py.File('./taxi.h5', 'r') as hf:
    # 读取数据
    in_data=np.array(hf['in_data'])
    in_area = in_data.transpose(2, 0, 1)  # 调整维度
    in_mean_day = in_area.mean(axis=2)
    sub_data=pd.Series(in_mean_day[116][0:100])      # 截取数据集中的100个作为样本
    draw_data(sub_data)
    
    # 确定模型阶数
    # p,q=get_order(sub_data)
    # print('阶数：',p,q) 
    
    # 利用ARIMA模型进行预测
    model=ARIMA(sub_data,order=(5,1,5)).fit()        # 传入参数，构建并拟合模型
    predict_data=model.predict(0,150)                # 预测数据
    forecast=model.forecast(21)                      # 预测未来数据
    
    # 绘制原数据和预测数据对比图
    plt.plot(sub_data,label='原数据')
    plt.plot(predict_data,label='预测数据')
    plt.legend()
    plt.show()

如下所示为拟合预测结果：

 数据预处理

稳定性

ARIMA所需要的数据是稳定的，这样才能进行规律的发掘和拟合，稳定的数据是指其均值和方差为一个常量、自协方差与时间无关，如下所示分别为稳定的数据和均值不稳定、方差不稳定、自协方差随时间变化

 

数据的稳定性可以利用单位根检验Dickey-Fuller Test 进行判断，即在一定置信水平下，假设时序数据Null hypothesis: 非稳定、序列具有单位根。因此对于一个平稳的时序数据，就需要在给定的置信水平上显著，拒绝原假设。如下通过statsmodels包中的adfuller来进行检测：

from statsmodels.tsa.stattools import adfuller

# 利用ADF检测数据是否平稳
def check_stable(data):
    adf_res=adfuller(data)
    print('t统计量',adf_res[0])
    print('t统计量的P值',adf_res[1])
    print('延迟阶数',adf_res[2])
    print('观测值的个数',adf_res[3])
    for key,value in adf_res[4].items():
        print('临界区间：',key,'，值：',value)

 统计结果如下：其中P值低于0.05，可以拒绝原假设，即数据是平稳的

 移动平均

那么如何使数据变得平稳，首先可以使用数据平滑技术消除数据的周期性起伏变化，常用的平滑技术有移动平均、加权移动平均，移动平均即利用一定时间间隔内的平均值作为某一期的估计值，而指数平均则是用变权的方法来计算均值。通过pandas的rolling()方法可以实现数据移动，ewa()实现加权移动，mean()实现平均。step平均的间隔。

# 对数据进行移动平均来平滑数据
def move_average(data,step):
    series=pd.Series(data)
    rol_mean=series.rolling(window=step).mean()     # 移动平均
    rol_weight_mean=pd.DataFrame.ewm(series,span=step).mean()       # 加权移动平均
    
    return rol_mean,rol_weight_mean

差分

通过将相同间隔的数据做差来消除周期性的因素。pandas的diff可以完成对序列相隔steps作差，例如当steps=7时，前7个数据是没法和前面的数据作差的，因此前七个数据为Nan，通过dropna()删除。

在预测操作之后，还需要将差分的数据还原回来。差分操作只是简单的相减，因此还原就是相加，将原数据通过shift向后移动7个位置然后相加，即可实现还原操作。

# 进行差分
diff=series.diff(steps)
diff=diff.dropna()

# 差分还原
diff_shift=sub_data.shift(7)
diff_recover=predict_data.add(diff_shift)