一、什么是時間序列
時間序列簡單的說就是各時間點上形成的數值序列，時間序列分析就是通過觀察歷史數據預測未來的值。
在這里需要強調一點的是，時間序列分析並不是關於時間的回歸，它主要是研究自身的變化規律的（這里不考慮含外生變量的時間序列）。
環境配置
python作為科學計算的利器，當然也有相關分析的包:statsmodels中tsa模塊，當然這個包和SAS、R是比不了，但是python有另一個神器：pandas！pandas在時間序列上的應用，能簡化我們很多的工作。這兩個包pip就能安裝。
數據准備
許多時間序列分析一樣，本文同樣使用航空乘客數據（AirPassengers.csv）作為樣例。下載鏈接。
用pandas操作時間序列
# -*- 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='Month') df.index = pd.to_datetime(df.index) # 將字符串索引轉換成時間索引 ts = df['Passengers'] # 生成pd.Series對象 # 查看數據格式 print(ts.head()) print(ts.head().index)
不知道為啥，時間隔1000，不過沒什么影響。
查看某日的值既可以使用字符串作為索引，又可以直接使用時間對象作為索引
ts['2049-01-01'] ts[datetime(2049,1,1)]
兩者的返回值都是第一個序列值：112
如果要查看某一年的數據，pandas也能非常方便的實現
ts['2049']
切片操作：
ts['2049-1' : '2049-6']
注意時間索引的切片操作起點和尾部都是包含的，這點與數值索引有所不同
pandas還有很多方便的時間序列函數，在后面的實際應用中在進行說明。
二、時間序列分析
1. 基本模型
自回歸移動平均模型(ARMA(p，q))是時間序列中最為重要的模型之一，它主要由兩部分組成： AR代表p階自回歸過程，MA代表q階移動平均過程，其公式如下：
依據模型的形式、特性及自相關和偏自相關函數的特征，總結如下：
在時間序列中，ARIMA模型是在ARMA模型的基礎上多了差分的操作。
2. 平穩性檢驗
我們知道序列平穩性是進行時間序列分析的前提條件，很多人都會有疑問，為什么要滿足平穩性的要求呢？在大數定理和中心定理中要求樣本同分布（這里同分布等價於時間序列中的平穩性），而我們的建模過程中有很多都是建立在大數定理和中心極限定理的前提條件下的，如果它不滿足，得到的許多結論都是不可靠的。以虛假回歸為例，當響應變量和輸入變量都平穩時，我們用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.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() 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，說明並不能拒絕原假設。
3. 平穩性處理
由前面的分析可知，該序列是不平穩的，然而平穩性是時間序列分析的前提條件，故我們需要對不平穩的序列進行處理將其轉換成平穩的序列。
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中，對差分的支持不是很好，它不支持高階和多階差分。
不過沒關系，我們有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 trend.plot() #seasonal.plot() #residual.plot() plt.show()
得到不同的分解成分后，就可以使用時間序列模型對各個成分進行擬合，當然也可以選擇其他預測方法。我曾經用過小波對時序數據進行過分解，然后分別采用時間序列擬合，效果還不錯。由於我對小波的理解不是很好，只能簡單的調用接口，如果有誰對小波、傅里葉、卡爾曼理解得比較透，可以將時序數據進行更加准確的分解，由於分解后的時序數據避免了他們在建模時的交叉影響，所以我相信它將有助於預測准確性的提高。
4. 模型識別
在前面的分析可知，該序列具有明顯的年周期與長期成分。對於年周期成分我們使用窗口為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')
5. 樣本擬合
模型擬合完后，我們就可以對其進行預測了。由於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，感覺還過得去。
6. 完善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]) # 預測 from statsmodels.tsa.arima_model import ARMA 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) # 繪制結果 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()
是不是發現這里的預測結果和上一篇的使用12階移動平均的預測結果一模一樣。這是因為12階移動平均加上一階差分與直接12階差分是等價的關系，后者是前者數值的12倍，這個應該不難推導。
7. BIC准則
8. 滾動預測
9. 模型序列化
三、總結
與其它統計語言（SAS和R）相比，python在統計分析這塊還顯得不那么“專業”。隨着numpy、pandas、scipy、sklearn、gensim、statsmodels等包的推動，我相信也祝願python在數據分析這塊越來越好。