一、什么是時(shí)間序列

時(shí)間序列簡單的說就是各時(shí)間點(diǎn)上形成的數(shù)值序列，時(shí)間序列分析就是通過觀察歷史數(shù)據(jù)預(yù)測未來的值。

在這里需要強(qiáng)調(diào)一點(diǎn)的是，時(shí)間序列分析并不是關(guān)于時(shí)間的回歸，它主要是研究自身的變化規(guī)律的(這里不考慮含外生變量的時(shí)間序列)。

環(huán)境配置

python作為科學(xué)計(jì)算的利器，當(dāng)然也有相關(guān)分析的包:statsmodels中tsa模塊，當(dāng)然這個(gè)包和SAS、R是比不了，但是python有另一個(gè)神器：pandas！pandas在時(shí)間序列上的應(yīng)用，能簡化我們很多的工作。這兩個(gè)包pip就能安裝。

數(shù)據(jù)準(zhǔn)備

許多時(shí)間序列分析一樣，本文同樣使用航空乘客數(shù)據(jù)(AirPassengers.csv)作為樣例。下載鏈接。

用pandas操作時(shí)間序列

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

importnumpy as npimportpandas as pdfrom datetime importdatetimeimportmatplotlib.pylab as plt#讀取數(shù)據(jù)，pd.read_csv默認(rèn)生成DataFrame對象，需將其轉(zhuǎn)換成Series對象

df = pd.read_csv('AirPassengers.csv', encoding='utf-8', index_col='Month')

df.index= pd.to_datetime(df.index) #將字符串索引轉(zhuǎn)換成時(shí)間索引

ts = df['Passengers'] #生成pd.Series對象#查看數(shù)據(jù)格式

print(ts.head())print(ts.head().index)

不知道為啥，時(shí)間隔1000，不過沒什么影響。

查看某日的值既可以使用字符串作為索引，又可以直接使用時(shí)間對象作為索引

ts['2049-01-01']

ts[datetime(2049,1,1)]

兩者的返回值都是第一個(gè)序列值：112

如果要查看某一年的數(shù)據(jù)，pandas也能非常方便的實(shí)現(xiàn)

ts['2049']

切片操作：

ts['2049-1' : '2049-6']

注意時(shí)間索引的切片操作起點(diǎn)和尾部都是包含的，這點(diǎn)與數(shù)值索引有所不同

pandas還有很多方便的時(shí)間序列函數(shù)，在后面的實(shí)際應(yīng)用中在進(jìn)行說明。

二、時(shí)間序列分析

1. 基本模型

自回歸移動平均模型(ARMA(p，q))是時(shí)間序列中最為重要的模型之一，它主要由兩部分組成： AR代表p階自回歸過程，MA代表q階移動平均過程，其公式如下：

依據(jù)模型的形式、特性及自相關(guān)和偏自相關(guān)函數(shù)的特征，總結(jié)如下：

在時(shí)間序列中，ARIMA模型是在ARMA模型的基礎(chǔ)上多了差分的操作。

2. 平穩(wěn)性檢驗(yàn)

我們知道序列平穩(wěn)性是進(jìn)行時(shí)間序列分析的前提條件，很多人都會有疑問，為什么要滿足平穩(wěn)性的要求呢？在大數(shù)定理和中心定理中要求樣本同分布(這里同分布等價(jià)于時(shí)間序列中的平穩(wěn)性)，而我們的建模過程中有很多都是建立在大數(shù)定理和中心極限定理的前提條件下的，如果它不滿足，得到的許多結(jié)論都是不可靠的。以虛假回歸為例，當(dāng)響應(yīng)變量和輸入變量都平穩(wěn)時(shí)，我們用t統(tǒng)計(jì)量檢驗(yàn)標(biāo)準(zhǔn)化系數(shù)的顯著性。而當(dāng)響應(yīng)變量和輸入變量不平穩(wěn)時(shí)，其標(biāo)準(zhǔn)化系數(shù)不在滿足t分布，這時(shí)再用t檢驗(yàn)來進(jìn)行顯著性分析，導(dǎo)致拒絕原假設(shè)的概率增加，即容易犯第一類錯(cuò)誤，從而得出錯(cuò)誤的結(jié)論。

