时间序列模型（ARIMA和ARMA）完整步骤详述_arma模型建模步骤

我于2019年发布此篇文章至今收获了许多人的指点，当时的代码的确晦涩难懂，近期有空，将代码重新整理了一遍，重新发送至此。希望能够帮助大家更好地理解。

建模步骤：

目录

数据包和版本申明

步骤一：数据准备与数据预处理

步骤二：数据重采样

步骤三：平滑处理

步骤四：平稳性检验

 步骤五： 时间序列定阶

（2）信息准则定阶

步骤六：模型构建

步骤七：模型评价

总结

---

数据包和版本申明

申明：本实验环境为python 3.7.4  statsmodels版文为：0.10.1

import pandas as pd 
import numpy as np
import seaborn as sns #热力图
import itertools 
import datetime
import matplotlib.pyplot as plt #画图
import statsmodels.api as sm 
from sklearn.metrics import mean_squared_error
from statsmodels.tsa.stattools import adfuller #ADF检验
from statsmodels.stats.diagnostic import acorr_ljungbox #白噪声检验
from statsmodels.graphics.tsaplots import plot_acf,plot_pacf #画图定阶
from statsmodels.tsa.arima_model import ARIMA #ARIMA模型
from statsmodels.tsa.arima_model import ARMA #ARMA模型
from statsmodels.stats.stattools import durbin_watson #DW检验
from statsmodels.graphics.api import qqplot #qq图

步骤一：数据准备与数据预处理

    自动生成2018年1月1日至2018年9月1日数据，数据服从标准正态分布，存入old_data.csv中。

#### Part:generate raw data and save in the old_data.csv 
#### 创建一个时间列表，从20180101到20180901数据，存入 old_data.csv 
def genertate_data():
    index = pd.date_range(start='2018-1-1',end = '2018-9-1',freq='10T') # 10分钟采样一次
    index = list(index)
    data_list = []
    for i in range(len(index)): 
        data_list.append(np.random.randn())  # 数据是符合标准正态分布的样本
    dataframe = pd.DataFrame({'time':index,'values':data_list})
    dataframe.to_csv('G:\\WX\\2\\old_data.csv',index=0)
    print('the data is existting')

故意去将文件夹中的某些值，改成了-10000，弄成了异常值，（因为老师说尽可能显得步骤完整，最后分数才会高-，-所以我自己手动添加异常）。这块的主要工作就是利用pandas里面的函数，去查看一下刚特殊操作后的数据。

#### Step 1 数据预处理
#### delete or revise some values in data and make data preprocessing
#### 删掉或者修改创建的数据后，进行简单数据预处理 
def data_preprocessing():
    data = pd.read_csv('G:\\WX\\2\\old_data.csv')
    #print(data.describe()) #查看统计信息,发现最小值有-10000的异常数据
    #print((data.isnull()).sum()) #查看是否存在缺失值
    #print((data.duplicated()).sum()) #重复值
    def change_zero(x):
        if x == -10000:
            return 0
        else :
            return x
    data['values'] = data['values'].apply(lambda x: change_zero(x))
 
    #利用均值填充缺失值
    mean = data['values'].mean()
    def change_mean(x):
        if x == 0:
            return mean
        else:
            return x
    data['values'] = data['values'].apply(lambda x: change_mean(x))
    #保存处理过的数据
    data.to_csv('G:\\WX\\2\\new_data.csv',index=0)
    print('new data is existing')

步骤二：数据重采样

为了得高分（-，-），做了很多个数据，然后一共有34992个数据，然后进行了一下重采样，数据以天进行重采样。

#### Step 2 重采样
#### Resample Data and Sampling frequency is days
#### 重采样，将采样频率换成以天为单位
def Resampling(): #重采样
    df = pd.read_csv('G:\\WX\\2\\new_data.csv')
     #将默认索引方式转换成时间索引
    df['time'] = pd.to_datetime(df['time'])
    df.set_index("time", inplace=True)
   
    train_data = df['2018-1-1':'2018-8-1'] ## 取到20180101 至 20180801 做训练 
    test = df['2018-8-1':'2018-9-1']       ## 取到20180801 至 20180901 做预测 
    train_data = train_data.resample('D').mean()  ## 以天为时间间隔取均值,重采样
    test_data = test.resample('D').mean()
 
    return train_data,test_data

步骤三：平滑处理

由于ARMA和ARIMA需要时间序列满足平稳性和非白噪声的要求，所以要用差分法和平滑法（滚动平均和滚动标准差）来实现序列的平稳性操作。一般情况下，对时间序列进行一阶差分法就可以实现序列的平稳性，有时需要二阶差分。

