《數據分析實戰-托馬茲.卓巴斯》讀書筆記第7章-時間序列技術（ARMA模型、ARIMA模型）

 python學習筆記-目錄索引

 
第7章探索瞭如何處理和理解時間序列數據，並創建ARMA模型以及ARIMA模型。注意：我在本章花的時間較長，主要是對dataframe結構不熟。

本章會介紹處理、分析和預測時間序列數據的各類技術。會學習如下技巧：
·在Python中如何處理日期對象
·理解時間序列數據
·平滑並轉換觀測值
·過濾時間序列數據
·移除趨勢和季節性
·使用ARMA和ARIMA模型預測將來

7.1導論

時間序列隨處可見；若是分析股票市場、太陽黑子，或河流，你就是在觀察隨時間延展的現象。數據科學家在職業生涯中處理時間序列數據老是不可避免的。本章中，咱們會遇到對時間序列進行處理、分析和構建模型的多種技巧。
本章中用到的數據集來自網上河流的Web文檔：http://ftp.uni-bayreuth.de/math/statlib/datasets/riverflow。這個文檔本質上就是一個shell腳本，爲本章建立數據集。要從文檔中建立原始文件，你可使用Windows下的Cygwin或者Mac/Linux下的Terminal，執行下述命令（假設你將文檔保存在riverflows.webarchive）：

/*  
sh riverflows.webarchive
*/

邀月建議：安裝cygwin巨麻煩，仍是用安裝好的CentOS虛擬機執行一下。

7.2在Python中如何處理日期對象

時間序列是以某個時間間隔進行採樣獲得的數據，例如，記錄每秒的車速。拿到這樣的數據，咱們能夠輕鬆估算通過的距離（假設觀測值加總併除以3600）或者汽車的加速度（計算兩個觀測值之間的差別）。能夠直接用pandas處理時間序列數據。
準備：需裝好pandas、NumPy和Matplotlib。

步驟：從Web文檔開始，咱們進行清理，並造成兩個數據集：美國河（http://www.theameri-canriver.com）和哥倫比亞河（http://www.ecy.wa.gov/programs/wr/cwp/cwpfactmap.html）。用pandas讀取時間序列數據集很簡單（ts_handlingData.py文件）：

 1 import numpy as np
 2 import pandas as pd
 3 import pandas.tseries.offsets as ofst
 4 import matplotlib
 5 import matplotlib.pyplot as plt
 6 
 7 # change the font size
 8 matplotlib.rc('xtick', labelsize=9)
 9 matplotlib.rc('ytick', labelsize=9)
10 matplotlib.rc('font', size=14)
11 
12 # files we'll be working with
13 files=['american.csv', 'columbia.csv']
14 
15 # folder with data
16 data_folder = '../../Data/Chapter07/'
17 
18 # colors
19 colors = ['#FF6600', '#000000', '#29407C', '#660000']
20 
21 # read the data
22 american = pd.read_csv(data_folder + files[0],
23     index_col=0, parse_dates=[0],
24     header=0, names=['','american_flow'])
25 
26 columbia = pd.read_csv(data_folder + files[1],
27     index_col=0, parse_dates=[0],
28     header=0, names=['','columbia_flow'])
29 
30 # combine the datasets
31 riverFlows = american.combine_first(columbia)
32 
33 # periods aren't equal in the two datasets so find the overlap
34 # find the first month where the flow is missing for american
35 idx_american = riverFlows \
36     .index[riverFlows['american_flow'].apply(np.isnan)].min()
37 
38 # find the last month where the flow is missing for columbia
39 idx_columbia = riverFlows \
40     .index[riverFlows['columbia_flow'].apply(np.isnan)].max()
41 
42 # truncate the time series
43 riverFlows = riverFlows.truncate(
44     before=idx_columbia + ofst.DateOffset(months=1),
45     after=idx_american - ofst.DateOffset(months=1))

Tips:

/*
 Traceback (most recent call last):
  File "D:\Java2018\practicalDataAnalysis\Codes\Chapter07\ts_handlingData.py", line 49, in <module>
    o.write(riverFlows.to_csv(ignore_index=True))
TypeError: to_csv() got an unexpected keyword argument 'ignore_index'

D:\Java2018\practicalDataAnalysis\Codes\Chapter07\ts_handlingData.py:80: FutureWarning: how in .resample() is deprecated
the new syntax is .resample(...).mean()
  year = riverFlows.resample('A', how='mean')
 
 */

解決方案：

/*
# year = riverFlows.resample('A', how='mean')
year = riverFlows.resample('A').mean()

# o.write(riverFlows.to_csv(ignore_index=True))
  o.write(riverFlows.to_csv(index=True)) */

原理：首先，咱們引入全部必需的模塊：pandas和NumPy。咱們將獲得兩個文件：american.csv和columbia.csv。它們都位於data_folder。
咱們使用已經熟悉的pandas的.read_csv（...）方法。先讀入american.csv文件。指定index_col=0，讓方法將第一列做爲索引。要讓pandas將某列當成日期處理，咱們顯式地命令.read_csv（...）方法將列做爲日期解析（parse_dates）。
將兩個文件合併成一個數據集。而後改動列名：咱們告訴方法，這裏沒有頭部，而且將自行提供名字。注意第一列不須要任何名字，它將被轉換成索引。咱們用一樣的方式讀入哥倫比亞河的數據。
讀入兩個文件後，將它們聯繫在一塊兒。pandas的.combine_first（...）方法操做第一個數據集，插入哥倫比亞河數據集的列。
若是沒有改變DataFrame的列名，.combine_first（...）方法將使用被調用DataFrame的數據來填充調用者DataFrame的空隙。
兩個文件的時期不一樣，可是有重疊部分：美國河數據從1906年到1960年，哥倫比亞河數據從1933年到1969年。查看一下重疊的時期；咱們鏈接的數據集只有1933年到1960年的數據。
首先，找到美國河沒有數據的最先日期（american_flow列）。檢查riverFlows的索引，選取american_flow值不是數字的全部日期；使用.apply（...）方法並使用NumPy的.isnan方法檢查DataFrame的元素。作完這個以後，選取序列中的最小日期。
而對於columbia_flow，咱們找的是沒有數據的最晚日期。和處理美國河數據相似，咱們先取出全部數據不是數字的日期，而後選取最大值。
.truncate（...）方法能夠根據DatetimeIndex從DataFrame中移除數據。
DatetimeIndex是數字組成的不可變數組。內部由大整數表示，但看上去是日期-時間對象：既有日期部分也有時間部分的對象。