平穩(wěn)時(shí)間序列有兩種定義：嚴(yán)平穩(wěn)和寬平穩(wěn)

嚴(yán)平穩(wěn)顧名思義，是一種條件非常苛刻的平穩(wěn)性，它要求序列隨著時(shí)間的推移，其統(tǒng)計(jì)性質(zhì)保持不變。對于任意的τ，其聯(lián)合概率密度函數(shù)滿足：

嚴(yán)平穩(wěn)的條件只是理論上的存在，現(xiàn)實(shí)中用得比較多的是寬平穩(wěn)的條件。

寬平穩(wěn)也叫弱平穩(wěn)或者二階平穩(wěn)(均值和方差平穩(wěn))，它應(yīng)滿足：

常數(shù)均值

常數(shù)方差

常數(shù)自協(xié)方差

平穩(wěn)性檢驗(yàn)：觀察法和單位根檢驗(yàn)法

基于此，我寫了一個(gè)名為test_stationarity的統(tǒng)計(jì)性檢驗(yàn)?zāi)K，以便將某些統(tǒng)計(jì)檢驗(yàn)結(jié)果更加直觀的展現(xiàn)出來。

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

from statsmodels.tsa.stattools importadfullerimportpandas as pdimportmatplotlib.pyplot as pltimportnumpy as npfrom statsmodels.graphics.tsaplots importplot_acf, plot_pacf#移動平均圖

defdraw_trend(timeSeries, size):

f= plt.figure(facecolor='white')#對size個(gè)數(shù)據(jù)進(jìn)行移動平均

rol_mean = timeSeries.rolling(window=size).mean()#對size個(gè)數(shù)據(jù)進(jìn)行加權(quán)移動平均

rol_weighted_mean = pd.DataFrame.ewm(timeSeries, span=size).mean()

timeSeries.plot(color='blue', label='Original')

rolmean.plot(color='red', label='Rolling Mean')

rol_weighted_mean.plot(color='black', label='Weighted Rolling Mean')

plt.legend(loc='best')

plt.title('Rolling Mean')

plt.show()defdraw_ts(timeSeries):

f= plt.figure(facecolor='white')

timeSeries.plot(color='blue')

plt.show()'''Unit Root Test

The null hypothesis of the Augmented Dickey-Fuller is that there is a unit

root, with the alternative that there is no unit root. That is to say the

bigger the p-value the more reason we assert that there is a unit root'''

deftestStationarity(ts):

dftest=adfuller(ts)#對上述函數(shù)求得的值進(jìn)行語義描述

dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])for key,value in dftest[4].items():

dfoutput['Critical Value (%s)'%key] =valuereturndfoutput#自相關(guān)和偏相關(guān)圖，默認(rèn)階數(shù)為31階

def draw_acf_pacf(ts, lags=31):

f= plt.figure(facecolor='white')

ax1= f.add_subplot(211)

plot_acf(ts, lags=31, ax=ax1)

ax2= f.add_subplot(212)

plot_pacf(ts, lags=31, ax=ax2)

plt.show()

觀察法，通俗的說就是通過觀察序列的趨勢圖與相關(guān)圖是否隨著時(shí)間的變化呈現(xiàn)出某種規(guī)律。所謂的規(guī)律就是時(shí)間序列經(jīng)常提到的周期性因素，現(xiàn)實(shí)中遇到得比較多的是線性周期成分，這類周期成分可以采用差分或者移動平均來解決，而對于非線性周期成分的處理相對比較復(fù)雜，需要采用某些分解的方法。下圖為航空數(shù)據(jù)的線性圖，可以明顯的看出它具有年周期成分和長期趨勢成分。平穩(wěn)序列的自相關(guān)系數(shù)會快速衰減，下面的自相關(guān)圖并不能體現(xiàn)出該特征，所以我們有理由相信該序列是不平穩(wěn)的。