#### Step 3  差分转平稳
def stationarity(timeseries): #平稳性处理（timeseries 时间序列）
    ## 差分法,保存成新的列
    diff1 = timeseries.diff(1).dropna()  # 1阶差分 dropna() 删除缺失值
    diff2 = diff1.diff(1) #在一阶差分基础上再做一次一阶差分，即二阶查分
    ## 画图
    #diff1.plot(color = 'red',title='diff 1',figsize=(10,4))
    #diff2.plot(color = 'black',title='diff 2',figsize=(10,4))
 
    
    ## 平滑法
    rollmean = timeseries.rolling(window=4,center = False).mean() ## 滚动平均
    rollstd = timeseries.rolling(window=4,center = False).std() ## 滚动标准差
    ## 画图 
    #rollmean.plot(color = 'yellow',title='Rolling Mean',figsize=(10,4))
    #rollstd.plot(color = 'blue',title='Rolling Std',figsize=(10,4))
    
    return diff1,diff2,rollmean,rollstd

差分法处理结果图：由图可以看出 一阶差分基本就满足了平稳性需要。

平滑法处理结果如图所示。

可以看出，平滑法不太适合我造出来的数据。一般情况下，平滑法更适合带有周期性稳步上升的数据类型。

步骤四：平稳性检验

利用ADF检验判断序列是否平稳，利用白噪声检验判断序列是否为随机性序列。

#### Step 4  平稳性检验
def ADF_test(timeseries): ## 用于检测序列是否平稳
    x = np.array(timeseries['values'])
    adftest = adfuller(x, autolag='AIC')
    #print (adftest) 
    if adftest[0] < adftest[4]["1%"] and adftest[1] < 10**(-8): 
    # 对比Adf结果和10%的时的假设检验 以及 P-value是否非常接近0(越小越好)
        print("序列平稳")
        return True 
    else:
        print("非平稳序列")
        return False
 
def random_test(timeseries) : #随机性检验（白噪声检验）
    p_value = acorr_ljungbox(timeseries, lags=1)  # p_value 返回二维数组，第二维为P值
    if p_value[1] < 0.05: 
        print("非随机性序列")
        return  True
    else:
        print("随机性序列,即白噪声序列")
        return False

（1）ADF检验结果如下：

如何确定该序列能否平稳呢？主要看：

（1）1%、%5、%10不同程度拒绝原假设的统计值和ADF Test result的比较，ADF Test result同时小于1%、5%、10%即说明非常好地拒绝该假设，本数据中，adf结果为-6.9， 小于三个level的统计值。

（2）P-value是否非常接近0.本数据中，P-value 为 7.9e-10，接近0。

ADF结果如何查看参考了这篇博客：

Python时间序列中ADF检验详解_学渣渣-CSDN博客_python进行adf检验

（2）白噪声结果如图：

                                                        

统计量的P值小于显著性水平0.05，则可以以95%的置信水平拒绝原假设，认为序列为非白噪声序列（否则，接受原假设，认为序列为纯随机序列。）         

由于P值为0.315远大于0.05所以接受原假设，认为时间序列是白噪声的，即是随机产生的序列，不具有时间上的相关性。（解释一下，由于老师没有给数据，所以只能硬着头皮，假设它是非白噪声的做）

 步骤五： 时间序列定阶

定阶方法主要为两种：

（1）ACF和PACF 利用拖尾和截尾来确定

（2）信息准则定阶（AIC、BIC、HQIC）

热力图为辅助方法使得p与q的取值更加明确,可视化,实际上依旧为方法二-信息准则定阶。

def determinate_order_acf(timeseries): #利用ACF和PACF判断模型阶数 
    plot_acf(timeseries,lags=40) #延迟数
    plot_pacf(timeseries,lags=40)
    plt.show()
     
def detetminante_order_AIC(timeseries): #信息准则定阶：AIC、BIC、HQIC
    #AIC
    AIC = sm.tsa.arma_order_select_ic(timeseries,\
        max_ar=6,max_ma=4,ic='aic')['aic_min_order']
    #BIC
    BIC = sm.tsa.arma_order_select_ic(timeseries,max_ar=6,\
           max_ma=4,ic='bic')['bic_min_order']
    #HQIC
    HQIC = sm.tsa.arma_order_select_ic(timeseries,max_ar=6,\
                 max_ma=4,ic='hqic')['hqic_min_order']
    print('the AIC is{},\nthe BIC is{}\n the HQIC is{}'.format(AIC,BIC,HQIC))
 
