本文將介紹使用Python來完成時間序列分析ARIMA模型的完整步驟與流程

時間序列分析概念

時間序列分析 是統計學中的一個非常重要的分支，是以概率論與數理統計為基礎、計算機應用為技術支撐，迅速發展起來的一種應用性很強的科學方法，時間序列是變數按時間間隔的順序而下形成的隨機變數序列，大量自然界、社會經濟等領域的統計指標都依年、季、月或日統計其指標值，隨著時間的推移，形成了統計指標的時間序列，例如，股價指數、物價指數、GDP和產品銷售量等等都屬于時間序列，

建立模型基本步驟

ARIMA模型建模實戰

匯入模塊

• 小白可以參考我哦

import sys
import os
import pandas as pd
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
import statsmodels.tsa.api as smt
from statsmodels.tsa.stattools import adfuller 
from statsmodels.stats.diagnostic import acorr_ljungbox 
from statsmodels.graphics.api import qqplot
import matplotlib.pylab as plt
from matplotlib.pylab import style
style.use('ggplot')
from arch.unitroot import ADF
import warnings
warnings.filterwarnings("ignore")
pd.set_option('display.float_format', lambda x: '%.5f' % x) 
np.set_printoptions(precision=5, suppress=True) 
pd.set_option('display.max_columns', 100)
pd.set_option('display.max_rows', 100)

%matplotlib inline
"""中文顯示問題"""
plt.rcParams['font.family'] = ['sans-serif']
plt.rcParams['font.sans-serif'] = ['SimHei']

加載資料

data = pd.read_excel("data.xlsx",index_col="年份",parse_dates=True)
data.head()

平穩性檢驗

時序圖

data["diff1"] = data["xt"].diff(1).dropna()
data["diff2"] = data["diff1"].diff(1).dropna()
"""由于差分后會有缺失值，使用該方法需要填充缺失值"""
#     data1 = data.loc[:,["xt","diff1","diff2"]]
#     data1.plot(subplots=True, figsize=(18, 12),title="差分圖")
fig = plt.figure(figsize=(20,6))
ax1 = fig.add_subplot(131)
ax1.plot(data['xt'],label="原序列")
plt.legend()
ax2 = fig.add_subplot(132)
ax2.plot(data['diff1'].dropna().dropna(),label="一階差分")
plt.legend()
ax3 = fig.add_subplot(133)
ax3.plot(data['diff2'].dropna(),label="二階差分")
plt.legend(loc="upper left")
plt.show()

時序圖檢驗 - - 全靠肉眼的判斷和判斷人的經驗，不同的人看到同樣的圖形，很可能會給出不同的判斷，因此我們需要一個更有說服力、更加客觀的統計方法來幫助我們檢驗時間序列的平穩性，這種方法，就是單位根檢驗，

單位根檢驗

print("單位根檢驗:\n")
print(ADF(data.diff1.dropna()))  

單位根檢驗:

   Augmented Dickey-Fuller Results   
=====================================
Test Statistic                 -3.156
P-value                         0.023
Lags                                0
-------------------------------------

Trend: Constant
Critical Values: -3.63 (1%), -2.95 (5%), -2.61 (10%)
Null Hypothesis: The process contains a unit root.
Alternative Hypothesis: The process is weakly stationary.

單位根檢驗 - -對其一階差分進行單位根檢驗，得到：1%、%5、%10不同程度拒絕原假設的統計值和ADF Test result的比較，本資料中，P-value 為 0.023,接近0，ADF Test result同時小于5%、10%即說明很好地拒絕該假設，本資料中，ADF結果為-3.156，拒絕原假設，即一階差分后資料是平穩的，

白噪聲檢驗

判斷序列是否為非白噪聲序列

from statsmodels.stats.diagnostic import acorr_ljungbox
acorr_ljungbox(data.diff1.dropna(), lags = [i for i in range(1,12)],boxpierce=True)

(array([11.30402, 13.03896, 13.37637, 14.24184, 14.6937 , 15.33042,
        16.36099, 16.76433, 18.15565, 18.16275, 18.21663]),
 array([0.00077, 0.00147, 0.00389, 0.00656, 0.01175, 0.01784, 0.02202,
        0.03266, 0.03341, 0.05228, 0.07669]),
 array([10.4116 , 11.96391, 12.25693, 12.98574, 13.35437, 13.85704,
        14.64353, 14.94072, 15.92929, 15.93415, 15.9696 ]),
 array([0.00125, 0.00252, 0.00655, 0.01135, 0.02027, 0.03127, 0.04085,
        0.06031, 0.06837, 0.10153, 0.14226]))

