時間序列（time series）是來自隨時間變化的系統的一系列度量。

本章使用的示例來自Zachary M. Jones。Jones的研究目的是調查像大麻合法化這樣的政策性決定會對市場產生何種影響。希望大家對本章內容感興趣，但藉此機會重申對資料分析保持專業性態度的重要性。藥品是否非法，哪些藥品應當屬於非法，這是很重要而又難以回答的公共政策問題，人們應當基於誠實準確的資料進行決策。

匯入和清洗資料

從Jones先生的網站下載資料，通過下面程式碼將資料讀取到一個pandas DataFrame中：

transactions = pandas.read_csv('mj-clean.csv', parse_dates=[
 
5])

這個DataFrame中每一行都代表一次交易，具有如下列：

• city 字串，代表城市名
• state 州名的兩字母縮寫
• price 交易價格，單位為美元
• amount 購買量，單位為克
• quality 由購買者彙報的貨品質量，值為高、中或低
• date 彙報日期，這個日期應該在購買日期後不久
• ppg 每克價格，單位為美元
• lat 交易發生地點的大致維度，由城市名獲得
• lon 交易發生地點的大致經度

每次交易都是一個時間事件，因此這個資料集可以看成一個事件序列。很多分析時間序列的方法要求度量是均勻分佈的，或者度量均勻分佈至少可以使分析更為簡單。

為了演示這些方法，將這個資料集按照彙報的質量分組，並計算每克的日均價格，將每組資料轉換為均勻分佈的序列。

def GroupByQualityAndDay(transactions):
    groups = transactions.groupby('quality')
    dailies = {}
    for name, group in groups:
        dailies[name] = GroupByDay(group)
    
    return dailies
    
def GroupByDay(transactions, func=np.mean):
    grouped = transactions[['date', 'ppg']].groupby('date'
 
)
    daily = grouped.aggregate(func)
    
    daily['date'] = daily.index
    start = daily.date[0]
    one_year = np.timedelta64(1, 'Y')
    daily['years'] = (daily.date - start) / one_year
    
    return daily

groupby是DataFrame的方法，返回一個GroupBy物件groups。得到的grouped是一個對映，將每個日期對應到包含該日期彙報價格的DataFrame。aggregate是GroupBy的方法，遍歷所有的組，對組中每列都執行一個函式。GroupByDay方法得到的DataFrame包含列ppg、date和years。

繪製圖形

GroupByQualityAndDay得到的結果是從每個質量級別到日均價格DataFrame的對映。繪製這3個時間序列的程式碼如下：

thinkplot.PrePlot(rows=3)
for i, (name, daily) in enumerate(dailies.items()):
    thinkplot.SubPlot(i+1)
    title = 'price per gram ($)' if i==0 else ''
    thinkplot.Config(ylim=[0, 20], title=title)
    thinkplot.Scatter(daily.index, daily.ppg, s=10, label=name)
    if i==2:
        pyplot.xticks(rotation=30)
    else:
        thinkplot.Config(xticks=[])

具體繪製圖形如下：
圖1 高、中、低質量大麻的每克每日價格時間序列

從圖中可以看到，這段時間內，高質量大麻價格似乎在降低；中等質量大麻價格在提高；低質量大麻的價格也在提高，但其價格變化過大，很難判斷趨勢。大麻質量的資料是由志願者提供的，因此，這些隨時間變化的趨勢可能反映了參與者判斷標準的變化。

線性迴歸

雖然時間序列分析有專門的方法，但對於很多問題，簡單的方法是使用通用的工具，如線性迴歸。下面的函式以每日價格的DataFrame為引數，計算一個最小二乘法擬合，返回使用StatsModels包得到的模型和結果物件。然後，遍歷dailies中的DataFrame，為每個DataFrame擬合模型。

def RunLinearModel(daily):
    model = smf.ols('ppg ~ years', data=daily)
    results = model.fit()
    return model, results
    
for name, daily in dailies.items():
    model, results = RunLinearModel(daily)
    print(name)
    regression.SummerizeResults(results)

結果如下：

估計所得斜率說明，在觀測區間內，高質量大麻的價格每年下降約71%；中等質量大麻的價格每年上漲約28%；低質量大麻價格每年上漲57%。這些估計值都是統計顯著的，p值很小。

下圖繪製高質量大麻資料點的散點圖及擬合值的線圖。

圖2 高質量大麻每克每日價格的時間序列及線性最小二乘法擬合