單位根檢驗(yàn)：ADF是一種常用的單位根檢驗(yàn)方法，他的原假設(shè)為序列具有單位根，即非平穩(wěn)，對于一個(gè)平穩(wěn)的時(shí)序數(shù)據(jù)，就需要在給定的置信水平上顯著，拒絕原假設(shè)。ADF只是單位根檢驗(yàn)的方法之一，如果想采用其他檢驗(yàn)方法，可以安裝第三方包arch，里面提供了更加全面的單位根檢驗(yàn)方法，個(gè)人還是比較鐘情ADF檢驗(yàn)。以下為檢驗(yàn)結(jié)果，其p值大于0.99，說明并不能拒絕原假設(shè)。

3. 平穩(wěn)性處理

由前面的分析可知，該序列是不平穩(wěn)的，然而平穩(wěn)性是時(shí)間序列分析的前提條件，故我們需要對不平穩(wěn)的序列進(jìn)行處理將其轉(zhuǎn)換成平穩(wěn)的序列。

a. 對數(shù)變換

對數(shù)變換主要是為了減小數(shù)據(jù)的振動幅度，使其線性規(guī)律更加明顯(我是這么理解的時(shí)間序列模型大部分都是線性的，為了盡量降低非線性的因素，需要對其進(jìn)行預(yù)處理，也許我理解的不對)。對數(shù)變換相當(dāng)于增加了一個(gè)懲罰機(jī)制，數(shù)據(jù)越大其懲罰越大，數(shù)據(jù)越小懲罰越小。這里強(qiáng)調(diào)一下，變換的序列需要滿足大于0，小于0的數(shù)據(jù)不存在對數(shù)變換。

ts_log =np.log(ts)

test_stationarity.draw_ts(ts_log)

b. 平滑法

根據(jù)平滑技術(shù)的不同，平滑法具體分為移動平均法和指數(shù)平均法。

移動平均即利用一定時(shí)間間隔內(nèi)的平均值作為某一期的估計(jì)值，而指數(shù)平均則是用變權(quán)的方法來計(jì)算均值

test_stationarity.draw_trend(ts_log, 12)

從上圖可以發(fā)現(xiàn)窗口為12的移動平均能較好的剔除年周期性因素，而指數(shù)平均法是對周期內(nèi)的數(shù)據(jù)進(jìn)行了加權(quán)，能在一定程度上減小年周期因素，但并不能完全剔除，如要完全剔除可以進(jìn)一步進(jìn)行差分操作。

c. 差分

時(shí)間序列最常用來剔除周期性因素的方法當(dāng)屬差分了，它主要是對等周期間隔的數(shù)據(jù)進(jìn)行線性求減。前面我們說過，ARIMA模型相對ARMA模型，僅多了差分操作，ARIMA模型幾乎是所有時(shí)間序列軟件都支持的，差分的實(shí)現(xiàn)與還原都非常方便。而statsmodel中，對差分的支持不是很好，它不支持高階和多階差分。

不過沒關(guān)系，我們有pandas。我們可以先用pandas將序列差分好，然后在對差分好的序列進(jìn)行ARIMA擬合，只不過這樣后面會多了一步人工還原的工作。

diff_12 = ts_log.diff(12)

diff_12.dropna(inplace=True)

diff_12_1= diff_12.diff(1)

diff_12_1.dropna(inplace=True)

test_stationarity.testStationarity(diff_12_1)

從上面的統(tǒng)計(jì)檢驗(yàn)結(jié)果可以看出，經(jīng)過12階差分和1階差分后，該序列滿足平穩(wěn)性的要求了。

d. 分解

所謂分解就是將時(shí)序數(shù)據(jù)分離成不同的成分。statsmodels使用的X-11分解過程，它主要將時(shí)序數(shù)據(jù)分離成長期趨勢、季節(jié)趨勢和隨機(jī)成分。與其它統(tǒng)計(jì)軟件一樣，statsmodels也支持兩類分解模型，加法模型和乘法模型，這里我只實(shí)現(xiàn)加法，乘法只需將model的參數(shù)設(shè)置為”multiplicative”即可。

