python實踐——時間序列分析建模理論及代碼實現

python進階教程

機器學習

深度學習關注

進入正文

聲明：本文所講的時間序列分析並不是指pandas的時間序列處理方法，pandas時間序列處理更缺確切地說時間序列的可視化、窗口移動操作等一些操作的統稱。本文所講的時間序列分析指的是一種算法，一種通過序列本身所潛在的規律去預測未來某個時刻可能發生的狀況。

python時間序列分析

目錄

一 什麼是時間序列

二 為什麼要用python

三 環境配置

四 時間序列分析

4.1 基本模型

4.2 pandas時間序列操作

4.3 平穩性檢驗

4.4 平穩性處理

4.5 模型識別

4.6 樣本擬合

4.7 完善ARIMA模型

4.8 滾動預測

4.9 模型序列化

五 總結

01

什麼是時間序列

時間序列簡單的說就是各時間點上形成的數值序列，時間序列分析就是通過觀察歷史數據預測未來的值。在這裡需要強調一點的是，時間序列分析並不是關於時間的迴歸，它主要是研究自身的變化規律的（這裡不考慮含外生變量的時間序列）。

02

為什麼用python

現在各類編程語言都有專門的時間序列分析開源框架或者是專門處理的相關程序庫，處理時間序列最多的是使用MATLAB，因為它現成的工具函數，只需要使用幾句話就可以實現一個時間序列建模模型。而我之所以選擇python來實現，用兩個字總結“情懷”，愛屋及烏，個人比較喜歡python，就用python擼了。能做時間序列的軟件很多，SAS、R、SPSS、Eviews甚至matlab等等，實際工作中應用得比較多的應該還是SAS和R，前者推薦王燕寫的《應用時間序列分析》，後者推薦“基於R語言的時間序列建模完整教程”這篇博文（翻譯版）。python作為科學計算的利器，當然也有相關分析的包:statsmodels中tsa模塊，當然這個包和SAS、R是比不了，但是python有另一個神器：pandas！pandas在時間序列上的應用，能簡化我們很多的工作。

03

環境配置

python推薦直接裝Anaconda，它集成了許多科學計算包，有一些包自己手動去裝還是挺費勁的。statsmodels需要自己去安裝，這裡我推薦使用0.6的穩定版，0.7及其以上的版本能在github上找到，該版本在安裝時會用C編譯好，所以修改底層的一些代碼將不會起作用。

04

時間序列分析

第四章

4.1 基本模型

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

依據模型的形式、特性及自相關和偏自相關函數的特徵，總結如下：

在時間序列中，ARIMA模型是在ARMA模型的基礎上多了差分的操作。

第四章

4.2 pandas時間序列操作

大熊貓真的很可愛，這裡簡單介紹一下它在時間序列上的可愛之處。和許多時間序列分析一樣，本文同樣使用航空乘客數據（AirPassengers.csv）作為樣例。

數據讀取：

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

import numpy as np

import pandas as pd

from datetime import datetime

import matplotlib.pylab as plt

# 讀取數據，pd.read_csv默認生成DataFrame對象，需將其轉換成Series對象

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

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

ts = df['x'] # 生成pd.Series對象

# 查看數據格式

ts.head()

ts.head().index

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

ts['1949-01-01']

ts[datetime(1949,1,1)]

兩者的返回值都是第一個序列值：112

如果要查看某一年的數據，pandas也能非常方便的實現

ts['1949']

切片操作：

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

注意時間索引的切片操作起點和尾部都是包含的，這點與數值索引有所不同

pandas還有很多方便的時間序列函數，在後面的實際應用中在進行說明。

第四章

我們知道序列平穩性是進行時間序列分析的前提條件，很多人都會有疑問，為什麼要滿足平穩性的要求呢？在大數定理和中心定理中要求樣本同分布（這裡同分布等價於時間序列中的平穩性），而我們的建模過程中有很多都是建立在大數定理和中心極限定理的前提條件下的，如果它不滿足，得到的許多結論都是不可靠的。以虛假迴歸為例，當響應變量和輸入變量都平穩時，我們用t統計量檢驗標準化係數的顯著性。而當響應變量和輸入變量不平穩時，其標準化係數不在滿足t分佈，這時再用t檢驗來進行顯著性分析，導致拒絕原假設的概率增加，即容易犯第一類錯誤，從而得出錯誤的結論。

平穩時間序列有兩種定義：嚴平穩和寬平穩

嚴平穩顧名思義，是一種條件非常苛刻的平穩性，它要求序列隨著時間的推移，其統計性質保持不變。對於任意的τ，其聯合概率密度函數滿足：

嚴平穩的條件只是理論上的存在，現實中用得比較多的是寬平穩的條件。

寬平穩也叫弱平穩或者二階平穩（均值和方差平穩），它應滿足：