參考http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DatetimeIndex.html。
咱們向.truncate（...）方法傳遞兩個參數。before參數指定了要捨棄哪一個日期以前的記錄，after參數指定了保留數據的最後一個日期。
idx_...對象保存了至少有一列沒有數據的日期的最小值和最大值。不過，若是將這些日期傳入.truncate（...）方法，咱們也會選出一樣沒有數據的極值點。應對這種狀況，咱們用.DateOffset（...）方法將日期移個位。咱們只挪一個月。
若是想更深刻了解.DateOffset（...）方法，可參考http://pandas.pydata.org/pandas-docs/stable/timeseries.html#dateoffset-objects。
最後將鏈接後的數據集保存到文件。（更多信息參考本書1.2節）

/*
Index of riverFlows
DatetimeIndex(['1933-01-31', '1933-02-28', '1933-03-31', '1933-04-30',
               '1933-05-31', '1933-06-30', '1933-07-31', '1933-08-31',
               '1933-09-30', '1933-10-31',
               ...
               '1960-03-31', '1960-04-30', '1960-05-31', '1960-06-30',
               '1960-07-31', '1960-08-31', '1960-09-30', '1960-10-31',
               '1960-11-30', '1960-12-31'],
              dtype='datetime64[ns]', name='', length=336, freq=None)

csv_read['1933':'1934-06']
            american_flow  columbia_flow
                                        
1933-01-31        10.7887          24.10
1933-02-28        14.6115          20.81
1933-03-31        19.6236          22.96
1933-04-30        21.9739          37.66
1933-05-31        28.0054         118.93
1933-06-30        66.0632         331.31
1933-07-31       113.4373         399.27
1933-08-31       162.0007         250.89
1933-09-30       156.6771         116.10
1933-10-31        17.9246          69.38
1933-11-30         7.0792          52.95
1933-12-31         4.0493          40.21
1934-01-31        11.5816          35.40
1934-02-28        18.5192          28.88
1934-03-31        53.8586          33.41
1934-04-30        75.8608         102.22
1934-05-31        89.3963         259.67
1934-06-30       116.2973         390.77

Shifting one month forward
            american_flow  columbia_flow
                                        
1933-02-28        10.7887          24.10
1933-03-31        14.6115          20.81
1933-04-30        19.6236          22.96
1933-05-31        21.9739          37.66
1933-06-30        28.0054         118.93
1933-07-31        66.0632         331.31

Shifting one year forward
            american_flow  columbia_flow
                                        
1934-01-31        10.7887          24.10
1934-02-28        14.6115          20.81
1934-03-31        19.6236          22.96
1934-04-30        21.9739          37.66
1934-05-31        28.0054         118.93
1934-06-30        66.0632         331.31

Averaging by quarter
            american_flow  columbia_flow
                                        
1933-03-31      15.007933      22.623333
1933-06-30      38.680833     162.633333

Averaging by half a year
            american_flow  columbia_flow
                                        
1933-01-31      10.788700      24.100000
1933-07-31      43.952483     155.156667

Averaging by year
            american_flow  columbia_flow
                                        
1933-12-31      51.852875     123.714167
1934-12-31      44.334742     128.226667 */

View Code

更多:從時間序列DataFrame中選取數據很直接：你仍然可使用已知的DataFrame取子集的方法（好比，使用元素索引的riverFlows[0：10]或者.head（...）方法）。不過，有了DataFrame和DatetimeIndex索引，你也能夠傳日期：

print(riverFlows['1933':'1934-06'])

這行命令會打印出1933年到1934年6月的全部記錄。注意，儘管有每個月的數據，咱們依然能夠傳入一個年份；pandas會選取相應的觀測值。
有時候，你須要時間平移後的觀測值，好比，測量儀的內部時鐘沒調整好，致使觀測值偏了一小時，或者一我的爲錯誤讓觀測值偏了一年。這能夠用.shift（...）方法輕鬆糾正：

1 # shifting the data
2 by_month = riverFlows.shift(1, freq='M')
3 by_year = riverFlows.shift(12, freq='M')

傳給.shift（...）方法的第一個參數指定了要偏移的數字，freq參數指定了頻率（本例中是一個月）。
有時，你想在季末或年底時計算平均值（或總和）。可經過.resample（...）方法輕鬆作到：

1 # averaging by quarter
2 quarter = riverFlows.resample('Q').mean()
3 # averaging by half a year
4 half = riverFlows.resample('6M').mean()
5 # averaging by year
6 year = riverFlows.resample('A').mean()