from statsmodels.tsa.seasonal importseasonal_decompose

decomposition= seasonal_decompose(ts_log, model="additive")

trend=decomposition.trend

seasonal=decomposition.seasonal

residual=decomposition.resid

trend.plot()#seasonal.plot()#residual.plot()

plt.show()

得到不同的分解成分后，就可以使用時(shí)間序列模型對各個(gè)成分進(jìn)行擬合，當(dāng)然也可以選擇其他預(yù)測方法。我曾經(jīng)用過小波對時(shí)序數(shù)據(jù)進(jìn)行過分解，然后分別采用時(shí)間序列擬合，效果還不錯(cuò)。由于我對小波的理解不是很好，只能簡單的調(diào)用接口，如果有誰對小波、傅里葉、卡爾曼理解得比較透，可以將時(shí)序數(shù)據(jù)進(jìn)行更加準(zhǔn)確的分解，由于分解后的時(shí)序數(shù)據(jù)避免了他們在建模時(shí)的交叉影響，所以我相信它將有助于預(yù)測準(zhǔn)確性的提高。

4. 模型識別

在前面的分析可知，該序列具有明顯的年周期與長期成分。對于年周期成分我們使用窗口為12的移動平進(jìn)行處理，對于長期趨勢成分我們采用1階差分來進(jìn)行處理。

rol_mean = ts_log.rolling(window=12).mean()

rol_mean.dropna(inplace=True)

ts_diff_1= rol_mean.diff(1)

ts_diff_1.dropna(inplace=True)

test_stationarity.testStationarity(ts_diff_1)

觀察其統(tǒng)計(jì)量發(fā)現(xiàn)該序列在置信水平為95%的區(qū)間下并不顯著，我們對其進(jìn)行再次一階差分。再次差分后的序列其自相關(guān)具有快速衰減的特點(diǎn)，t統(tǒng)計(jì)量在99%的置信水平下是顯著的，這里我不再做詳細(xì)說明。

ts_diff_2 = ts_diff_1.diff(1)

ts_diff_2.dropna(inplace=True)

數(shù)據(jù)平穩(wěn)后，需要對模型定階，即確定p、q的階數(shù)。觀察上圖，發(fā)現(xiàn)自相關(guān)和偏相系數(shù)都存在拖尾的特點(diǎn)，并且他們都具有明顯的一階相關(guān)性，所以我們設(shè)定p=1, q=1。下面就可以使用ARMA模型進(jìn)行數(shù)據(jù)擬合了。這里我不使用ARIMA(ts_diff_1, order=(1, 1, 1))進(jìn)行擬合，是因?yàn)楹胁罘植僮鲿r(shí)，預(yù)測結(jié)果還原老出問題，至今還沒弄明白。

from statsmodels.tsa.arima_model importARMA

model= ARMA(ts_diff_2, order=(1, 1))

result_arma= model.fit( disp=-1, method='css')

5. 樣本擬合

模型擬合完后，我們就可以對其進(jìn)行預(yù)測了。由于ARMA擬合的是經(jīng)過相關(guān)預(yù)處理后的數(shù)據(jù)，故其預(yù)測值需要通過相關(guān)逆變換進(jìn)行還原。

predict_ts =result_arma.predict()#一階差分還原

diff_shift_ts = ts_diff_1.shift(1)

diff_recover_1=predict_ts.add(diff_shift_ts)#再次一階差分還原

rol_shift_ts = rol_mean.shift(1)

diff_recover=diff_recover_1.add(rol_shift_ts)#移動平均還原

rol_sum = ts_log.rolling(window=11).sum()

rol_recover= diff_recover*12 - rol_sum.shift(1)#對數(shù)還原

log_recover =np.exp(rol_recover)

log_recover.dropna(inplace=True)

我們使用均方根誤差(RMSE)來評估模型樣本內(nèi)擬合的好壞。利用該準(zhǔn)則進(jìn)行判別時(shí)，需要剔除“非預(yù)測”數(shù)據(jù)的影響。