圖中模型似乎很好地擬合了資料，然而，線性迴歸並不是擬合這種資料的最適宜方法，原因如下：

• 預期長期趨勢為線性的或符合其他簡單函式，這是沒有依據的。通常，價格由供需關係決定，而供給和需求都會隨著時間發生不可預期的變化。
• 線性迴歸模型對所有資料使用通用的權重，不管是過去的還是近期資料。在進行預測是，近期資料也許應該得到更多的權重。
• 線性迴歸假設殘差是無關的噪音。而在時間序列資料中，相鄰資料是相關的，這項假設通常不能成立。

移動平均值

大部分時間序列分析基於的模型假設都認為，觀察序列是三部分的總和：

• 趨勢 描述持續變化的一個平滑函式
• 季節性 週期性的變化，可能包括每日、每週、每月或每年週期
• 噪音 長期趨勢周圍的隨機變化

迴歸是從一個序列中獲取趨勢的方法。但是，如果這個趨勢不是一個簡單函式，那麼一個很好的方法是移動平均值（moving average）。移動平均值將序列分為相互重疊的區域（稱為視窗），計算每個視窗的平均值。

最簡單的移動平均值是滾動均值（rolling mean），滾動均值計算每個視窗中值的均值。pandas提供rolling_mean方法，引數為Series和視窗大小，返回一個新的Series。

在對大麻資料使用rolling_mean方法之前，必須首先處理缺失資料。目前使用的DataFrame缺失這些資料，索引跳過了沒有資料的日期。為了隨後的分析，需要明確標明這些缺失的資料，為此，可以對DataFrame進行“重建索引”。

下面程式碼進行索引重建和滾動均值的繪製。

dates = pandas.date_range(daily.index.minx(), daily.index.max())
reindexed = daily.reindex(dates)

roll_mean = pandas.rolling_mean(reindexed.ppg, 30)
thinkplot.Plot(roll_mean.index, roll_mean)

另一種移動平均值是指數權重移動平均（exponentially-weighted moving average, EWMA)。它有兩個優點：指數權重移動平均計算加權平均值，最近的值具有最高的權重；pandas的EWMA實現可以較好地處理缺失值。

下圖展示了兩種平均值計算的效果：

圖3 每日價格、滾動均值（左）和指數權重移動均值（右）

缺失值

填充缺失資料，一個簡單的方法是使用移動平均值。Series的方法fillna就實現了這一功能。這一方法的缺點是弱化了序列中的噪音，這一問題可以通過新增重抽樣的殘差解決。

resid = (reindexed.ppg - ewma).dropna()
fake_data = ewma + Resample(resid, len(reindexed))
reindexed.ppg.fillna(fake_data, inplace=True)

resid包含殘差，但不包含ppg為nan的日期。fake_data包含移動平均值的結果與殘差的一個隨機樣本的和。最後，fillna將nan替換為fake_data中的值。

下圖展示了填充資料後的結果。

圖4 填充資料後的每日價格

序列相關

隨著價格每日變化，你可能會期望看到一些模式。如果週一價格很高，可能隨後幾天都會較高；如果價格很低，可能會保持低位。在這種模式中，每個值都與序列中的下一個值相關，因此稱為序列相關（serial correlation）。

要計算序列相關，我們可以將時間序列移動一個稱為滯後（lag）間隔，然後計算移動後的序列與原序列的相關性。

def SerialCorr(series, lag=1):
    xs = series[lag:]
    ys = series.shift(lag)[lag:]
    corr = Corr(xs, ys)
    return corr

使用不同的滯後值，可以檢驗序列的每週、每月和每年的季節性特徵，具體結果如下：

下一節將檢驗這些相關是否統計顯著（結果不是），目前暫時認為，這些序列沒有顯著的季節性模式，至少在使用上述滯後值時沒有。

自相關

自相關函式（autocorrelaton funciton）將滯後值對映到使用該值得到的序列相關。“自相關”是序列相關的另一個名字，常用於滯後值不為1時。

StatsModels軟體包提供了時間序列分析函式，其中有計算自相關函式的acf。acf使用從0到nlags的滯後值計算序列相關。引數unbiased為True，告訴acf要為樣本規模校正估計值。

如果選擇高質量大麻的每日價格，抽取滯後為1、7、30和365的相關值，可以驗證acf和SerialCorr得到的結果大致相同。

import statsmodels.tsa.stattools as smtsa
acf = smtsa.acf(filled.resid, nlags=365, unbiased=True)