第一個參數定義了頻率：Q意味着季末（三月、六月、九月和十二月），how指定了怎麼處理觀測值（本例中使用mean，由於對咱們的數據來講平均值比較有意義，但若是要處理銷售額數據，你也許會指定sum）。
A是年底，6M是半年底，即十二月和六月。
另外，pandas可輕鬆繪製時間序列的圖：

 1 import matplotlib
 2 import matplotlib.pyplot as plt
 3 matplotlib.rc('font', size=14)
 4 
 5 # colors
 6 colors = ['#FF6600', '#000000', '#29407C', '#660000']
 7 
 8 # plot time series
 9 # monthly time series
10 riverFlows.plot(title='Monthly river flows', color=colors)
11 plt.savefig(data_folder + '/charts/monthly_riverFlows.png',
12     dpi=300)
13 plt.close() 

先引入Matplotlib模塊，更改x軸和y軸的數字與標題。
Matplotlib的配置項可參考http://matplotlib.org/users/customizing.html。
而後指定顏色；這個常量在全部圖表中都會用到。
.plot（...）方法繪圖。咱們定義標題，容許方法從顏色列表中挑選顏色。最後，將圖保存到文件並關閉，以釋放內存。
月度和季度圖表以下：

查看月度圖表，你也許會以爲兩條河的流動沒什麼變化；而季度圖則展現了不一樣之處：

理解時間序列數據
但願從數據集中看出什麼時，一件基本事實是：若是不理解你要處理的數據，你不可能構建一個成功的統計模型。
準備：需裝好pandas、Statsmodels和Matplotlib。
步驟：
處理時間序列的基本統計方法是ACF（autocorrelation function，自相關函數）、PACF（partial autocorrelation function，偏自相關函數）和譜密度（ts_timeSeriesFunctions.py文件）：

 1 # time series tools
 2 import statsmodels.api as sm
 3 
 4 # read the data
 5 riverFlows = pd.read_csv(data_folder + 'combined_flow.csv',
 6     index_col=0, parse_dates=[0])
 7 
 8 # autocorrelation function
 9 acf = {}    # to store the results
10 f = {}
11 
12 for col in riverFlows.columns:
13     acf[col] = sm.tsa.stattools.acf(riverFlows[col])
14 
15 # partial autocorrelation function
16 pacf = {}
17 
18 for col in riverFlows.columns:
19     pacf[col] = sm.tsa.stattools.pacf(riverFlows[col])
20 
21 # periodogram (spectral density)
22 sd = {}
23 
24 for col in riverFlows.columns:
25     sd[col] = sm.tsa.stattools.periodogram(riverFlows[col])

原理:爲了節省篇幅，咱們這裏省略了引入和聲明變量的步驟。細節可查看源文件。
首先，以相似的方式讀入（咱們在以前的技巧裏建立的）數據集，不過不指定變量的名字，而是用.read_csv（...）方法從文件的第一行獲取。
而後計算ACF。ACF體現了t時刻的觀測值與t+lag時刻的觀測值的關聯有多強，這裏lag是時間間隔。在本書後續部分（使用ARMA和ARIMA模型預測將來），咱們會使用ACF函數計算ARMA（或ARIMA）模型的MA（moving average，移動平均）。在這裏用Statsmodels模型的.acf（...）方法計算時間序列的ACF值。.acf（...）方法的nlags參數，默認值是40，咱們就用默認值。
在給定了以往延遲的條件下，PCAF能夠當作是對當前觀測值的迴歸。咱們用PACF曲線獲得ARMA（或ARIMA）模型的AR（auto-regressive，自迴歸）。使用Statsmodels模型的.pacf（...）方法計算時間序列的PACF值。
若是對ACF和PACF感興趣的話，你能夠參考https://onlinecourses.science.psu.edu/stat510/node/62。
最後計算週期圖法，這個方法能夠幫咱們找到數據中的基礎頻率，即數據中波峯和波谷的主頻率。
讓咱們來繪製該圖：

 1 # plot the data
 2 fig, ax = plt.subplots(2, 3) # 2 rows and 3 columns
 3 
 4 # set the size of the figure explicitly
 5 fig.set_size_inches(12, 7)
 6 
 7 # plot the charts for American
 8 ax[0, 0].plot(acf['american_flow'], colors[0])
 9 ax[0, 1].plot(pacf['american_flow'],colors[1])
10 ax[0, 2].plot(sd['american_flow'],  colors[2])
11 ax[0, 2].yaxis.tick_right() # shows the numbers on the right
12 
13 # plot the charts for Columbia
14 ax[1, 0].plot(acf['columbia_flow'], colors[0])
15 ax[1, 1].plot(pacf['columbia_flow'],colors[1])
16 ax[1, 2].plot(sd['columbia_flow'],  colors[2])
17 ax[1, 2].yaxis.tick_right()
18 
19 # set titles for columns
20 ax[0, 0].set_title('ACF')
21 ax[0, 1].set_title('PACF')
22 ax[0, 2].set_title('Spectral density')
23 
24 # set titles for rows
25 ax[0, 0].set_ylabel('American')
26 ax[1, 0].set_ylabel('Columbia')

首先，咱們建立了一個圖表，將以兩行三列的方式放6個子圖。咱們顯式地指定圖表的大小。
下面畫圖。ax是放座標軸對象的NumPy數組；在咱們的例子中，這是個二維數組。經過指定對象的座標，咱們指定圖表的位置（0，0是左上角）。
行中的最後一個圖表，咱們經過調用.yaxis.tick_right（）方法設置y軸的數字，展現在圖表的右側。
咱們只在繪製ACF、PACF和譜密度時設定圖表標題。對於行，咱們指定y軸的標籤，以區分不一樣的河流。
繪製出的圖表以下：

Tips:

你會發現ACF圖表中有一個重複的模式。這暗示着咱們的數據（如預期般）是週期性的（年度週期模式）。這也暗示着咱們的處理過程是不平穩的。
平穩過程指的是，方差和聯合機率不隨時間變化（http://www.itl.nist.gov/div898/handbook/pmc/section4/pmc442.htm）。
PACF展現了，時刻t的觀測值強烈依賴於時刻t-1和t-2的觀測值（對更以前的快照依賴較少）。分析譜密度可知基礎頻率大約在29；即，每29個月會重複一個基本的模式。這出乎咱們的意料，咱們本覺得序列重複的間隔會是12的倍數。這多是多種因素的影響：錯誤記錄的觀測值，錯誤的設備，或者某些觀測年發生的基本偏移模式。
更多：以前ACF和PACF的圖表並不容易獲得置信區間的主要誤差：這些誤差讓咱們很容易辨別出時間序列過程的AR和MA。
不過，Statsmodels提供了便捷的函數.plot_acf（...）和.plot_pacf（...）（ts_timeSeries Functions_alternative.py文件）：

 1 # plot the data
 2 fig, ax = plt.subplots(2, 2, sharex=True)
 3 
 4 # set the size of the figure explicitly
 5 fig.set_size_inches(8, 7)
 6 
 7 # plot the charts for american
 8 sm.graphics.tsa.plot_acf(
 9     riverFlows['american_flow'].squeeze(),
10     lags=40, ax=ax[0, 0])
11 
12 sm.graphics.tsa.plot_pacf(
13     riverFlows['american_flow'].squeeze(),
14     lags=40, ax=ax[0, 1])
15 
16 # plot the charts for columbia
17 sm.graphics.tsa.plot_acf(
18     riverFlows['columbia_flow'].squeeze(),
19     lags=40, ax=ax[1, 0])
20 
21 sm.graphics.tsa.plot_pacf(
22     riverFlows['columbia_flow'].squeeze(),
23     lags=40, ax=ax[1, 1])
24 
25 # set titles for rows
26 ax[0, 0].set_ylabel('American')
27 ax[1, 0].set_ylabel('Columbia')

首先建立圖表，可容納4個子圖（放在2×2的網格上）。將sharex參數設置爲True意味着圖表的座標軸對象有相同的時域。咱們也顯式指定圖表大小。
.plot_acf（...）方法和.plot_pacf（...）方法都以用來繪圖的數據做爲第一個參數。.squeeze（）方法將單列的DataFrame轉換爲一個序列對象。這裏，咱們顯式指定時間間隔。ax參數指定了圖表放置的位置。


7.4平滑並轉換觀測值

顧名思義，平滑數據就是移除數據中的噪音，使得圖表更爲平滑。
那你該這麼作嗎？從展現的角度說，固然！從建模的角度看，不是必需的——參考William Briggs在http://wmbriggs.com/post/195/的論述。
準備：需裝好pandas、NumPy和Matplotlib。
步驟：本技巧中，咱們將介紹如何計算移動平均和用對數函數轉換數據（ts_smoothing.py文件）：

1 ma_transform12  = pd.rolling_mean(riverFlows, window=12)
2 ma_transformExp = pd.ewma(riverFlows, span=3)
3 log_transfrom   = riverFlows.apply(np.log)

原理：pandas的.rolling_mean（...）方法計算移動平均數。這個方法對時刻>t<與時刻>（t+window-1）<之間的觀測值計算平均數，並用這個平均數替代時刻>t<的觀測值。你也能夠不取window-1個以前的觀測值，而是指定center參數，這樣會在時刻t兩邊取一樣數目的觀測值。
.ewma（...）方法使用一個指數衰減的函數來平滑數據。span參數控制影響窗口——計算加權平均時還有多少相關的觀測值。
若是要更瞭解這些技巧，能夠學習pandas的文檔：http://pandas.pydata.org/pandas-docs/stable/computation.html#exponentially-weighted-moment-functions。
對數據作對數轉換有利於計算，有時候轉換後也能看出更多模式。和平滑不一樣，這是一個可逆的過程，因此咱們不會損失精確度（只要你不作進一步的轉換，好比，計算對數轉換後的移動平均）。這裏，咱們使用NumPy庫裏的.log方法。它取用時刻t的觀測值，和window-1個以前的觀測值，計算其平均數，並用獲得的值替換時刻t的觀測值。
獲得的圖表展現了技巧將原始的數據集變得多麼不一樣：


MA平滑技巧移除了不少噪音，揭示了數據中的局部趨勢；對美國河而言，水位從開始到1942年左右是上升的，而後直到1950年都是降低的，以後又是有升有降。而對於哥倫比亞河，你能夠看到一開始的降低趨勢，1945年左右才翻轉過來。
指數平滑不像MA這麼自然，它只移除數據中最大的波峯，同時保留其總體形狀。log變換將數據的幅度正規化。前面提到過，這是惟一一個徹底可逆的技巧，不用擔憂數據丟失。
更多：Holt變換是另外一個指數平滑的例子。區別在於它不對smoothing參數取冪（ts_smoothing_alternative.py文件）：

 1 import pandas as pd
 2 import numpy as np
 3 
 4 def holt_transform(column, alpha):
 5     '''
 6         Method to apply Holt transform
 7 
 8         The transform is given as
 9         y(t) = alpha * x(t) + (1-alpha) y(t-1)
10     '''
11     # create an np.array from the column
12     original = np.array(column)
13 
14     # starting point for the transformation
15     transformed = [original[0]]
16 
17     # apply the transform to the rest of the data
18     for i in range(1, len(original)):
19         transformed.append(
20             original[i] * alpha +
21             (1-alpha) * transformed[-1])
22 
23     return transformed

 