ts = ts[log_recover.index] #過濾沒有預(yù)測的記錄

plt.figure(facecolor='white')

log_recover.plot(color='blue', label='Predict')

ts.plot(color='red', label='Original')

plt.legend(loc='best')

plt.title('RMSE: %.4f'% np.sqrt(sum((log_recover-ts)**2)/ts.size))

plt.show()

觀察上圖的擬合效果，均方根誤差為11.8828，感覺還過得去。

6. 完善ARIMA模型

前面提到statsmodels里面的ARIMA模塊不支持高階差分，我們的做法是將差分分離出來，但是這樣會多了一步人工還原的操作。基于上述問題，我將差分過程進(jìn)行了封裝，使序列能按照指定的差分列表依次進(jìn)行差分，并相應(yīng)的構(gòu)造了一個(gè)還原的方法，實(shí)現(xiàn)差分序列的自動還原。

#差分操作

defdiff_ts(ts, d):globalshift_ts_list#動態(tài)預(yù)測第二日的值時(shí)所需要的差分序列

globallast_data_shift_list

shift_ts_list=[]

last_data_shift_list=[]

tmp_ts=tsfor i ind:

last_data_shift_list.append(tmp_ts[-i])print(last_data_shift_list)

shift_ts=tmp_ts.shift(i)

shift_ts_list.append(shift_ts)

tmp_ts= tmp_ts -shift_ts

tmp_ts.dropna(inplace=True)returntmp_ts#還原操作

defpredict_diff_recover(predict_value, d):ifisinstance(predict_value, float):

tmp_data=predict_valuefor i inrange(len(d)):

tmp_data= tmp_data + last_data_shift_list[-i-1]elifisinstance(predict_value, np.ndarray):

tmp_data=predict_value[0]for i inrange(len(d)):

tmp_data= tmp_data + last_data_shift_list[-i-1]else:

tmp_data=predict_valuefor i inrange(len(d)):try:

tmp_data= tmp_data.add(shift_ts_list[-i-1])except:raise ValueError('What you input is not pd.Series type!')

tmp_data.dropna(inplace=True)return tmp_data

現(xiàn)在我們直接使用差分的方法進(jìn)行數(shù)據(jù)處理，并以同樣的過程進(jìn)行數(shù)據(jù)預(yù)測與還原。

#差分

diffed_ts = diff_ts(ts_log, d=[12, 1])#預(yù)測

from statsmodels.tsa.arima_model importARMA

model= ARMA(diffed_ts, order=(1, 1))

result_arma= model.fit( disp=-1, method='css')

predict_ts=result_arma.predict()#還原

predict_ts =model.properModel.predict()

diff_recover_ts= predict_diff_recover(predict_ts, d=[12, 1])

log_recover=np.exp(diff_recover_ts)#繪制結(jié)果

ts = ts[log_recover.index] #過濾沒有預(yù)測的記錄

plt.figure(facecolor='white')

log_recover.plot(color='blue', label='Predict')

ts.plot(color='red', label='Original')

plt.legend(loc='best')

plt.title('RMSE: %.4f'% np.sqrt(sum((log_recover-ts)**2)/ts.size))

plt.show()

是不是發(fā)現(xiàn)這里的預(yù)測結(jié)果和上一篇的使用12階移動平均的預(yù)測結(jié)果一模一樣。這是因?yàn)?2階移動平均加上一階差分與直接12階差分是等價(jià)的關(guān)系，后者是前者數(shù)值的12倍，這個(gè)應(yīng)該不難推導(dǎo)。

7. BIC準(zhǔn)則

8. 滾動預(yù)測

9. 模型序列化

三、總結(jié)

與其它統(tǒng)計(jì)語言(SAS和R)相比，python在統(tǒng)計(jì)分析這塊還顯得不那么“專業(yè)”。隨著numpy、pandas、scipy、sklearn、gensim、statsmodels等包的推動，我相信也祝愿python在數(shù)據(jù)分析這塊越來越好。

總結(jié)

以上是生活随笔為你收集整理的python时间序列分析航空旅人_python时间序列分析的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。