通過P<α,拒絕原假設，故差分后的序列是平穩的非白噪聲序列，可以進行下一步建模

模型定階

現在我們已經得到一個平穩的時間序列，接來下就是選擇合適的ARIMA模型，即ARIMA模型中合適的p,q，
第一步我們要先檢查平穩時間序列的自相關圖和偏自相關圖，通過sm.graphics.tsa.plot_acf和sm.graphics.tsa.plot_pacf得到圖形

從一階差分序列的自相關圖和偏自相關圖可以發現:

1. 自相關圖拖尾或一階截尾
2. 偏自相關圖一階截尾,

• 所以我們可以建立ARIMA(1,1,0)、ARIMA(1,1,1)、ARIMA(0,1,1)模型，

模型優化

• 其中L是在該模型下的最大似然，n是資料數量，k是模型的變數個數，
  python代碼如下：

arma_mod20 = sm.tsa.ARIMA(data["xt"],(1,1,0)).fit()
arma_mod30 = sm.tsa.ARIMA(data["xt"],(0,1,1)).fit()
arma_mod40 = sm.tsa.ARIMA(data["xt"],(1,1,0)).fit()
values = [[arma_mod20.aic,arma_mod20.bic,arma_mod20.hqic],[arma_mod30.aic,arma_mod30.bic,arma_mod30.hqic],[arma_mod40.aic,arma_mod40.bic,arma_mod40.hqic]]
df = pd.DataFrame(values,index=["AR(1,1,0)","MA(0,1,1)","ARMA(1,1,0)"],columns=["AIC","BIC","hqic"])
df

• 構造這些統計量所遵循的統計思想是一致的，就是在考慮擬合殘差的同時，依自變數個數施加“懲罰”，但要注意的是，這些準則不能說明某一個模型的精確度，也即是說，對于三個模型Ａ，Ｂ，Ｃ，我們能夠判斷出B模型是最好的，但不能保證B模型能夠很好地刻畫資料

引數估計

from statsmodels.tsa.arima_model import ARIMA
model = ARIMA(data["xt"], order=(0,1,1))
result = model.fit()
print(result.summary())

                             ARIMA Model Results                              
==============================================================================
Dep. Variable:                   D.xt   No. Observations:                   36
Model:                 ARIMA(0, 1, 1)   Log Likelihood                -122.987
Method:                       css-mle   S.D. of innovations              7.309
Date:                Tue, 22 Dec 2020   AIC                            251.973
Time:                        09:11:55   BIC                            256.724
Sample:                    01-01-1953   HQIC                           253.631
                         - 01-01-1988                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          4.9956      2.014      2.481      0.013       1.048       8.943
ma.L1.D.xt     0.6710      0.165      4.071      0.000       0.348       0.994
                                    Roots                                    
=============================================================================
                  Real          Imaginary           Modulus         Frequency
-----------------------------------------------------------------------------
MA.1           -1.4902           +0.0000j            1.4902            0.5000
-----------------------------------------------------------------------------

模型檢驗

引數的顯著性檢驗

P<α，拒絕原假設，認為該引數顯著非零MA(2)模型擬合該序列，殘差序列已實作白噪聲

模型的顯著性檢驗

resid = result.resid#殘差
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
fig = qqplot(resid, line='q', ax=ax, fit=True)

qq圖顯示，我們看到紅色的KDE線與N(0,1)平行，這是殘留物正太分布的良好指標，說明殘差序列是白噪聲序列，模型的資訊的提取充分，當讓大家也可以使用前面介紹的檢驗白噪聲的方法LB統計量來檢驗

ARIMA(0,1,1)模型擬合該序列，殘差序列已實作白噪聲，且引數均顯著非零，說明AR(0,11)模型是該序列的有效擬合模型

模型預測

pred = result.predict('1988', '1990',dynamic=True, typ='levels')
print (pred)

1988-01-01   278.35527
1989-01-01   283.35088
1990-01-01   288.34649
Freq: AS-JAN, dtype: float64

plt.figure(figsize=(12, 8))
plt.xticks(rotation=45)
plt.plot(pred)
plt.plot(data.xt)
plt.show()

到這里就結束了，如果對你有幫助你，歡迎點贊關注，你的點贊對我很重要