介紹
時間序列（簡稱TS）被認為是分析領域比較少人知道的技能。（我也是幾天前才知道它）。但是你一定知道最近的小型編程馬拉松就是基於時間序列發展起來的，我參加了這項活動去學習了解決時間序列問題的基本步驟，在這兒我要分享給大家。這絕對能幫助你在編程馬拉松中獲得一個合適的模型。
文章之前，我極力推薦大家閱讀《基於R語言的時間序列建模完整教程》A Complete Tutorial on Time Series Modeling in R，它就像這篇文章的前篇。它關注基本概念和基於R語言，我將重點使用這些概念來解決Python編程里面端到端的問題。R語言存在許多關於時間序列的資源，但是很少關於Python的，所以本文將使用Python。
36大數據專稿，原文作者：AARSHAY JAIN 本文由36大數據翻譯，任何不標明譯者和出處以及本文鏈接http://www.36dsj.com/archives/43811的均為侵權。
我們的過程包括下面幾步：
1、時間序列有什么特別之處？
2、在Pandas上傳和加載時間序列（pandas 是基於 Numpy 構建的含有更高級數據結構和工具的數據分析包，類似於 Numpy 的核心是 ndarray，pandas 也是圍繞着 Series 和 DataFrame 兩個核心數據結構展開的 。）
3、如何檢驗時間序列的穩定性？
4、如何令時間序列穩定？
5、時間序列預測。
1、時間序列有什么特別之處？
顧名思義，時間序列是時間間隔不變的情況下收集的時間點集合。這些集合被分析用來了解長期發展趨勢，為了預測未來或者表現分析的其他形式。但是是什么令時間序列與常見的回歸問題的不同？
有兩個原因：
1、時間序列是跟時間有關的。所以基於線性回歸模型的假設：觀察結果是獨立的在這種情況下是不成立的。
2、隨着上升或者下降的趨勢，更多的時間序列出現季節性趨勢的形式，如：特定時間框架的具體變化。即：如果你看到羊毛夾克的銷售上升，你就一定會在冬季做更多銷售。
因為時間序列的固有特性，有各種不同的步驟可以對它進行分析。下文將詳細分析。通過在Python上傳時間序列對象開始。我們將使用飛機乘客數據集。
請記住本文的目的是希望使你熟悉關於時間序列的不同使用方法。本文的例子只是用來方便解釋時間序列對象，我重點關注題目的廣泛性，不會做非常精確的預測。
2、在pandas上傳和加載時間序列
Pandas有專門處理時間序列對象的庫，特別是可以存儲時間信息和允許人們執行快速合作的datatime64(ns)類。從激發所需的庫開始。
import pandas as pd import numpy as np import matplotlib.pylab as plt %matplotlib inline from matplotlib.pylab import rcParams rcParams['figure.figsize'] = 15, 6
現在，我們可以上傳數據集和查看一些最初的行以及列的數據類型。
data = pd.read_csv('AirPassengers.csv')
print data.head()
print '\n Data Types:'
print data.dtypes
數據包含了指定的月份和該月的游客數量。但是時間序列對象的讀取和數據類型的“對象”和“整數類型”的讀取是不一樣的。為了將讀取的數據作為時間序列，我們必須通過特殊的參數讀取csv指令。
dateparse = lambda dates: pd.datetime.strptime(dates, '%Y-%m')
data = pd.read_csv('AirPassengers.csv', parse_dates='Month', index_col='Month',date_parser=dateparse)
print data.head()
我們逐個解釋這些參數：
1、parse_dates:這是指定含有時間數據信息的列。正如上面所說的，列的名稱為“月份”。
2、index_col:使用pandas 的時間序列數據背后的關鍵思想是：目錄成為描述時間數據信息的變量。所以該參數告訴pandas使用“月份”的列作為索引。
3、date_parser：指定將輸入的字符串轉換為可變的時間數據。Pandas默認的數據讀取格式是‘YYYY-MM-DD HH:MM:SS’。如需要讀取的數據沒有默認的格式，就要人工定義。這和dataparse的功能部分相似，這里的定義可以為這一目的服務。
現在我們看到數據有作為索引的時間對象和作為列的乘客（#Passengers）。我們可以通過以下指令再次檢查索引的數據類型。
data.index
注意： dtype=’datetime[ns]’ 確認它是一個時間數據對象。個人而言，我會將列轉換為序列對象，這樣當我每次使用時間序列的時候，就不需要每次都要提及列名稱。當然，這因人而異，如果能令你更好工作，可以使用它作為數據框架。
ts = data[‘#Passengers’] ts.head(10)
在進一步深入前，我會探討一些關於時間序列數據的索引技術。先在序列對象選擇一個特殊值。可以通過以下兩種方式實現：
#1. Specific the index as a string constant: ts['1949-01-01'] #2. Import the datetime library and use 'datetime' function: from datetime import datetime ts[datetime(1949,1,1)]
兩種方法都會返回值“112”，這可以通過先前的結果確認。希望我們可以獲得1949年5月（包括1949年5月）之前的數據。這又可以用以下兩種方法實現：
#1. Specify the entire range: ts['1949-01-01':'1949-05-01'] #2. Use ':' if one of the indices is at ends: ts[:'1949-05-01']
兩種方法都會輸出以下結果：
我們要注意兩點：
- 跟數值索引不一樣，結束索引在這兒是被包含的。比如，如果令一個列表成為索引作為一個數組[:5]，它將在索引返回值-[0,1,2,3,4].在這里索引‘1949-05-01’是包含在結果輸出里面的。
- 目錄必須為了工作區間而進行分類。如果你隨意打亂這些索引，將不能工作。
考慮到可以使用另外一個例子：你需要1949年所有的值。可以這樣做：
ts['1949']
可見，月份部分已經省略。如果你要獲得某月所有的日期，日期部分也可以省略。
現在，讓我們開始分析時間序列。
3、如何檢驗時間序列的穩定性？
如果一個時間序列的統計特征如平均數，方差隨着時間保持不變，我們就可以認為它是穩定的。為什么時間序列的穩定性這么重要？大部分時間序列模型是在假設它是穩定的前提下建立的。直觀地說，我們可以這樣認為，如果一個時間序列隨着時間產生特定的行為，就有很高的可能性認為它在未來的行為是一樣的。同時，根據穩定序列得出的理論是更加成熟的， 也是更容易實現與非穩定序列的比較。
穩定性的確定標准是非常嚴格的。但是，如果時間序列隨着時間產生恆定的統計特征，根據實際目的我們可以假設該序列是穩定的。如下：
- 恆定的平均數
- 恆定的方差
- 不隨時間變化的自協方差
接下來，是關於測試穩定性的方法。首先是簡化坐標和數據，進行可視分析。數據可以通過以下指令定位。
ts['1949']
非常清晰的看到，隨着季節性的變動，飛機乘客的數量總體上是在不斷增長的。但是，不是經常都可以獲得這樣清晰的視覺推論（下文給出相應例子）。所以，更正式的講，我們可以通過下面的方法測試穩定性。
1、繪制滾動統計：我們可以繪制移動平均數和移動方差，觀察它是否隨着時間變化。隨着移動平均數和方差的變化，我認為在任何“t”瞬間，我們都可以獲得去年的移動平均數和方差。如：上一個12個月份。但是，這更多的是一種視覺技術。
2、DF檢驗：這是一種檢查數據穩定性的統計測試。無效假設：時間序列是不穩定的。測試結果由測試統計量和一些置信區間的臨界值組成。如果“測試統計量”少於“臨界值”，我們可以拒絕無效假設，並認為序列是穩定的。
這些概念不是憑直覺得出來的。我推薦大家瀏覽前篇文章。如果你對一些理論數據感興趣，你可以參考Brockwell 和 Davis關於時間序列和預測的介紹的書。這本書的數據很多，但是如果你能讀懂言外之意，你會明白這些概念和正面接觸這些數據。
回到檢查穩定性這件事上，我們將使用滾動數據坐標連同許多DF測試結果，我已經定義了一個需要時間序列作為輸入的函數，為我們生成結果。請注意,我已經繪制標准差來代替方差，為了保持單元和平均數相似。
from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries):
#Determing rolling statistics
rolmean = pd.rolling_mean(timeseries, window=12)
rolstd = pd.rolling_std(timeseries, window=12)
#Plot rolling statistics:
orig = plt.plot(timeseries, color='blue',label='Original')
mean = plt.plot(rolmean, color='red', label='Rolling Mean')
std = plt.plot(rolstd, color='black', label = 'Rolling Std')
plt.legend(loc='best')
plt.title('Rolling Mean & Standard Deviation')
plt.show(block=False)
#Perform Dickey-Fuller test:
print 'Results of Dickey-Fuller Test:'
dftest = adfuller(timeseries, autolag='AIC')
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
print dfoutput
輸入序列開始運行：
基本上是不可能使序列完全穩定，我們只能努力讓它盡可能的穩定。
先讓我們弄明白是什么導致時間序列不穩定。這兒有兩個主要原因。
- 趨勢-隨着時間產生不同的平均值。舉例：在飛機乘客這個案例中，我們看到總體上，飛機乘客的數量是在不斷增長的。
- 季節性-特定時間框架內的變化。舉例：在特定的月份購買汽車的人數會有增加的趨勢，因為車價上漲或者節假日到來。
模型的根本原理或者預測序列的趨勢和季節性，從序列中刪除這些因素，將得到一個穩定的序列。然后統計預測技術可以在這個序列上完成。最后一步是通過運用趨勢和季節性限制倒回到將預測值轉換成原來的區間。
注意：我將探討一系列方法。可能有些對文中情況有用，有些不能。但是我的目的是得到一系列可用方法，而不是僅僅關注目前的問題。
讓我們通過分析趨勢的一部分開始工作吧。
預測&消除趨勢
消除趨勢的第一個方法是轉換。例如,在本例中,我們可以清楚地看到,有一個顯著的趨勢。所以我們可以通過變換，懲罰較高值而不是較小值。這可以采用日志,平方根,立方跟等等。讓我們簡單在這兒轉換一個對數。
ts_log = np.log(ts) plt.plot(ts_log)
在這個簡單的例子中,很容易看到一個向前的數據趨勢。但是它表現的不是很直觀。所以我們可以使用一些技術來估計或對這個趨勢建模,然后將它從序列中刪除。這里有很多方法,最常用的有:
- 聚合-取一段時間的平均值（月/周平均值）
- 平滑-取滾動平均數
- 多項式回歸分析-適合的回歸模型
我在這兒討論將平滑,你也應該嘗試其他可以解決的問題的技術。平滑是指采取滾動估計,即考慮過去的幾個實例。有各種方法可以解決這些問題，但我將主要討論以下兩個。
移動平均數
在這個方法中，根據時間序列的頻率采用“K”連續值的平均數。我們可以采用過去一年的平均數，即過去12個月的平均數。關於確定滾動數據，pandas有特定的功能定義。
moving_avg = pd.rolling_mean(ts_log,12) plt.plot(ts_log) plt.plot(moving_avg, color='red')
紅色線表示了滾動平均數。讓我們從原始序列中減去這個平均數。注意，從我們采用過去12個月的值開始，滾動平均法還沒有對前11個月的值定義。我們可以看到：
ts_log_moving_avg_diff = ts_log - moving_avg ts_log_moving_avg_diff.head(12)
注意前11個月是非數字的，現在讓我們對這11個月降值和檢查這些模塊去測試穩定性。
ts_log_moving_avg_diff.dropna(inplace=True) test_stationarity(ts_log_moving_avg_diff)
這看起來像個更好的序列。滾動平均值出現輕微的變化，但是沒有明顯的趨勢。同時，檢驗統計量比5%的臨界值小，所以我們在95%的置信區間認為它是穩定序列。
但是，這個方法有一個缺陷：要嚴格定義時段。在這種情況下，我們可以采用年平均數，但是對於復雜的情況的像預測股票價格，是很難得到一個數字的。所以，我們采取“加權移動平均法”可以對最近的值賦予更高的權重。關於指定加重這兒有很多技巧。
指數加權移動平均法是很受歡迎的方法，所有的權重被指定給先前的值連同衰減系數。這可以通過pandas實現：
expwighted_avg = pd.ewma(ts_log, halflife=12) plt.plot(ts_log) plt.plot(expwighted_avg, color='red')
注意,這里使用了參數“半衰期”來定義指數衰減量。這只是一個假設,將很大程度上取決於業務領域。其他參數,如跨度和質心也可以用來定義衰減，正如上面鏈接分享的探討。現在讓我們從這個序列轉移，繼續檢查穩定性。
ts_log_ewma_diff = ts_log - expwighted_avg test_stationarity(ts_log_ewma_diff)
這個時間序列有更少的平均值變化和標准差大小變化。同時，檢驗統計量小於1%的臨界值,這比以前的情況好。請注意在這種情況下就不會有遺漏值因為所有的值在一開始就被賦予了權重。所以在運行的時候，它沒有先前的值參與。
消除趨勢和季節性
之前討論來了簡單的趨勢減少技術不能在所有情況下使用，特別是在高季節性情況下。讓我們談論一下兩種消除趨勢和季節性的方法。
- 差分–采用一個特定時間差的差值
- 分解——建立有關趨勢和季節性的模型和從模型中刪除它們。
差分
處理趨勢和季節性的最常見的方法之一就是差分法。在這種方法中,我們采用特定瞬間和它前一個瞬間的不同的觀察結果。這主要是在提高平穩性。pandas可以實現一階差分：
ts_log_diff = ts_log - ts_log.shift() plt.plot(ts_log_diff)
表中可以看出很大程度上減少了趨勢。讓我們通過模塊驗證一下：
ts_log_diff.dropna(inplace=True) test_stationarity(ts_log_diff)
我們可以看到平均數和標准差隨着時間有小的變化。同時，DF檢驗統計量小於10% 的臨界值，因此該時間序列在90%的置信區間上是穩定的。我們同樣可以采取二階或三階差分在具體應用中獲得更好的結果。這些方法你可以自己嘗試。
分解
在這種方法中,趨勢和季節性是分別建模的並倒回到序列的保留部分。我將跳過統計數據，直接給出結果:
from statsmodels.tsa.seasonal import seasonal_decompose decomposition = seasonal_decompose(ts_log) trend = decomposition.trend seasonal = decomposition.seasonal residual = decomposition.resid plt.subplot(411) plt.plot(ts_log, label='Original') plt.legend(loc='best') plt.subplot(412) plt.plot(trend, label='Trend') plt.legend(loc='best') plt.subplot(413) plt.plot(seasonal,label='Seasonality') plt.legend(loc='best') plt.subplot(414) plt.plot(residual, label='Residuals') plt.legend(loc='best') plt.tight_layout()
在這里我們可以看到趨勢,季節性從數據分離，我們可以建立殘差的模型，讓我們檢查殘差的穩定性：
ts_log_decompose = residual ts_log_decompose.dropna(inplace=True) test_stationarity(ts_log_decompose)
DF測試統計量明顯低於1%的臨界值，這樣時間序列是非常接近穩定。你也可以嘗試高級的分解技術產生更好的結果。同時,你應該注意到, 在這種情況下將殘差轉換為原始值對未來數據不是很直觀。
預測時間序列
我們看到不同的技術和它們有效的工作使得時間序列得以穩定。讓我們建立差分后的時間序列模型，因為它是很受歡迎的技術，也相對更容易添加噪音和季節性倒回到預測殘差。在執行趨勢和季節性評估技術上,有兩種情況:
- 不含依賴值的嚴格穩定系列。簡單的情況下,我們可以建立殘差模型作為白噪音（指功率譜密度在整個頻域內均勻分布的噪聲）。但這是非常罕見的。
- 序列含有明顯的依賴值。在這種情況下,我們需要使用一些統計模型像ARIMA（差分自回歸移動平均模型）來預測數據。
讓我給你簡要介紹一下ARIMA，我不會介紹技術細節,但如果你希望更有效地應用它們，你應該理解這些概念的細節。ARIMA代表自回歸整合移動平均數。平穩時間序列的ARIMA預測的只不過是一個線性方程(如線性回歸)。預測依賴於ARIMA模型參數(p d q)。
- 自回歸函數（AR）的條件（p）：AR條件僅僅是因變量的滯后。如：如果P等於5，那么預測x（t）將是x（t-1）。。。（t-5）。
- 移動平均數(MA)的條件(q):MA條件是預測方程的滯后預測錯誤。如：如果q等於5，預測x（t）將是e（t-1）。。。e(t-5),e(i)是移動平均叔在第ith個瞬間和實際值的差值。
- 差分（d）：有非季節性的差值，即這種情況下我們采用一階差分。所以傳遞變量，令d=0或者傳遞原始變量，令d=1。兩種方法得到的結果一樣。
在這里一個重要的問題是如何確定“p”和“q”的值。我們使用兩個坐標來確定這些數字。我們來討論它們。
- 、自相關函數（ACF）：這是時間序列和它自身滯后版本之間的相關性的測試。比如在自相關函數可以比較時間的瞬間‘t1’…’t2’以及序列的瞬間‘t1-5’…’t2-5’ (t1-5和t2 是結束點)。
- 部分自相關函數(PACF):這是時間序列和它自身滯后版本之間的相關性測試，但是是在預測（已經通過比較干預得到解釋）的變量后。如：滯后值為5，它將檢查相關性，但是會刪除從滯后值1到4得到的結果。
時間序列的自回歸函數和部分自回歸函數可以在差分后繪制為:
#ACF and PACF plots: from statsmodels.tsa.stattools import acf, pacf
lag_acf = acf(ts_log_diff, nlags=20) lag_pacf = pacf(ts_log_diff, nlags=20, method='ols')
#Plot ACF:
plt.subplot(121)
plt.plot(lag_acf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.title('Autocorrelation Function')
#Plot PACF:
plt.subplot(122)
plt.plot(lag_pacf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.title('Partial Autocorrelation Function')
plt.tight_layout()
在這個點上,0的每一條邊上的兩條虛線之間是置信區間。這些可以用來確定“p”和“q”的值:
1、p-部分自相關函數表第一次截斷的上層置信區間是滯后值。如果你仔細看，該值是p=2。
2、q- 自相關函數表第一次截斷的上層置信區間是滯后值。如果你仔細看，該值是q=2。
現在，考慮個體以及組合效應建立3個不同的ARIMA模型。我也會發布各自的RSS（是一種描述和同步網站內容的格式，是使用最廣泛的XML應用）。請注意,這里的RSS是指殘差值,而不是實際序列。
首先，我們需要上傳ARIMA模型。
from statsmodels.tsa.arima_model import ARIMA
p,d,q值可以指定使用ARIMA的命令參數即采用一個元(p,d,q)。建立三種情況下的模型：
自回歸（AR）模型：
model = ARIMA(ts_log, order=(2, 1, 0))
results_AR = model.fit(disp=-1)
plt.plot(ts_log_diff)
plt.plot(results_AR.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_AR.fittedvalues-ts_log_diff)**2))
組合模型
model = ARIMA(ts_log, order=(0, 1, 2))
results_MA = model.fit(disp=-1)
plt.plot(ts_log_diff)
plt.plot(results_MA.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_MA.fittedvalues-ts_log_diff)**2))
移動平均數（MA ）模型
model = ARIMA(ts_log, order=(2, 1, 2))
results_ARIMA = model.fit(disp=-1)
plt.plot(ts_log_diff)
plt.plot(results_ARIMA.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_ARIMA.fittedvalues-ts_log_diff)**2))
在這里我們可以看到,自回歸函數模型和移動平均數模型幾乎有相同的RSS，但相結合效果顯著更好。現在,我們只剩下最后一步,即把這些值倒回到原始區間。
倒回到原始區間
既然組合模型獲得更好的結果，讓我們將它倒回原始值，看看它如何執行。第一步是作為一個獨立的序列，存儲預測結果,觀察它。
predictions_ARIMA_diff = pd.Series(results_ARIMA.fittedvalues, copy=True) print predictions_ARIMA_diff.head()
注意,這些是從‘1949-02-01’開始,而不是第一個月。為什么?這是因為我們將第一個月份取為滯后值，一月前面沒有可以減去的元素。將差分轉換為對數尺度的方法是這些差值連續地添加到基本值。一個簡單的方法就是首先確定索引的累計總和,然后將其添加到基本值。累計總和可以在下面找到:
predictions_ARIMA_diff_cumsum = predictions_ARIMA_diff.cumsum() print predictions_ARIMA_diff_cumsum.head()
你可以在頭腦使用之前的輸出結果進行回算，檢查這些是否正確的。接下來我們將它們添加到基本值。為此我們將使用所有的值創建一個序列作為基本值，並添加差值。我們這樣做：
predictions_ARIMA_log = pd.Series(ts_log.ix[0], index=ts_log.index) predictions_ARIMA_log = predictions_ARIMA_log.add(predictions_ARIMA_diff_cumsum,fill_value=0) predictions_ARIMA_log.head()
第一個元素是基本值本身，從基本值開始值累計添加。最后一步是將指數與原序列比較。
predictions_ARIMA = np.exp(predictions_ARIMA_log)
plt.plot(ts)
plt.plot(predictions_ARIMA)
plt.title('RMSE: %.4f'% np.sqrt(sum((predictions_ARIMA-ts)**2)/len(ts)))
最后我們獲得一個原始區間的預測結果。雖然不是一個很好的預測。但是你獲得了思路對嗎?現在,我把它留個你去進一步改進，做一個更好的方案。
最后注意
在本文中,我試圖提供你們一個標准方法去解決時間序列問題。這個不可能達到一個更好的時間,因為今天是我們的小型編程馬拉松，挑戰你們是否可以解決類似的問題。我們廣泛的討論了穩定性的概念和最終的預測殘差。這是一個漫長的過程,我跳過了一些統計細節,我鼓勵大家使用這些作為參考材料。如果你不想復制粘貼,你可以從我的GitHub庫下載iPython筆記本的代碼。