• 常數均值
• 常數方差
• 常數自協方差

平穩性檢驗：觀察法和單位根檢驗法

基於此，我寫了一個名為test_stationarity的統計性檢驗模塊，以便將某些統計檢驗結果更加直觀的展現出來。

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

from statsmodels.tsa.stattools import adfuller

import pandas as pd

import matplotlib.pyplot as plt

import numpy as np

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# 移動平均圖

def draw_trend(timeSeries, size):

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

# 對size個數據進行移動平均

rol_mean = timeSeries.rolling(window=size).mean()

# 對size個數據進行加權移動平均

rol_weighted_mean = pd.ewma(timeSeries, span=size)

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()

def draw_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

'''

def testStationarity(ts):

dftest = adfuller(ts)

# 對上述函數求得的值進行語義描述

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] = value

return dfoutput

# 自相關和偏相關圖，默認階數為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()

觀察法，通俗的說就是通過觀察序列的趨勢圖與相關圖是否隨著時間的變化呈現出某種規律。所謂的規律就是時間序列經常提到的週期性因素，現實中遇到得比較多的是線性週期成分，這類週期成分可以採用差分或者移動平均來解決，而對於非線性週期成分的處理相對比較複雜，需要採用某些分解的方法。下圖為航空數據的線性圖，可以明顯的看出它具有年週期成分和長期趨勢成分。平穩序列的自相關係數會快速衰減，下面的自相關圖並不能體現出該特徵，所以我們有理由相信該序列是不平穩的。

單位根檢驗：ADF是一種常用的單位根檢驗方法，他的原假設為序列具有單位根，即非平穩，對於一個平穩的時序數據，就需要在給定的置信水平上顯著，拒絕原假設。ADF只是單位根檢驗的方法之一，如果想採用其他檢驗方法，可以安裝第三方包arch，裡面提供了更加全面的單位根檢驗方法，個人還是比較鍾情ADF檢驗。以下為檢驗結果，其p值大於0.99，說明並不能拒絕原假設。

第四章

4.4 平穩性處理

由前面的分析可知，該序列是不平穩的，然而平穩性是時間序列分析的前提條件，故我們需要對不平穩的序列進行處理將其轉換成平穩的序列。

a. 對數變換

對數變換主要是為了減小數據的振動幅度，使其線性規律更加明顯（我是這麼理解的時間序列模型大部分都是線性的，為了儘量降低非線性的因素，需要對其進行預處理，也許我理解的不對）。對數變換相當於增加了一個懲罰機制，數據越大其懲罰越大，數據越小懲罰越小。這裡強調一下，變換的序列需要滿足大於0，小於0的數據不存在對數變換。

ts_log = np.log(ts)

test_stationarity.draw_ts(ts_log)

b. 平滑法

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

移動平均即利用一定時間間隔內的平均值作為某一期的估計值，而指數平均則是用變權的方法來計算均值。

test_stationarity.draw_trend(ts_log, 12)

從上圖可以發現窗口為12的移動平均能較好的剔除年週期性因素，而指數平均法是對週期內的數據進行了加權，能在一定程度上減小年週期因素，但並不能完全剔除，如要完全剔除可以進一步進行差分操作。

c. 差分

時間序列最常用來剔除週期性因素的方法當屬差分了，它主要是對等週期間隔的數據進行線性求減。前面我們說過，ARIMA模型相對ARMA模型，僅多了差分操作，ARIMA模型幾乎是所有時間序列軟件都支持的，差分的實現與還原都非常方便。而statsmodel中，對差分的支持不是很好，它不支持高階和多階差分，為什麼不支持，這裡引用作者的說法：

作者大概的意思是說預測方法中並沒有解決高於2階的差分，有沒有感覺很牽強，不過沒關係，我們有pandas。我們可以先用pandas將序列差分好，然後在對差分好的序列進行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)

從上面的統計檢驗結果可以看出，經過12階差分和1階差分後，該序列滿足平穩性的要求了。

d. 分解

所謂分解就是將時序數據分離成不同的成分。statsmodels使用的X-11分解過程，它主要將時序數據分離成長期趨勢、季節趨勢和隨機成分。與其它統計軟件一樣，statsmodels也支持兩類分解模型，加法模型和乘法模型，這裡我只實現加法，乘法只需將model的參數設置為"multiplicative"即可。

from statsmodels.tsa.seasonal import seasonal_decompose

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

trend = decomposition.trend

seasonal = decomposition.seasonal

residual = decomposition.resid

得到不同的分解成分後，就可以使用時間序列模型對各個成分進行擬合，當然也可以選擇其他預測方法。我曾經用過小波對時序數據進行過分解，然後分別採用時間序列擬合，效果還不錯。由於我對小波的理解不是很好，只能簡單的調用接口，如果有誰對小波、傅里葉、卡爾曼理解得比較透，可以將時序數據進行更加準確的分解，由於分解後的時序數據避免了他們在建模時的交叉影響，所以我相信它將有助於預測準確性的提高。

第四章

4.5 模型識別

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

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)

觀察其統計量發現該序列在置信水平為95%的區間下並不顯著，我們對其進行再次一階差分。再次差分後的序列其自相關具有快速衰減的特點，t統計量在99%的置信水平下是顯著的，這裡我不再做詳細說明。

ts_diff_2 = ts_diff_1.diff(1)

ts_diff_2.dropna(inplace=True)

數據平穩後，需要對模型定階，即確定p、q的階數。觀察上圖，發現自相關和偏相係數都存在拖尾的特點，並且他們都具有明顯的一階相關性，所以我們設定p=1, q=1。下面就可以使用ARMA模型進行數據擬合了。這裡我不使用ARIMA(ts_diff_1, order=(1, 1, 1))進行擬合，是因為含有差分操作時，預測結果還原老出問題，至今還沒弄明白。

from statsmodels.tsa.arima_model import ARMA

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

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

第四章

4.6 樣本擬合

模型擬合完後，我們就可以對其進行預測了。由於ARMA擬合的是經過相關預處理後的數據，故其預測值需要通過相關逆變換進行還原。

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)

# 對數還原

log_recover = np.exp(rol_recover)

log_recover.dropna(inplace=True)

我們使用均方根誤差（RMSE）來評估模型樣本內擬合的好壞。利用該準則進行判別時，需要剔除“非預測”數據的影響。

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

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，感覺還過得去。

第四章

4.7 完善ARIMA模型

前面提到statsmodels裡面的ARIMA模塊不支持高階差分，我們的做法是將差分分離出來，但是這樣會多了一步人工還原的操作。基於上述問題，我將差分過程進行了封裝，使序列能按照指定的差分列表依次進行差分，並相應的構造了一個還原的方法，實現差分序列的自動還原。

# 差分操作

def diff_ts(ts, d):

global shift_ts_list

# 動態預測第二日的值時所需要的差分序列

global last_data_shift_list

shift_ts_list = []

last_data_shift_list = []

tmp_ts = ts

for i in d:

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)

return tmp_ts

# 還原操作

def predict_diff_recover(predict_value, d):

if isinstance(predict_value, float):

tmp_data = predict_value

for i in range(len(d)):

tmp_data = tmp_data + last_data_shift_list[-i-1]

elif isinstance(predict_value, np.ndarray):

tmp_data = predict_value[0]

for i in range(len(d)):

tmp_data = tmp_data + last_data_shift_list[-i-1]

else:

tmp_data = predict_value

for i in range(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

現在我們直接使用差分的方法進行數據處理，並以同樣的過程進行數據預測與還原。

diffed_ts = diff_ts(ts_log, d=[12, 1])

model = arima_model(diffed_ts)

model.certain_model(1, 1)

predict_ts = model.properModel.predict()

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

log_recover = np.exp(diff_recover_ts)

是不是發現這裡的預測結果和上一篇的使用12階移動平均的預測結果一模一樣。這是因為12階移動平均加上一階差分與直接12階差分是等價的關係，後者是前者數值的12倍，這個應該不難推導。

對於個數不多的時序數據，我們可以通過觀察自相關圖和偏相關圖來進行模型識別，倘若我們要分析的時序數據量較多，例如要預測每隻股票的走勢，我們就不可能逐個去調參了。這時我們可以依據BIC準則識別模型的p, q值，通常認為BIC值越小的模型相對更優。這裡我簡單介紹一下BIC準則，它綜合考慮了殘差大小和自變量的個數，殘差越小BIC值越小，自變量個數越多BIC值越大。個人覺得BIC準則就是對模型過擬合設定了一個標準（過擬合這東西應該以辯證的眼光看待）。

def proper_model(data_ts, maxLag):

init_bic = sys.maxint

init_p = 0

init_q = 0

init_properModel = None

for p in np.arange(maxLag):

for q in np.arange(maxLag):

model = ARMA(data_ts, order=(p, q))

try:

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

except:

continue

bic = results_ARMA.bic

if bic < init_bic:

init_p = p

init_q = q

init_properModel = results_ARMA

init_bic = bic

return init_bic, init_p, init_q, init_properModel

相對最優參數識別結果：BIC: -1090.44209358 p: 0 q: 1 ， RMSE:11.8817198331。我們發現模型自動識別的參數要比我手動選取的參數更優。

第四章

4.8 滾動預測

所謂滾動預測是指通過添加最新的數據預測第二天的值。對於一個穩定的預測模型，不需要每天都去擬合，我們可以給他設定一個閥值，例如每週擬合一次，該期間只需通過添加最新的數據實現滾動預測即可。基於此我編寫了一個名為arima_model的類，主要包含模型自動識別方法，滾動預測的功能，詳細代碼可以查看附錄。數據的動態添加：

from dateutil.relativedelta import relativedelta

def _add_new_data(ts, dat, type='day'):

if type == 'day':

new_index = ts.index[-1] + relativedelta(days=1)

elif type == 'month':

new_index = ts.index[-1] + relativedelta(months=1)

ts[new_index] = dat

def add_today_data(model, ts, data, d, type='day'):

_add_new_data(ts, data, type) # 為原始序列添加數據

# 為滯後序列添加新值

d_ts = diff_ts(ts, d)

model.add_today_data(d_ts[-1], type)

def forecast_next_day_data(model, type='day'):

if model == None:

raise ValueError('No model fit before')

fc = model.forecast_next_day_value(type)

return predict_diff_recover(fc, [12, 1])

現在我們就可以使用滾動預測的方法向外預測了，取1957年之前的數據作為訓練數據，其後的數據作為測試，並設定模型每第七天就會重新擬合一次。這裡的diffed_ts對象會隨著add_today_data方法自動添加數據，這是由於它與add_today_data方法中的d_ts指向的同一對象，該對象會動態的添加數據。

ts_train = ts_log[:'1956-12']

ts_test = ts_log['1957-1':]

diffed_ts = diff_ts(ts_train, [12, 1])

forecast_list = []

for i, dta in enumerate(ts_test):

if i%7 == 0:

model = arima_model(diffed_ts)

model.certain_model(1, 1)

forecast_data = forecast_next_day_data(model, type='month')

forecast_list.append(forecast_data)

add_today_data(model, ts_train, dta, [12, 1], type='month')

predict_ts = pd.Series(data=forecast_list, index=ts['1957-1':].index)

log_recover = np.exp(predict_ts)

original_ts = ts['1957-1':]

動態預測的均方根誤差為：14.6479，與前面樣本內擬合的均方根誤差相差不大，說明模型並沒有過擬合，並且整體預測效果都較好。

第四章

4.9 模型化預測

在進行動態預測時，我們不希望將整個模型一直在內存中運行，而是希望有新的數據到來時才啟動該模型。這時我們就應該把整個模型從內存導出到硬盤中，而序列化正好能滿足該要求。序列化最常用的就是使用json模塊了，但是它是時間對象支持得不是很好，有人對json模塊進行了拓展以使得支持時間對象，這裡我們不採用該方法，我們使用pickle模塊，它和json的接口基本相同，有興趣的可以去查看一下。我在實際應用中是採用的面向對象的編程，它的序列化主要是將類的屬性序列化即可，而在面向過程的編程中，模型序列化需要將需要序列化的對象公有化，這樣會使得對前面函數的參數改動較大，我不在詳細闡述該過程。

04

總結

與其它統計語言相比，python在統計分析這塊還顯得不那麼“專業”。隨著numpy、pandas、scipy、sklearn、gensim、statsmodels等包的推動，我相信也祝願python在數據分析這塊越來越好。與SAS和R相比，python的時間序列模塊還不是很成熟，我這裡僅起到拋磚引玉的作用，希望各位能人志士能貢獻自己的力量，使其更加完善。實際應用中我全是面向過程來編寫的，為了闡述方便，我用面向過程重新羅列了一遍，實在感覺很不方便。原本打算分三篇來寫的，還有一部分實際應用的部分，不打算再寫了，還請大家原諒。實際應用主要是具體問題具體分析，這當中第一步就是要查詢問題，這步花的時間往往會比較多，然後再是解決問題。以我前面項目遇到的問題為例，當時遇到了以下幾個典型的問題：1.週期長度不恆定的週期成分，例如每月的1號具有周期性，但每月1號與1號之間的時間間隔是不相等的；2.含有缺失值以及含有記錄為0的情況無法進行對數變換；3.節假日的影響等等。！

看了這篇文章是不是對時間序列分析有著跟深入的理解和認知呢？其實時間序列分析是一個非常龐大的課題，有很多專業的書籍專門用一本書來講解時間序列分析理論，本文只是一個簡單的實現，如果有興趣，小夥伴可以好好學一學時間序列分析哦。後面還會更新一些系列文章，包括python的設計模式、深度神經網絡的原理講解、數據結構的Python實現，想要每天不斷學習的小夥伴都可以關注一下哦，認真看完本系列文章，我相信一定會大有收穫的。歡迎小夥伴們繼續關注和支持。如果你有需要，就可以私信（

感謝大家持續關注！並且可以轉發文章讓大家受益！

獲取方式：私信老師(學習)就可以拿到了！

閱讀更多 馬士兵尚學堂 的文章

關鍵字: 歷史 R語言 機器學習