def heatmap_AIC(timeseries):
    #设置遍历循环的初始条件，以热力图的形式展示，原理同AIC，BIC，HQIC定阶
    p_min = 0
    q_min = 0
    p_max = 5
    q_max = 5
    d_min = 0
    d_max = 5
    # 创建Dataframe,以BIC准则
    results_aic = pd.DataFrame(index=['AR{}'.format(i) \
                               for i in range(p_min,p_max+1)],\
            columns=['MA{}'.format(i) for i in range(q_min,q_max+1)])
    # itertools.product 返回p,q中的元素的笛卡尔积的元组
    for p,d,q in itertools.product(range(p_min,p_max+1),\
                                   range(d_min,d_max+1),range(q_min,q_max+1)):
        if p==0 and q==0:
            results_aic.loc['AR{}'.format(p), 'MA{}'.format(q)] = np.nan
            continue
        try:
            model = sm.tsa.ARIMA(timeseries, order=(p, d, q))
            results = model.fit()
            #返回不同pq下的model的BIC值
            results_aic.loc['AR{}'.format(p), 'MA{}'.format(q)] = results.aic
        except:
            continue
    results_aic = results_aic[results_aic.columns].astype(float)
    #print(results_bic)
    
    fig, ax = plt.subplots(figsize=(10, 8))
    ax = sns.heatmap(results_aic,
                 #mask=results_aic.isnull(),
                 ax=ax,
                 annot=True, #将数字显示在热力图上
                 fmt='.2f',
                 )
    ax.set_title('AIC')
    plt.show() 

直接利用步骤3的一阶差分来进行定阶，结果如图所示：

                                          

                                              

上面分别是ACF和PACF的图，至于如何定阶不详细叙述了。一般是通过截尾和拖尾来确定阶数。

热力图定阶结果如下所示：

                                   

黑色的位置最好，可以看出p，q取（1,1）（3,1）（1,4）都可以。一般情况下是越小越好。

热力图实现过程参考了了一篇文章，但是博客链接丢失，如果有侵权，请告知，会删除相关部分。

步骤六：模型构建

def ARMA_model(train_data,order): # 训练数据,测试数据,定阶
    arma_model = ARMA(train_data,order) #ARMA模型
    arma = arma_model.fit()#激活模型
    #print(result.summary()) #给出一份模型报告
    
    ############ in-sample ############ 样本内预测 
    in_sample_pred = arma.predict()
 
    ############ out-sample ########## 样本外预测
    #### 样本外预测需要从train_data 样本内的某一个时间节点开始
    #### 利用start和end 控制样外预测 起止时间
    out_sample_pred = arma.predict(start=len(train_data)-2,end = len(train_data)+30, \
                              dynamic=True) 
    #in_sample_pred.plot()
    #train_data.plot()
    return arma,in_sample_pred,out_sample_pred
 
def ARIMA_model(train_data,order):
    arima_model = ARIMA(train_data,order) #ARIMA模型
    arima = arima_model.fit()
    #print(result.summary()) #给出一份模型报告
    
    ########样本内预测#########
    in_sample_pred = arima.predict()
 
 
    ####### 样本外预测##########
    out_sample_pred = arima.predict(start=len(train_data)-2,end = len(train_data)+30, \
                              dynamic=True)
 
    return arima,in_sample_pred,out_sample_pred

预测过程有两种预测方式，一种是样本内的预测（in_sample_pred），一种是样本外的预测(out_sample_pred)。样本内预测就是的是2018-1-1到2018-8-1的。但是要预测的是8-1到9-1的情况，是out-sample预测，一般情况下，out-sample是我们想要的，而不是样本内的预测。

样本外预测是由dynamic参数决定的，特别注意：样本外的预测也要从样本内的某一个时间点开始才能进行预测。因此样本外的预测开始时间要从train_data长度内的某一个时间节点开始。

步骤七：模型评价

主要分为四种方法：（1）QQ图检验残差是否满足正态分布（2）利用D-W检验,检验残差的自相关性（3）计算预测值和真实值的标准差，误差相关等 （4）还原预测序列和测试序列，用图来直观评价模型

def evaluate_model(model,train_data,predict_data):
    
    ###（1）利用QQ图检验残差是否满足正态分布
    resid = model.resid  # 求解模型残差
    plt.figure(figsize=(12,8))
    qqplot(resid,line='q',fit=True)
 
    ###（2）利用D-W检验,检验残差的自相关性
    print('D-W检验值为{}'.format(durbin_watson(resid.values)))
 
    ###（3）利用预测值和真实值的误差检测，这里用的是标准差
    #row_train_data 是从 2018-1-1开始的,经过差分后train_data发生变化
    print('标准差为{}'.format(mean_squared_error(train_data,predict_data,sample_weight=None,\
        multioutput='uniform_average'))) #标准差（均方差）
 
 
def string_toDatetime(string): # 截取时间
    return  datetime.datetime.strptime(string, "%Y-%m-%d %H:%M:%S")
 
