托馬茲·卓巴斯的《資料分析實戰》，2018年6月出版，本系列為讀書筆記。主要是為了系統整理，加深記憶。 第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軸的标簽，以區分不同的河流。
繪制出的圖表如下：
你會發現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/的論述。
步驟：本技巧中，我們将介紹如何計算移動*均和用對數函數轉換資料（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章完。
随書源碼官方下載下傳：
http://www.hzcourse.com/web/refbook/detail/7821/92
邀月注：本文版權由邀月和部落格園共同所有，轉載請注明出處。
助人等于自助! [email protected]