時間序列簡介

時間序列 是指將同一統(tǒng)計指標(biāo)的數(shù)值按其先后發(fā)生的時間順序排列而成的數(shù)列氯哮。時間序列分析的主要目的是根據(jù)已有的歷史數(shù)據(jù)對未來進(jìn)行預(yù)測者甲。

常用的時間序列模型

常用的時間序列模型有四種：自回歸模型 AR(p)、移動平均模型 MA(q)慎冤、自回歸移動平均模型 ARMA(p,q)疼燥、自回歸差分移動平均模型 ARIMA(p,d,q), 可以說前三種都是 ARIMA(p,d,q)模型的特殊形式。模型的具體方程可以查找相關(guān)的專業(yè)書籍及網(wǎng)上的資料蚁堤。

AIRIMA模型的建立與預(yù)測

ARIMA 模型是在平穩(wěn)的時間序列基礎(chǔ)上建立起來的醉者，因此時間序列的平穩(wěn)性是建模的重要前提。檢驗時間序列模型平穩(wěn)的方法一般采用 ADF 單位根檢驗?zāi)Ｐ腿z驗披诗。當(dāng)然如果時間序列不穩(wěn)定撬即，也可以通過一些操作去使得時間序列穩(wěn)定（比如取對數(shù)，差分）呈队，然后進(jìn)行 ARIMA 模型預(yù)測剥槐，得到穩(wěn)定的時間序列的預(yù)測結(jié)果，然后對預(yù)測結(jié)果進(jìn)行之前使序列穩(wěn)定的操作的逆操作（取指數(shù)宪摧，差分的逆操作）粒竖，就可以得到原始數(shù)據(jù)的預(yù)測結(jié)果颅崩。

ARIMA 模型實踐

模型具體的理論知識就不再做過多說明了，來個實際的例子吧蕊苗。

ARIMA 模型對湖北省 GDP 的實證分析及預(yù)測

這里的例子是采用了一篇論文的數(shù)據(jù)沿后，【ARIMA模型在湖北省GDP預(yù)測中的應(yīng)用】，可以去中國知網(wǎng)搜索篇名進(jìn)行下載朽砰。

數(shù)據(jù)的平穩(wěn)性處理及檢驗

這里我們用 Python 對數(shù)據(jù)進(jìn)行分析處理建模尖滚。

畫出原始數(shù)據(jù)的時間路徑圖

#-*- coding:utf-8 -*-

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.graphics.tsaplots import acf,pacf,plot_acf,plot_pacf
from statsmodels.tsa.arima_model import ARMA

time_series = pd.Series([151.0, 188.46, 199.38, 219.75, 241.55, 262.58, 328.22, 396.26, 442.04, 517.77, 626.52, 717.08, 824.38, 913.38, 1088.39, 1325.83, 1700.92, 2109.38, 2499.77, 2856.47, 3114.02, 3229.29, 3545.39, 3880.53, 4212.82, 4757.45, 5633.24, 6590.19, 7617.47, 9333.4, 11328.92, 12961.1, 15967.61])
time_series.index = pd.Index(sm.tsa.datetools.dates_from_range('1978','2010'))
time_series.plot(figsize=(12,8))
plt.show()

原始數(shù)據(jù).png

由上圖我們可以看出，這個時間序列是呈指數(shù)形式的瞧柔，波動性比較大漆弄，不是穩(wěn)定的時間序列，一般對于這種指數(shù)形式的數(shù)據(jù)非剃，可以對其取對數(shù)置逻，將其轉(zhuǎn)化為線性趨勢。

time_series = np.log(time_series)
time_series.plot(figsize=(8,6))
plt.show()

對原始數(shù)據(jù)取對數(shù).png

由上圖可以看出备绽，去了對數(shù)之后的時間路徑圖明顯具有線性趨勢券坞，為了確定其穩(wěn)定性，對取對數(shù)后的數(shù)據(jù)進(jìn)行 adf 檢驗

t=sm.tsa.stattools.adfuller(time_series, )
output=pd.DataFrame(index=['Test Statistic Value', "p-value", "Lags Used", "Number of Observations Used","Critical Value(1%)","Critical Value(5%)","Critical Value(10%)"],columns=['value'])
output['value']['Test Statistic Value'] = t[0]
output['value']['p-value'] = t[1]
output['value']['Lags Used'] = t[2]
output['value']['Number of Observations Used'] = t[3]
output['value']['Critical Value(1%)'] = t[4]['1%']
output['value']['Critical Value(5%)'] = t[4]['5%']
output['value']['Critical Value(10%)'] = t[4]['10%']
print(output)

檢驗結(jié)果如下

由上表可知肺素，t 統(tǒng)計量要大于任何置信度的臨界值恨锚，因此認(rèn)為該序列是非平穩(wěn)的，所以再對序列進(jìn)行差分處理倍靡，發(fā)現(xiàn)差分之后的序列基本達(dá)到穩(wěn)定猴伶，如下圖所示，并且通過了 ADF 檢驗塌西，檢驗結(jié)果見下表他挎。

time_series = time_series.diff(1)
time_series = time_series.dropna(how=any)
time_series.plot(figsize=(8,6))
plt.show()
t=sm.tsa.stattools.adfuller(time_series)
output=pd.DataFrame(index=['Test Statistic Value', "p-value", "Lags Used", "Number of Observations Used","Critical Value(1%)","Critical Value(5%)","Critical Value(10%)"],columns=['value'])
output['value']['Test Statistic Value'] = t[0]
output['value']['p-value'] = t[1]
output['value']['Lags Used'] = t[2]
output['value']['Number of Observations Used'] = t[3]
output['value']['Critical Value(1%)'] = t[4]['1%']
output['value']['Critical Value(5%)'] = t[4]['5%']
output['value']['Critical Value(10%)'] = t[4]['10%']
print(output)