acf[0], acf[1], acf[7], acf[30], acf[365]
1.000, -0.029, 0.020, 0.014, 0.044

下圖（左）展示了nlags=40時，3種質量分類的自相關函式。圖中的灰色區域是不存在自相關是的正態可變性，位於這個區域之外的都是統計顯著的。

圖5 每日價格的自相關函式（左）及模擬的每週季節性的每日價格（右）

為了展示存在季節性因素的自相關函式，在資料中加入一個每週的迴圈進行模擬。假設在週末是大麻的需求量較大，那麼價格可能會較高。為了模擬這個效果，選擇週五或週六的日期，在價格上新增一個隨機量，這個隨機量選自從0~2美分的均勻分佈。

上圖（右）展示出添加了模擬季節性的價格自相關函式。如我們所預期的，當滯後為7的倍數時，相關性較高。對於高質量和中等質量大麻，新得到的相關是統計顯著的。而低質量大麻則不然，因為這一分類的殘差最大，必須使用更大的模擬值。

預測

前面提到的線性迴歸可以用於預測，RegressionResults類提供predic方法，以包含解釋變數的DataFrame為引數，返回一個預測序列。如果我們希望得到的只是單一的、最佳推測預測，那麼predict就可以滿足需求。但大部分時候需要對誤差進行量化，即希望得知預測結果的準確性如何。

我們需要考慮三種誤差來源：

• 抽樣誤差 預測基於估計引數，而估計引數依賴樣本中的隨機變異。
• 隨機變異 即使估計引數是完美的，觀測資料也會在長期趨勢附近隨機變動。
• 建模誤差

另一種誤差來源於無法預期的未來事件。農產品價格受天氣影響，所有價格都會受政策和法律影響。

建模誤差和無法預期的未來事件很難量化，而抽樣誤差和隨機變異較容易處理，因此我們先處理後兩種誤差。

我們依然使用重抽樣方法對抽樣誤差進行量化，重抽樣的目的是使用實際觀測來模擬重複進行實驗時可能得到的資料。這種模擬基於的假設是，估計引數是正確的，但隨機殘差可能會不同。實現這種模擬的函式如下：

def SimulateResults(daily, iters=101):
    model, results = RunLinearModel(daily)
    fake = daily.copy()
    
    result_seq = []
    for i in range(iters):
        fake.pgg = results.fittedvalues + Resample(results.resid)
        _, fake_results = RunLinearModel(fake)
        result_seq.append(fake_results)
        
    return result_seq

在每次迴圈中，程式碼對殘差進行重抽樣，將其附加在擬合值上，生成一個“偽”資料集，然後對偽資料擬合一個線性模型，將結果儲存在RegressionResults物件中。

下一步是使用模擬結果生成預測：

def GeneragePredictions(result_seq, years, add_resid=False):
    n = len(years)
    d = dict(Intercept=np.ones(n), years=years, years2=years**2)
    predict_df = pandas.DataFrame(d)
    
    predict_seq = []
    for fake_results in result_seq:
        predict = fake_results.predict(predict_df)
        if add_resid:
            predict += Resample(fake_results.resid, n)
        predict_seq.append(predict)
        
    return predict_seq

最後一步是繪製預測結果的90%置信區間,下圖展示程式碼執行結果。圖中深灰色區域代表取樣誤差的90%置信區間，即由取樣導致的估計斜率和截距不確定性。

圖6 基於線性擬合的預測，展示了由抽樣誤差和預測誤差導致的變異

迴歸模型基於的假設是，系統是平穩的（stationary），即模型引數不會隨著時間發生變化。具體來說，迴歸模型假設模型的斜率和截距是常數，殘差分佈也不變。

但前面的移動平均值圖中，斜率在觀測區間中似乎至少變化了一次，而且前半部分的殘差方差似乎也比後半部分大。因此，我們得到的引數依賴觀測區間。為了檢驗這一現象對預測結果的影響，我們可以對SimulateResult進行擴充套件，使用具有不同開始和結束日期的觀測區間進行擬合。

下圖展示了中等質量分類的預測結果。圖中淺灰色區域代表各種誤差的置信區間，包括取樣誤差、隨機變異和觀測區間變化。

圖7 基於線性擬合的預測，展示了由觀測間隔導致的變異

基於整個區間的模型的斜率為正數，說明價格在上漲。而最近的區間顯示出價格有下跌跡象，因此基於最近資料的模型斜率為負數。

參考文獻：
    統計思維. Allen B.Downey. 金迎 譯