平滑後時刻t的值由其觀測值與以前平滑後的值決定；每一個觀測值的影響由alpha參數控制。
咱們先建立一個NumPy數組。變換的起點是初始的觀測值。而後遍歷數組中全部的觀測值，應用處理邏輯。
這裏使用Holt平滑法，alpha設爲0.5：

# transformations
ma_transformHolt = riverFlows.apply(
    lambda col: holt_transform(col, 0.5), axis=0)

axis設爲0時，.apply（...）方法將列逐一傳入holt_transform（...）方法。
咱們也能夠進行差分處理，這樣可快速移除數據的趨勢，使其穩定；差分是計算時刻t及其前導（時刻t-1）觀測值之差的方法：

difference = riverFlows - riverFlows.shift()

這裏，咱們使用了以前在7.2節中解釋過的.shift（...）方法。
獲得的變換以下所示：

真心贊一個！
能夠看出，如同指數平滑同樣，這個方法也移除了大部分噪音，但保留了數據的形狀。若是相比於時間序列中的值自己，你更關心預測變化，那麼差分是頗有用的；預測股市變化就是一個應用場景。

7.5過濾時間序列數據

經過平滑移除噪音只是技巧之一。本技巧中，咱們看看如何經過卷積及其餘濾波器從數據中提取某些頻率。
準備：需裝好pandas、Statsmodels、NumPy、Scipy和Matplotlib。

步驟：
簡單來講，卷積就是f函數（咱們的時間序列）和g函數（咱們的濾波器）的交疊。卷積模糊了時間序列（從這個角度，也能夠將其當作一個平滑技巧）。
這裏有一個不錯的對卷積的介紹：http://www.songho.ca/dsp/convolution/convolution.html。
下面的代碼在ts_filtering.py中：

 1 # prepare different filters準備不一樣的濾波器
 2 MA_filter = [1] * 12
 3 linear_filter = [d * (1 / 12) for d in range(0, 13)]
 4 gaussian = [0] * 6 + list(sc.signal.gaussian(13, 2)[:7])
 5 
 6 # convolve卷積
 7 conv_ma = riverFlows.apply(
 8     lambda col: sm.tsa.filters.convolution_filter(
 9         col, MA_filter), axis=0).dropna()
10 
11 conv_linear = riverFlows.apply(
12     lambda col: sm.tsa.filters.convolution_filter(
13         col, linear_filter), axis=0).dropna()
14 
15 conv_gauss = riverFlows.apply(
16     lambda col: sm.tsa.filters.convolution_filter(
17         col, gaussian), axis=0).dropna()

原理：
第一個濾波器，咱們經過卷積能夠實現一個移動平均（MA）濾波器：這個濾波器是由12個元素組成的列表，每一個元素都是1。
第二個濾波器在12個週期內線性增加；這個濾波器逐步下降輸出值中舊觀測值的重要性。
最後一個濾波器使用高斯函數過濾數據。
這些濾波器的響應看上去像這樣：

咱們使用Statsmodels的.convolution_filter（...）方法過濾數據集。這個方法將數據做爲第一個參數及過濾數組。
對有些觀測值，方法自己並不能爲其生成任何結果，因此咱們使用.dropna（）移除缺失的觀測值。過濾後的數據集以下所示：

經過卷積的MA產生了和以前的.rolling_mean（...）方法一樣的結果。線性濾波器從數據集中移除波峯。最後，高斯模糊不只減小了觀測值的幅度，並且讓其更平滑。
更多：對於經濟數據（或者像咱們的例子，天然中的數據），某些濾波器會更合適：好比BK（Baxter-King）、HP（Hodrick-Prescott）與CF（Christiano-Fitzgerald）。
若是有興趣，能夠查看相關的學術論文：http://www.nber.org/papers/w5022.pdf（Baxter-King），

https://www0.gsb.columbia.edu/faculty/rhodrick/prescott-hodrick1997.pdf（Hodrick-Prescott）

與http://www.nber.org/papers/w7257.pdf（Christiano-Fitzgerald）。
咱們使用Statsmodels實現這些濾波器：

bkfilter = sm.tsa.filters.bkfilter(riverFlows, 18, 96, 12)
hpcycle, hptrend = sm.tsa.filters.hpfilter(riverFlows, 129600)
cfcycle, cftrend = sm.tsa.filters.cffilter(riverFlows,
    18, 96, False) 

BK過濾器是一個帶通濾波器；從時間序列中移除高頻和低頻。.bkfilter（...）方法以數據做爲第一個參數。下一個參數是振動的最小週期；對於月度數據，推薦使用18，季度推薦使用6，年度推薦使用1.5。下一個參數定義了振動的最大週期：對於月度數據，推薦使用96。最後一個參數決定了濾波器的超前-滯後（或者說，一個窗口）。
HP過濾器經過解決一個最小化問題，將初始的時間序列拆成趨勢和業務週期組件。.hpfilter（...）方法以數據做爲第一個參數及平滑的參數。一般，對季度數據使用1600做爲頻率，年度數據使用6.25，月度數據（即本例中）使用129600。
CF過濾器將初始的時間序列拆成趨勢和業務週期組件，在這個意義上來講，和HP過濾器相似。.cffilter（...）也是以數據做爲第一個參數。與BK過濾器相似，接下來兩個參數指定了振動的最小週期和最大週期。最後一個參數drift指定了是否要從數據中移除趨勢。
過濾後的數據以下：

BK過濾器從數據中移除幅度，並使其靜止。分析HP過濾器的輸出能夠看出，哥倫比亞河的長期趨勢幾乎是恆定的，而美國河的趨勢一直在變。美國河的週期組件也體現了相似的模式。CF過濾器的輸出證實了這一點：美國河的趨勢組件比哥倫比亞河的更多變。