差分.png

確定自相關(guān)系數(shù)和平均移動系數(shù)（p,q）

根據(jù)時間序列的識別規(guī)則，采用 ACF 圖捡需、PAC 圖办桨，AIC 準(zhǔn)則（赤道信息量準(zhǔn)則）和 BIC 準(zhǔn)則（貝葉斯準(zhǔn)則）相結(jié)合的方式來確定 ARMA 模型的階數(shù), 應(yīng)當(dāng)選取 AIC 和 BIC 值達(dá)到最小的那一組為理想階數(shù)。

plot_acf(time_series)
plot_pacf(time_series)
plt.show()

r,rac,Q = sm.tsa.acf(time_series, qstat=True)
prac = pacf(time_series,method='ywmle')
table_data = np.c_[range(1,len(r)), r[1:],rac,prac[1:len(rac)+1],Q]
table = pd.DataFrame(table_data, columns=['lag', "AC","Q", "PAC", "Prob(>Q)"])

print(table）

偏自相關(guān)圖.png

自相關(guān)圖.png

ARIMA.png

根據(jù)上面的幾個圖站辉，我們可以先取 p=1, q=2呢撞。進(jìn)行模型估計,結(jié)果見下圖。

p,d,q = (1,1,2)
arma_mod = ARMA(time_series,(p,d,q)).fit(disp=-1,method='mle')
summary = (arma_mod.summary2(alpha=.05, float_format="%.8f"))
print(summary)

模型結(jié)果.png

這里的 p和q 參數(shù)可以調(diào)整饰剥，然后找出最佳的（AIC最小殊霞，BIC最小）汰蓉，經(jīng)過比較绷蹲， p=0,q=1 為理想階數(shù)。
這里有一個自動取 p和q 的函數(shù)顾孽，如果要自動定階的話瘸右，可以采用

(p, q) =(sm.tsa.arma_order_select_ic(dta,max_ar=3,max_ma=3,ic='aic')['aic_min_order'])
#這里需要設(shè)定自動取階的 p和q 的最大值娇跟，即函數(shù)里面的max_ar,和max_ma。ic 參數(shù)表示選用的選取標(biāo)準(zhǔn)太颤，這里設(shè)置的為aic,當(dāng)然也可以用bic苞俘。然后函數(shù)會算出每個 p和q 組合(這里是(0,0)~(3,3)的AIC的值，取其中最小的,這里的結(jié)果是(p=0,q=1)龄章。

殘差和白噪聲檢驗

個人感覺這個就是對模型 ARIMA(0,1,1) 的殘差序列 arma_mod.resid 進(jìn)行 ADF 檢驗吃谣。

arma_mod = ARMA(time_series,(0,1,1)).fit(disp=-1,method='mle')
resid = arma_mod.resid
t=sm.tsa.stattools.adfuller(resid)
output=pd.DataFrame(index=['Test Statistic Value', "p-value", "Lags Used", "Number of Observations Used","Critical Value(1%)","Critical Value(5%)","Critical Value(10%)"],columns=['value'])
output['value']['Test Statistic Value'] = t[0]
output['value']['p-value'] = t[1]
output['value']['Lags Used'] = t[2]
output['value']['Number of Observations Used'] = t[3]
output['value']['Critical Value(1%)'] = t[4]['1%']
output['value']['Critical Value(5%)'] = t[4]['5%']
output['value']['Critical Value(10%)'] = t[4]['10%']
print(output)
#結(jié)果如下
Test Statistic Value           -3.114
p-value                      0.025534
Lags Used                           1
Number of Observations Used        30
Critical Value(1%)           -3.66992
Critical Value(5%)           -2.96407
Critical Value(10%)          -2.62117

當(dāng)然這里也可以畫出 acf 圖和 pacf 圖。

模型預(yù)測

arma_model = sm.tsa.ARMA(time_series,(0,1)).fit(disp=-1,maxiter=100)

predict_data = arma_model.predict(start=str(1979), end=str(2010+3), dynamic = False)

預(yù)測結(jié)果.png

預(yù)測結(jié)果還原

對預(yù)測出來的數(shù)據(jù)做裙，進(jìn)行逆差分操作（由原始數(shù)據(jù)取對數(shù)后的數(shù)據(jù)加上預(yù)測出來的數(shù)據(jù)）岗憋，然后再取指數(shù)即可還原。

預(yù)測結(jié)果還原.png

上圖最后3個為預(yù)測值锚贱，然后查詢2011年到2013年湖北GDP的實際值仔戈，可以進(jìn)行對照

原始數(shù)據(jù)與預(yù)測結(jié)果.png

小結(jié)

從預(yù)測對結(jié)果看，2011年到2013年的預(yù)測結(jié)果和實際的差別不大拧廊。這個模型在短期預(yù)測結(jié)果比較好监徘。模型處理主要還是應(yīng)用了Python 第三方庫 statsmodels 中的模型算法，其中還有很多細(xì)節(jié)吧碾，可以查閱相關(guān)文檔凰盔，這里只是簡單的應(yīng)用了一下，由于代碼都是一小段一小段寫的倦春，很亂户敬，只提供了一些片段供參考。

參考資料

python時間序列分析
時間序列實戰(zhàn)（一）
用ARIMA模型做需求預(yù)測
python 時間序列分析之ARIMA
ARIMA模型文檔