#### 绘制图像,查看预测效果
def draw_picture(row_train_data,out_sample_pred,test_data): 
    #print(out_sample_pred)
    # 样本外预测传入 test_data,out_sample_pred
    # 由于预测都是由差分后的平稳序列得出,因此需要对差分后的数据进行还原
    # 还原后绘制同一起点的曲线
 
    #######还原 out_sample_pred #########
    #### out_sample 2018-07-31 
    #### test_data 2018-8-1
    
    ##2018-8-1 00:00 到 2018-9-1 00:00 ###
    #将差分后的序列还原,re_out_sample_pred为还原之后
    re_out_sample_pred = pd.Series(np.array(row_train_data)[-2][0],\
         index=[row_train_data.index[-2]]).append(out_sample_pred[1:]).cumsum()
   
    
    #### 横坐标 
    x = []
    for i in range(32):
        x.append(i+1)
    x = np.array(x)
    
    #### 纵坐标
    y1 = np.array(test_data)
    y2 = np.array(re_out_sample_pred[1:])
    
    #### 画图
    plt.plot(x,y1,color='blue')
    plt.plot(x,y2,color='red')
    plt.show()

（1）qq图如下所示：                   

                                                 

           通过qq图可以看出，残差基本满足了正态分布。

（2）D-W检验结果为：

                                                        

         当D-W检验值接近于2时，不存在自相关性，说明模型较好。

         D-W检验如何数学说明，可以参考下面链接。

         DW值判断准则 - 百度文库

（3）利用标准差来评价模型时，尤其为样本外预测时，注意时间序列的时间对齐。

在利用图来还原预测数据的过程中，主要利用cumsum（）函数，主要作用是累加操作。

  re_out_sample_pred = pd.Series(np.array(row_train_data)[-2][0],\
         index=[row_train_data.index[-2]]).append(out_sample_pred[1:]).cumsum()

   调用以上步骤函数代码如下：

if __name__ == "__main__":
    ## Step 1 and Step 2 都只运行一遍
    genertate_data() # 生成数据
    data_preprocessing() # 1：数据预处理 
    train_data,test_data = Resampling() #:2：数据重采样,返回训练数据和测试数据
    row_train_data = train_data # 保存差分前的序列,为了后面做评估
    ### 差分或者滚动平均,利用步骤4函数确定
    Smooth_data = stationarity(train_data) # 4 差分
 
    for data in zip(Smooth_data,range(4)):# range(4) 用于判断哪种方法 满足平稳性和白噪声  
        if ADF_test(data[0])  and  random_test(data[0]) : # 平稳性和白噪声检测
            train_data = data[0]    # 先用差分，再用平滑,分别对应4个序列
            method = data[1]
            print(method)  #### 如果是差分做的,那么后面ARIMA模型中要使用这个参数
            break 
    ## 三种选择一种即可
    determinate_order_acf(train_data) # ACF定阶
    detetminante_order_AIC(train_data) # BIC 定阶
    heatmap_AIC(train_data) # 热力图 显示
    #### 模型建议和模型评价
    #### order 由差分和定阶给出
    order = (1,1) ## ARMA p,q
    order = (2,1,0) ## ARIMA  p,d,q
    #### 调用模型
    arma,in_sample_pred,out_sample_pred = ARMA_model(train_data,order)
    arima,in_sample_pred,out_sample_pred = ARIMA_model(train_data,order)
    
    #### 模型评价（样本内外均可,此处只用于样本内）
    evaluate_model(arma,train_data,in_sample_pred) # 样本内预测
    evaluate_model(arima,train_data,in_sample_pred) # 样本内预测
   
    #### 画图-差异比较（样本外预测）
    draw_picture(row_train_data,out_sample_pred,test_data)

总结

关于ARMA和ARIMA模型，从数据处理到最后建模实现，就完成了。

但是，里面其实有一个很大的问题，就是当数据不是平稳性的数据的时候，用到了差分法进行处理，用到了dropna()这个函数，这个函数的意思是去掉序列中nan（在这个了里面是0）。因此当序列中两列相邻值相等时，就会去掉前面那一列，因此处理后的数据可能不是按照每一天的数据分布的，但是预测出来的是每一天都存在的。

如果不加dropna（）这个函数的话，定阶那些都会报错（错误信息是存在nan值），但是模型不会报错。因此这块是一个存在的问题，还亟待处理，但是我看了很多的文章，对于这里好像没有过多深入的研究。

整篇博客都是代码实现的，具体的数学公式，自行百度吧~~

2022年1月5日更新：

此版本代码基于statsmodels 0.10.1实现，目前最新statsmodels版本，已经不能像下面这样调用ARMA 和ARIMA模型了。

from statsmodels.tsa.arima_model import ARIMA #ARIMA模型
from statsmodels.tsa.arima_model import ARMA #ARMA模型

新版本为：

from statsmodels.tsa.arima.model import ARIMA #ARIMA模型

关于ARMA好像封装到了 其他模块中，如果想要利用ARMA模型去讨论，还需去官网自行查看。

由于文章过长，完整版代码，会同步到下面公众号，支持提问并第一时间回复。欢迎关注。