7.6移除趨勢和季節性

如同以前提到的，時間序列是平穩的等價條件，而且其平均值、方差和自相關都不隨時間變化。這意味着，有趨勢和季節性的時間過程就是不平穩的。
下一個技巧要介紹的ARMA模型和ARIMA模型要求數據是平穩的（或接近平穩）。因此，本技巧中，咱們學習如何從河水流動數據中移除趨勢和季節性。
準備：需裝好pandas、NumPy、Statsmodels和Matplotlib。

步驟：Statsmodels提供了方便的移除趨勢和季節性的方法（ts_detrendAndRemoveSeasonality.py文件）：

 1 import pandas as pd
 2 import numpy as np
 3 import matplotlib
 4 import matplotlib.pyplot as plt
 5 
 6 # change the font size
 7 matplotlib.rc('xtick', labelsize=9)
 8 matplotlib.rc('ytick', labelsize=9)
 9 matplotlib.rc('font', size=14)
10 
11 # time series tools
12 import statsmodels.api as sm
13 
14 def period_mean(data, freq):
15     '''
16         Method to calculate mean for each frequency
17     '''
18     return np.array(
19         [np.mean(data[i::freq]) for i in range(freq)])
20 
21 # folder with data
22 data_folder = '../../Data/Chapter07/'
23 
24 # colors
25 colors = ['#FF6600', '#000000', '#29407C', '#660000']
26 
27 # read the data
28 # riverFlows = pd.read_csv(data_folder + 'combined_flow.csv',
29 #     index_col=0, parse_dates=[0])
30 riverFlows = pd.read_csv(data_folder + 'combined_flow.csv')
31 
32 
33 # detrend the data
34 detrended = sm.tsa.tsatools.detrend(riverFlows,
35     order=1, axis=0)
36 
37 # create a data frame with the detrended data
38 detrended = pd.DataFrame(detrended, index=riverFlows.index,
39     columns=['american_flow_d', 'columbia_flow_d'])
40 
41 # join to the main dataset
42 riverFlows = riverFlows.join(detrended)
43 
44 # calculate trend
45 riverFlows['american_flow_t'] = riverFlows['american_flow'] \
46     - riverFlows['american_flow_d']
47 riverFlows['columbia_flow_t'] = riverFlows['columbia_flow'] \
48     - riverFlows['columbia_flow_d']
49 
50 # number of observations and frequency of seasonal component
51 nobs = len(riverFlows)
52 freq = 12   # yearly seasonality
53 
54 # remove the seasonality
55 for col in ['american_flow_d', 'columbia_flow_d']:
56     period_averages = period_mean(riverFlows[col], freq)
57     riverFlows[col[:-2]+'_s'] = np.tile(period_averages,
58         nobs // freq + 1)[:nobs]
59     riverFlows[col[:-2]+'_r'] = np.array(riverFlows[col]) \
60         - np.array(riverFlows[col[:-2]+'_s']) 

原理：首先，咱們讀入數據。
使用Statsmodels的.detrend（...）方法，咱們從數據中移除趨勢。order參數指定了趨勢的類型：0表明常數，1表明趨勢是線性的，2表明要移除的是一個二次趨勢。
.detrend（...）方法返回一個NumPy數組；在咱們的例子中，這是一個二維數組，由於初始的數據集中有兩列。因此，爲了便於加入初始數據集，下一步咱們用移除趨勢的數據建立一個DataFrame。咱們重用了初始DataFrame的索引：riverFlows。並且，爲了不和初始數據集衝突，列名加上後綴_d。
pandas的.join（...）方法（默認狀況下）根據索引合併兩個DataFrame：匹配兩個DataFrame的索引，返回對應的值。不過，.join（...）方法給你控制權，容許你指定要鏈接到的列；鍵-列在兩個數據集中都要有。它也容許你指定如何鏈接DataFrame：默認狀況下，這是個左鏈接，從調用方DataFrame的每行記錄出發，鏈接傳入DataFrame中鍵名匹配的全部值；因爲兩個DataFrame中索引相同，咱們只須要調用這個方法便可。
要理解其餘類型的鏈接，參考http://www.sql-join.com。原做者也推薦查看.join（...）方法的文檔：http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.join.html。
獲得的時間序列和初始的並無太大差異：

/* 錯，生成空白！！！！ */

解決方案：

/*
# create a data frame with the detrended data
#用移除趨勢的數據建立一個DateFrame
# detrended = pd.DataFrame(detrended, index=riverFlows.index,
#     columns=['american_flow_d', 'columbia_flow_d'])
#原方法中列'american_flow_d', 'columbia_flow_d'都是NaN
detrended.rename(columns={'american_flow':'american_flow_d', 'columbia_flow':'columbia_flow_d'}, inplace = True)
  */

 

趨勢就是初始值與去趨勢化以後的值之間的差異。
去趨勢化時間序列以後，咱們能夠計算季節性組件。咱們指望獲得的是一個年度季節性，因此freq被設爲12。這意味着咱們指望每12個月會重複一個模式。period_mean（...）方法就是計算這個的。data[i：：freq]選取第i個觀測值，並在第i個以後，每隔12個（freq）選取觀測值，動態建立列表。
子集語法d[f：l：freq]指定了第一個位置f，最後一個位置l，以及取樣的頻率freq。假設d=[1，2，3，4，5，6，7，8，9，10]和d[1：：2]，咱們會獲得[2，4，6，8，10]。
NumPy的.mean（...）方法計算數組中全部元素的平均值。對於每條河流，咱們獲得下面的結果：

美國河看起來更柔順：八月份左右逐漸漲到頂峯，而後越到年底越回落。相反，哥倫比亞河一年中大部分都很平靜，到夏季忽然漲得厲害。
計算了季節性的平均，咱們能夠從去趨勢化後的數據中提取出來。咱們使用NumPy的.tile（...）方法針對時間序列的長度重複季節模式。這個方法以要重複的模式做爲第一個參數。第二個參數指定模式要重複多少次。
//操做符作的是整數除法，將獲得的商圓整到最近的整數。
最後，咱們計算移除季節性後的殘差：

這裏，咱們將時間序列分解成一個線性趨勢、季節性組件以及殘差。能夠看出，兩條河的流量都隨着時間增加，可是增加可忽略。咱們能夠清楚地看出（指望中的）年度模式，哥倫比亞河的方差更大。對於殘差，美國河變化的幅度要顯著高於哥倫比亞河。
更多：
Statsmodels有另外一個分解時間序列的有用方法：.seasonal_decompose（...）。
不幸的是，Statsmodels的文檔是缺失的；項目只在發佈筆記中提到了.seasonal_decompose（...）方法，舉了一個例子，而文檔是找不到的。學習GitHub上的源代碼，能夠發現更多內容，做者推薦你也這麼作（https://github.com/statsmodels/statsmodels/blob/master/statsmodels/tsa/seasonal.py。
這個方法以要分解的數據做爲第一個參數。你能夠指定基礎模型的類型：可加的（咱們的例子）或可乘的（傳入m）。你也能夠指定freq參數；這個參數可不傳，頻率也能經過數據推斷出。
關於可加和可乘的分解模型的區別，推薦閱讀https://onlinecourses.science.psu.edu/stat510/node/69。
方法返回一個DecomposeResult對象。要訪問分解的組件，可調用對象的屬性（ts_seasonalDecomposition.py）：

 1  for col in riverFlows.columns:
 2     # seasonal decomposition of the data
 3     sd = sm.tsa.seasonal_decompose(riverFlows[col],
 4         model='a', freq=12)
 5 
 6     riverFlows[col + '_resid'] = sd.resid \
 7         .fillna(np.mean(sd.resid))
 8 
 9     riverFlows[col + '_trend'] = sd.trend \
10         .fillna(np.mean(sd.trend))
11 
12     riverFlows[col + '_seas'] = sd.seasonal \
13         .fillna(np.mean(sd.seasonal))

獲得的分解和以前獲得的看上去不一樣：

差別來自於移除趨勢的方式不一樣：咱們假設趨勢在整個時域上是線性的，而.seasonal_decompose（...）方法使用卷積濾波以發現基礎趨勢。仔細觀察的話，你會發現這個趨勢圖和過濾時間序列數據中的MA圖的類似之處；區別在於.seasonal_decompose（...）方法使用的權重。

7.7使用ARMA和ARIMA模型預測將來

ARMA（autoregressive moving average，自迴歸移動平均）模型及其泛化——ARIMA（autoregressive integrated moving average，差分自迴歸移動平均）模型——是從時間序列預測將來時經常使用的兩個模型。ARIMA的泛化是指這是一個總體：模型的第一步是在估算AR和MA以前作差分。
準備:需裝好pandas、NumPy、Statsmodels和Matplotlib。你還須要前面技巧裏準備的數據。
步驟：咱們將處理過程包裝在方法裏，以便大部分建模能夠自動化（ts_arima.py文件）：

 1 def plot_functions(data, name):
 2     '''
 3         Method to plot the ACF and PACF functions
 4     '''
 5     # create the figure
 6     fig, ax = plt.subplots(2)
 7 
 8     # plot the functions
 9     sm.graphics.tsa.plot_acf(data, lags=18, ax=ax[0])
10     sm.graphics.tsa.plot_pacf(data, lags=18, ax=ax[1])
11 
12     # set titles for charts
13     ax[0].set_title(name.split('_')[-1])
14     ax[1].set_title('')
15 
16     # set titles for rows
17     ax[0].set_ylabel('ACF')
18     ax[1].set_ylabel('PACF')
19 
20     # save the figure
21     plt.savefig(data_folder+'/charts/'+name+'.png',
22         dpi=300)
23 
24 def fit_model(data, params, modelType, f, t):
25     '''
26         Wrapper method to fit and plot the model
27     '''
28     # create the model object
29     model = sm.tsa.ARIMA(data, params)
30 
31     # fit the model
32     model_res = model.fit(maxiter=600, trend='nc',
33         start_params=[.1] * (params[0]+params[2]), tol=1e-06)
34 
35     # create figure
36     fig, ax = plt.subplots(1, figsize=(12, 8))
37     
38     e = model.geterrors(model_res.params)
39     ax.plot(e, colors[3])
40 
41     chartText = '{0}: ({1}, {2}, {3})'.format(
42         modelType.split('_')[0], params[0],
43         params[1], params[2])
44 
45     # and put it on the chart
46     ax.text(0.1, 0.95, chartText, transform=ax.transAxes)
47 
48     # and save the figure
49     plt.savefig(data_folder+'/charts/'+modelType+'_errors.png',
50         dpi=300)
51 
52 
53     # plot the model
54     plot_model(data['1950':], model_res, params,
55         modelType, f, t)
56 
57     # and save the figure
58     plt.savefig(data_folder+'/charts/'+modelType+'.png',
59         dpi=300)
60 
61 def plot_model(data, model, params, modelType, f, t):
62     '''
63         Method to plot the predictions of the model
64     '''
65     # create figure
66     fig, ax = plt.subplots(1, figsize=(12, 8))
67     
68     # plot the data
69     data.plot(ax=ax, color=colors[0])
70 
71     # plot the forecast
72     model.plot_predict(f, t, ax=ax, plot_insample=False)
73 
74     # define chart text
75     chartText = '{0}: ({1}, {2}, {3})'.format(
76         modelType.split('_')[0], params[0],
77         params[1], params[2])
78 
79     # and put it on the chart
80     ax.text(0.1, 0.95, chartText, transform=ax.transAxes)

原理：讀入數據後，咱們首先查看殘差的ACF和PACF函數：

1 #plot the ACF and PACF functions
2 plot_functions(riverFlows['american_flow_r'],
3     'ACF_PACF_American')
4 plot_functions(riverFlows['columbia_flow_r'],
5     'ACF_PACF_Columbia') 

分析圖表有助於咱們決定模型的AR和MA組件：

前一張圖和下一張圖分別表明美國河和哥倫比亞河的ACF和PACF函數：


plot_functions（...）方法很像咱們在7.3節中寫的代碼，因此這裏咱們不討論它。
查看圖表，咱們能夠看到自迴歸組件和移動平均組件，繼而能夠在模型構建階段用其做爲起始點：ACF函數有助於決定MA的順序，而PACF函數容許咱們決定AR部分。首要規則是查看顯著落在置信區間以外的觀測值。
從圖表中能夠推斷：對於美國河模型咱們以AR（2）和MA（4）開始，而對哥倫比亞河模型，咱們從AR（3）和MA（2）開始：

 1 # fit american models
 2 fit_model(riverFlows['american_flow_r'], (2, 0, 4),
 3     'ARMA_American', '1960-11-30', '1962')
 4 fit_model(riverFlows['american_flow_r'], (2, 1, 4),
 5     'ARIMA_American', '1960-11-30', '1962')
 6 
 7 # fit colum models
 8 fit_model(riverFlows['columbia_flow_r'], (3, 0, 2),
 9     'ARMA_Columbia', '1960-09-30', '1962')
10 fit_model(riverFlows['columbia_flow_r'], (3, 1, 2),
11     'ARIMA_Columbia', '1960-09-30', '1962')

fit_model（...）方法封裝了模型構建須要的全部步驟。它以要估算模型的數據做爲第一個參數。params參數用於定義模型的參數。modelType參數用於在咱們調用plot_model（...）方法時裝飾圖表；f（從）參數和t（到）參數也傳給了plot_model（...）方法。
首先，咱們建立模型對象。咱們選擇.ARIMA（...）模型，由於它能夠特化爲ARMA模型。.ARIMA（...）方法以數據做爲第一個參數，以參數元組做爲第二個參數。元組的第一個元素描述了模型的AR部分，第二個元素描述了（ARIMA模型中）用於差分的滯後，最後一個元素描述了MA組件。
將差分部分設爲0，ARIMA模型就成了ARMA模型。
而後，咱們調整模型。咱們將循環的最大次數設爲300，趨勢不帶常量。咱們也定義start_params；以一個六元素的列表，全部AR與MR組件的和開始。列表中每一個元素都設成了起始點0.1。容忍度設爲0.000001；若是循環間偏差的差異小於容忍度，就中止估算。
估算出模型後，咱們看一下它是如何完成的。plot_model（...）方法繪製觀測的數據。從圖表的簡明考慮，咱們調用方法時限制數據從1950開始。.plot_predict（...）方法使用估算的模型參數預測將來的觀測值。頭兩個參數指定預測的上下界限：咱們能夠選擇起始點與結束點。咱們將plot_insample設爲False；這個參數的文檔極爲缺少，咱們就經驗主義一把，認爲這個參數避免了方法用不一樣的顏色畫與預測交疊的觀測值。
不要預測太遠的將來，由於預測的置信區間沒理由就會變寬。
咱們在圖表的左上角加上文字；用.text（...）方法作到這一點。頭兩個參數是圖表上的座標：用.transformAxes轉換軸座標，咱們知道座標（0，0）是左下角，（1，1）是右上角。第三個參數是但願加的文字。
估算出了模型，看下圖表：

預測的美國河流量一開始遵循一個短時間趨勢，可是接下來（因爲咱們要從最後一個觀測值往前走更遠），置信區間激增，預測值就成了一個常數：


對於哥倫比亞河，一開始預測的流量看上去偏離了觀測值，可是很快就平坦了：

ARIMA模型的預測和ARMA模型的相差不大；預測一開始遵循觀測數據，後來平坦了：

模型本質上是預測了殘差的均值。
查看John Wittenauer對標普500指數的預測，可明白預測時間序列不是一件小事（http://www.johnwittenauer.net/a-simple-time-series-analysis-of-the-sp-500-index/）。
ACF和PACF函數獲得的AR和MA參數應該只用做起始點。若是模型表現很好，就保留。不然，你能夠微調。
咱們也應該查看模型偏差項的分佈：


美國河的ARMA模型的偏差項看上去很隨機，這是預測時間序列時咱們要追求的：


對於哥倫比亞河，偏差項能看到一些重複的模式，這給模型預測將來流量的能力打上了一個問號：


美國河的ARIMA模型也生成了相似隨機殘差的結果：


最終，偏差項看上去應該像是白噪音。對於哥倫比亞河模型，這不必定對，由於咱們也看到出現了模式。你能夠試試不一樣的AR和MA參數，看模型是否有不一樣的結果。
參考

McKinney、Perktold和Seabold對時間序列作出了一個很好的展現/教程。值得一看：http://conference.scipy.org/scipy2011/slides/mckinney_time_series.pdf。

 

第7章完。

 python學習筆記-目錄索引

 

隨書源碼官方下載：
http://www.hzcourse.com/web/refbook/detail/7821/92