機器學習開放課程:九、基於Python分析真實手遊時序數據

編者按：Zeptolab數據科學家Dmitriy Sergeev介紹了分析和預測時序數據的主要方法。

大家好！

這次的開放機器學習課程的內容是時序數據。

我們將查看如何使用Python處理時序數據，哪些方法和模型可以用來預測；什麼是雙指數平滑和三指數平滑；如果平穩（stationarity）不是你的菜，該怎麼辦；如何創建SARIMA並且活下來；如何使用XGBoost做出預測。所有這些都將應用於（嚴酷的）真實世界例子。

導言

平移、平滑、評估

滾動窗口估計

指數平滑，Holt-Winters模型

時序交叉驗證，參數選擇

計量經濟學方法

時序數據上的線性（以及不那麼線性）模型

特徵提取

線性模型，特徵重要性

正則化，特徵選取

XGBoost

相關資源

在我的工作中，我幾乎每天都會碰到和時序有關的任務。最頻繁的問題是——明天/下一周/下個月/等等，我們的指標將是什麼樣——有多少玩家會安裝應用，他們的在線時長會是多少，他們會進行多少次操作，取決於預測所需的質量，預測周期的長度，以及時刻，我們需要選擇特徵，調整參數，以取得所需結果。

時序的簡單定義：

時序——一系列以時間順序為索引（或列出、繪出）的數據點。

因此，數據以相對確定的時刻組織。所以，和隨機樣本相比，可能包含我們將嘗試提取的額外信息。

讓我們導入一些庫。首先我們需要statsmodels庫，它包含了一大堆統計學建模函數，包括時序。對不得不遷移到Python的R粉來說，絕對會感到statsmodels很熟悉，因為它支持類似Wage ~ Age + Education這樣的模型定義。

import numpy as np                               # 向量和矩陣

import pandas as pd                              # 表格和數據處理

import matplotlib.pyplot as plt                  # 繪圖

import seaborn as sns                            # 更多繪圖功能

from dateutil.relativedelta import relativedelta # 處理不同格式的時間日期

from scipy.optimize import minimize              # 最小化函數

import statsmodels.formula.api as smf            # 統計學和計量經濟學

import statsmodels.tsa.api as smt

import statsmodels.api as sm

import scipy.stats as scs

from itertools import product                    # 一些有用的函數

from tqdm import tqdm_notebook

import warnings                                  # 勿擾模式

warnings.filterwarnings('ignore')

%matplotlib inline

作為例子，讓我們使用一些真實的手遊數據，玩家每小時觀看的廣告，以及每天遊戲幣的花費：

ads = pd.read_csv('../../data/ads.csv', index_col=['Time'], parse_dates=['Time'])

currency = pd.read_csv('../../data/currency.csv', index_col=['Time'], parse_dates=['Time'])

plt.figure(figsize=(15, 7))

plt.plot(ads.Ads)

plt.title('Ads watched (hourly data)')

plt.grid(True)

plt.show()

plt.figure(figsize=(15, 7))

plt.plot(currency.GEMS_GEMS_SPENT)

plt.title('In-game currency spent (daily data)')

plt.grid(True)

plt.show()

預測質量指標

在實際開始預測之前，先讓我們理解下如何衡量預測的質量，查看下最常見、使用最廣泛的測度：

R2，決定係數（在經濟學中，可以理解為模型能夠解釋的方差比例），(-inf, 1]sklearn.metrics.r2_score

平均絕對誤差（Mean Absolute Error），這是一個易於解釋的測度，因為它的計量單位和初始序列相同，[0, +inf)sklearn.metrics.mean_absolute_error

中位絕對誤差（Median Absolute Error），同樣是一個易於解釋的測度，對離群值的魯棒性很好，[0, +inf)sklearn.metrics.median_absolute_error

均方誤差（Mean Squared Error），最常用的測度，給較大的錯誤更高的懲罰，[0, +inf)sklearn.metrics.mean_squared_error

均方對數誤差（Mean Squared Logarithmic Error），和MSE差不多，只不過先對序列取對數，因此能夠照顧到較小的錯誤，通常用於具有指數趨勢的數據，[0, +inf)sklearn.metrics.mean_squared_log_error

平均絕對百分誤差（Mean Absolute Percentage Error），類似MAE不過基於百分比——當你需要向管理層解釋模型的質量時很方便——[0, +inf)，sklearn中沒有實現。

# 引入上面提到的所有測度

from sklearn.metrics import r2_score, median_absolute_error, mean_absolute_error

from sklearn.metrics import median_absolute_error, mean_squared_error, mean_squared_log_error

# 自行實現sklearn沒有提供的平均絕對百分誤差很容易

def mean_absolute_percentage_error(y_true, y_pred):

   return np.mean(np.abs((y_true - y_pred) / y_true)) * 100

棒極了，現在我們知道如何測量預測的質量了，可以使用哪些測度，以及如何向老闆翻譯結果。剩下的就是創建模型了。

讓我們從一個樸素的假設開始——「明天會和今天一樣」，但是我們並不使用類似y^t=y(t-1)這樣的模型（這其實是一個適用於任意時序預測問題的很好的基線，有時任何模型都無法戰勝這一模型），相反，我們將假定變量未來的值取決於前n個值的平均，所以我們將使用的是移動平均（moving average）。

def moving_average(series, n):

   """

       計算前n項觀測的平均數

   """

   return np.average(series[-n:])

# 根據前24小時的數據預測

moving_average(ads, 24)

結果：116805.0

不幸的是這樣我們無法做出長期預測——為了預測下一步的數據我們需要實際觀測的之前的數據。不過移動平均還有一種用途——平滑原時序以顯示趨勢。pandas提供了實現DataFrame.rolling(window).mean()。窗口越寬，趨勢就越平滑。遇到噪聲很大的數據時（財經數據十分常見），這一過程有助於偵測常見模式。

def plotMovingAverage(series, window, plot_intervals=False, scale=1.96, plot_anomalies=False):

   """

       series - 時序dateframe

       window - 滑窗大小

       plot_intervals - 顯示置信區間

       plot_anomalies - 顯示異常值

   """

   rolling_mean = series.rolling(window=window).mean()

   plt.figure(figsize=(15,5))

   plt.title("Moving average\n window size = {}".format(window))

   plt.plot(rolling_mean, "g", label="Rolling mean trend")

   # 繪製平滑後的數據的置信區間

   if plot_intervals:

       mae = mean_absolute_error(series[window:], rolling_mean[window:])

       deviation = np.std(series[window:] - rolling_mean[window:])

       lower_bond = rolling_mean - (mae + scale * deviation)

       upper_bond = rolling_mean + (mae + scale * deviation)

       plt.plot(upper_bond, "r--", label="Upper Bond / Lower Bond")

       plt.plot(lower_bond, "r--")

       # 得到區間後，找出異常值

       if plot_anomalies:

           anomalies = pd.DataFrame(index=series.index, columns=series.columns)

           anomalies[series<lower_bond] = series[series<lower_bond]

           anomalies[series>upper_bond] = series[series>upper_bond]

           plt.plot(anomalies, "ro", markersize=10)

   plt.plot(series[window:], label="Actual values")

   plt.legend(loc="upper left")

plt.grid(True)

平滑（窗口大小為4小時）：

plotMovingAverage(ads, 4)

平滑（窗口大小為12小時）：

plotMovingAverage(ads, 12)

平滑（窗口大小為24小時）：

plotMovingAverage(ads, 24)

如你所見，在小時數據上按日平滑讓我們可以清楚地看到瀏覽廣告的趨勢。周末數值較高（周末是娛樂時間），工作日一般數值較低。

我們可以同時繪製平滑值的置信區間：

plotMovingAverage(ads, 4, plot_intervals=True)

現在讓我們在移動平均的幫助下創建一個簡單的異常檢測系統。不幸的是，在這段時序數據中，一切都比較正常，所以讓我們故意弄出點異常來：

ads_anomaly = ads.copy()

# 例如廣告瀏覽量下降了20%

ads_anomaly.iloc[-20] = ads_anomaly.iloc[-20] * 0.2

讓我們看看這個簡單的方法能不能捕獲異常：

plotMovingAverage(ads_anomaly, 4, plot_intervals=True, plot_anomalies=True)

酷！按周平滑呢？

plotMovingAverage(currency, 7, plot_intervals=True, plot_anomalies=True)

不好！這是簡單方法的缺陷——它沒能捕捉月度數據的季節性，幾乎將所有30天出現一次的峰值當作異常值。如果你不想有這麼多虛假警報，最好考慮更複雜的模型。

順便提下移動平均的一個簡單修正——加權平均（weighted average）。其中不同的觀測具有不同的權重，所有權重之和為一。通常最近的觀測具有較高的權重。

def weighted_average(series, weights):

   result = 0.0

   weights.reverse()

   for n in range(len(weights)):

       result += series.iloc[-n-1] * weights[n]

   return float(result)

weighted_average(ads, [0.6, 0.3, 0.1])

結果：98423.0

指數平滑

那麼，如果我們不只加權最近的幾項觀測，而是加權全部現有的觀測，但對歷史數據的權重應用指數下降呢？

這一模型的值是當前觀測和歷史觀測的加權平均。權重α稱為平滑因子，定義多快「遺忘」之前的觀測。α越小，之前的值的影響就越大，序列就越平滑。

指數隱藏在函數的遞歸調用之中，y^t-1本身包含(1-α)y^t-2，以此類推，直到序列的開始。

def exponential_smoothing(series, alpha):

   """

       series - 時序數據集

       alpha - 浮點數，範圍[0.0, 1.0]，平滑參數

   """

   result = [series[0]] # 第一項和序列第一項相同

   for n in range(1, len(series)):

       result.append(alpha * series[n] + (1 - alpha) * result[n-1])

   return result

def plotExponentialSmoothing(series, alphas):

   with plt.style.context('seaborn-white'):    

       plt.figure(figsize=(15, 7))

       for alpha in alphas:

           plt.plot(exponential_smoothing(series, alpha), label="Alpha {}".format(alpha))

       plt.plot(series.values, "c", label = "Actual")

       plt.legend(loc="best")

       plt.axis('tight')

       plt.title("Exponential Smoothing")

       plt.grid(True);

plotExponentialSmoothing(ads.Ads, [0.3, 0.05])

plotExponentialSmoothing(currency.GEMS_GEMS_SPENT, [0.3, 0.05])

雙指數平滑

目前我們的方法能給出的只是單個未來數據點的預測（以及一些良好的平滑），這很酷，但還不夠，所以讓我們擴展下指數平滑以預測兩個未來數據點（當然，同樣經過平滑）。

序列分解應該能幫到我們——我們得到兩個分量：截距（也叫水平）l和趨勢（也叫斜率）b。我們使用之前提到的方法學習預測截距（或期望的序列值），並將同樣的指數平滑應用於趨勢（假定時序未來改變的方向取決於之前的加權變化）。

上面的第一個函數描述截距，和之前一樣，它取決於序列的當前值，只不過第二項現在分成水平和趨勢兩個分量。第二個函數描述趨勢，它取決於當前一步的水平變動，以及之前的趨勢值。這裡β係數是指數平滑的權重。最後的預測為模型對截距和趨勢的預測之和。

def double_exponential_smoothing(series, alpha, beta):

   result = [series[0]]

   for n in range(1, len(series)+1):

       if n == 1:

           level, trend = series[0], series[1] - series[0]

       if n >= len(series):

           value = result[-1]

       else:

           value = series[n]

       last_level, level = level, alpha*value + (1-alpha)*(level+trend)

       trend = beta*(level-last_level) + (1-beta)*trend

       result.append(level+trend)

return result

def plotDoubleExponentialSmoothing(series, alphas, betas):

   with plt.style.context('seaborn-white'):    

       plt.figure(figsize=(20, 8))

       for alpha in alphas:

           for beta in betas:

               plt.plot(double_exponential_smoothing(series, alpha, beta), label="Alpha {}, beta {}".format(alpha, beta))

       plt.plot(series.values, label = "Actual")

       plt.legend(loc="best")

       plt.axis('tight')

       plt.title("Double Exponential Smoothing")

       plt.grid(True)

 plotDoubleExponentialSmoothing(ads.Ads, alphas=[0.9, 0.02], betas=[0.9, 0.02])

plotDoubleExponentialSmoothing(currency.GEMS_GEMS_SPENT, alphas=[0.9, 0.02], betas=[0.9, 0.02])

現在我們有兩個可供調節的參數——α和β。前者根據趨勢平滑序列，後者平滑趨勢本身。這兩個參數越大，最新的觀測的權重就越高，建模的序列就越不平滑。這兩個參數的組合可能產生非常怪異的結果，特別是手工設置時。我們很快將查看自動選擇參數的方法，在介紹三次指數平滑之後。

Holt-Winters模型

好哇！讓我們看下一個指數平滑的變體，這次是三次指數平滑。

這一方法的思路是我們加入第三個分量——季節性。這意味著，如果我們的時序不具有季節性（我們之前的例子就不具季節性），我們就不應該使用這一方法。模型中的季節分量將根據季節長度解釋截距和趨勢上的重複波動，季節長度也就是波動重複的周期。季節中的每項觀測有一個單獨的分量，例如，如果季節長度為7（按周計的季節），我們將有7個季節分量，每個季節分量對應一天。

現在我們得到了一個新系統：

現在，截斷取決於時序的當前值減去相應的季節分量，趨勢沒有變動，季節分量取決於時序的當前值減去截斷，以及前一個季節分量的值。注意分量在所有現有的季節上平滑，例如，周一分量會和其他所有周一平均。關於如何計算平均以及趨勢分量和季節分量的初始逼近，可以參考工程統計手冊6.4.3.5：

https://www.itl.nist.gov/div898/handbook/pmc/section4/pmc435.htm

具備季節分量後，我們可以預測任意m步未來，而不是一步或兩步，非常鼓舞人心。

下面是三次指數平滑模型的代碼，也稱Holt-Winters模型，得名於發明人的姓氏——Charles Holt和他的學生Peter Winters。此外，模型中還引入了Brutlag方法，以創建置信區間：

其中T為季節的長度，d為預測偏差。你可以參考以下論文了解這一方法的更多內容，以及它在時序異常檢測中的應用：

https://annals-csis.org/proceedings/2012/pliks/118.pdf

class HoltWinters:

   """

   Holt-Winters模型，使用Brutlag方法檢測異常

   # series - 初始時序

   # slen - 季節長度

   # alpha, beta, gamma - Holt-Winters模型參數

   # n_preds - 預測視野

   # scaling_factor - 設置Brutlag方法的置信區間（通常位於2到3之間）

   """

   def __init__(self, series, slen, alpha, beta, gamma, n_preds, scaling_factor=1.96):

       self.series = series

       self.slen = slen

       self.alpha = alpha

       self.beta = beta

       self.gamma = gamma

       self.n_preds = n_preds

       self.scaling_factor = scaling_factor

   def initial_trend(self):

       sum = 0.0

       for i in range(self.slen):

           sum += float(self.series[i+self.slen] - self.series[i]) / self.slen

       return sum / self.slen  

   def initial_seasonal_components(self):

       seasonals = {}

       season_averages = []

       n_seasons = int(len(self.series)/self.slen)

       # 計算季節平均

       for j in range(n_seasons):

           season_averages.append(sum(self.series[self.slen*j:self.slen*j+self.slen])/float(self.slen))

       # 計算初始值

       for i in range(self.slen):

           sum_of_vals_over_avg = 0.0

           for j in range(n_seasons):

               sum_of_vals_over_avg += self.series[self.slen*j+i]-season_averages[j]

           seasonals[i] = sum_of_vals_over_avg/n_seasons

       return seasonals  

   def triple_exponential_smoothing(self):

       self.result = []

       self.Smooth = []

       self.Season = []

       self.Trend = []

       self.PredictedDeviation = []

       self.UpperBond = []

       self.LowerBond = []

       seasonals = self.initial_seasonal_components()

       for i in range(len(self.series)+self.n_preds):

           if i == 0: # 成分初始化

               smooth = self.series[0]

               trend = self.initial_trend()

               self.result.append(self.series[0])

               self.Smooth.append(smooth)

               self.Trend.append(trend)

               self.Season.append(seasonals[i%self.slen])

               self.PredictedDeviation.append(0)

               self.UpperBond.append(self.result[0] +

                                     self.scaling_factor *

                                     self.PredictedDeviation[0])

               self.LowerBond.append(self.result[0] -

                                     self.scaling_factor *

                                     self.PredictedDeviation[0])

               continue

           if i >= len(self.series): # 預測

               m = i - len(self.series) + 1

               self.result.append((smooth + m*trend) + seasonals[i%self.slen])

               # 預測時在每一步增加不確定性

               self.PredictedDeviation.append(self.PredictedDeviation[-1]*1.01)

           else:

               val = self.series[i]

               last_smooth, smooth = smooth, self.alpha*(val-seasonals[i%self.slen]) + (1-self.alpha)*(smooth+trend)

               trend = self.beta * (smooth-last_smooth) + (1-self.beta)*trend

               seasonals[i%self.slen] = self.gamma*(val-smooth) + (1-self.gamma)*seasonals[i%self.slen]

               self.result.append(smooth+trend+seasonals[i%self.slen])

               # 據Brutlag算法計算偏差

               self.PredictedDeviation.append(self.gamma * np.abs(self.series[i] - self.result[i])

                                              + (1-self.gamma)*self.PredictedDeviation[-1])

           self.UpperBond.append(self.result[-1] +

                                 self.scaling_factor *

                                 self.PredictedDeviation[-1])

           self.LowerBond.append(self.result[-1] -

                                 self.scaling_factor *

                                 self.PredictedDeviation[-1])

           self.Smooth.append(smooth)

           self.Trend.append(trend)

self.Season.append(seasonals[i%self.slen])

時序交叉驗證

現在我們該兌現之前的承諾，討論下如何自動估計模型參數。

這裡沒什麼不同尋常的，我們需要選擇一個適合任務的損失函數，以了解模型逼近數據的程度。接著，我們通過交叉驗證為給定的模型參數評估選擇的交叉函數，計算梯度，調整模型參數，等等，勇敢地下降到誤差的全局最小值。

問題在於如何在時序數據上進行交叉驗證，因為，你知道，時序數據確實具有時間結構，不能在一折中隨機混合數據（不保留時間結構），否則觀測間的所有時間相關性都會丟失。這就是為什麼我們將使用技巧性更強的方法來優化模型參數的原因。我不知道這個方法是否有正式的名稱，但是在CrossValidated上（在這個網站上你可以找到所有問題的答案，生命、宇宙以及任何事情的終極答案除外），有人提出了「滾動式交叉驗證」（cross-validation on a rolling basis）這一名稱。

這一想法很簡單——我們在時序數據的一小段上訓練模型，從時序開始到某一時刻t，預測接下來的t+n步並計算誤差。接著擴張訓練樣本至t+n個值，並預測從t+n到t+2×n的數據。持續擴張直到窮盡所有觀測。初始訓練樣本到最後的觀測之間可以容納多少個n，我們就可以進行多少折交叉驗證。

了解了如何設置交叉驗證，我們將找出Holt-Winters模型的最優參數，回憶一下，我們的廣告數據有按日季節性，所以我們有slen=24。

from sklearn.model_selection import TimeSeriesSplit

def timeseriesCVscore(params, series, loss_function=mean_squared_error, slen=24):

   errors = []

   values = series.values

   alpha, beta, gamma = params

   # 設定交叉驗證折數

   tscv = TimeSeriesSplit(n_splits=3)

   for train, test in tscv.split(values):

       model = HoltWinters(series=values[train], slen=slen,

                           alpha=alpha, beta=beta, gamma=gamma, n_preds=len(test))

       model.triple_exponential_smoothing()

       predictions = model.result[-len(test):]

       actual = values[test]

       error = loss_function(predictions, actual)

       errors.append(error)

return np.mean(np.array(errors))

和其他指數平滑模型一樣，Holt-Winters模型中，平滑參數的取值範圍在0到1之間，因此我們需要選擇一種支持給模型參數添加限制的算法。我們選擇了截斷牛頓共軛梯度（Truncated Newton conjugate gradient）。

data = ads.Ads[:-20] # 留置一些數據用於測試

# 初始化模型參數alpha、beta、gamma

x = [0, 0, 0]

# 最小化損失函數

opt = minimize(timeseriesCVscore, x0=x,

              args=(data, mean_squared_log_error),

              method="TNC", bounds = ((0, 1), (0, 1), (0, 1))

)

# 取最優值……

alpha_final, beta_final, gamma_final = opt.x

print(alpha_final, beta_final, gamma_final)

# ……並據此訓練模型，預測接下來50個小時的數據

model = HoltWinters(data, slen = 24,

                   alpha = alpha_final,

                   beta = beta_final,

                   gamma = gamma_final,

                   n_preds = 50, scaling_factor = 3)

model.triple_exponential_smoothing()

最優參數：

0.11652680227350454 0.002677697431105852 0.05820973606789237

繪圖部分的代碼：

def plotHoltWinters(series, plot_intervals=False, plot_anomalies=False):

   """

       series - 時序數據集

       plot_intervals - 顯示置信區間

       plot_anomalies - 顯示異常值

   """

   plt.figure(figsize=(20, 10))

   plt.plot(model.result, label = "Model")

   plt.plot(series.values, label = "Actual")

   error = mean_absolute_percentage_error(series.values, model.result[:len(series)])

   plt.title("Mean Absolute Percentage Error: {0:.2f}%".format(error))

   if plot_anomalies:

       anomalies = np.array([np.NaN]*len(series))

       anomalies[series.values<model.LowerBond[:len(series)]] = \

           series.values[series.values<model.LowerBond[:len(series)]]

       anomalies[series.values>model.UpperBond[:len(series)]] = \

           series.values[series.values>model.UpperBond[:len(series)]]

       plt.plot(anomalies, "o", markersize=10, label = "Anomalies")

   if plot_intervals:

       plt.plot(model.UpperBond, "r--", alpha=0.5, label = "Up/Low confidence")

       plt.plot(model.LowerBond, "r--", alpha=0.5)

       plt.fill_between(x=range(0,len(model.result)), y1=model.UpperBond,

                        y2=model.LowerBond, alpha=0.2, color = "grey")    

   plt.vlines(len(series), ymin=min(model.LowerBond), ymax=max(model.UpperBond), linestyles='dashed')

   plt.axvspan(len(series)-20, len(model.result), alpha=0.3, color='lightgrey')

   plt.grid(True)

   plt.axis('tight')

   plt.legend(loc="best", fontsize=13);

plotHoltWinters(ads.Ads)

plotHoltWinters(ads.Ads, plot_intervals=True, plot_anomalies=True)

上面的圖形表明，我們的模型能夠很好地逼近初始時序，捕捉每日季節性，總體的下降趨勢，甚至一些異常。如果我們查看下建模偏差（見下圖），我們將很明顯地看到，模型對序列結構的改變反應相當鮮明，但接著很快偏差就回歸正常值，「遺忘」了過去。模型的這一特性讓我們甚至可以為相當噪雜的序列快速構建異常檢測系統，而無需花費過多時間和金錢準備數據和訓練模型。

plt.figure(figsize=(25, 5))

plt.plot(model.PredictedDeviation)

plt.grid(True)

plt.axis('tight')

plt.title("Brutlag's predicted deviation");

遇到異常時偏差會增加

我們將在第二個序列上應用相同的算法，我們知道，第二個序列具有趨勢和每月季節性。

data = currency.GEMS_GEMS_SPENT[:-50]

slen = 30

x = [0, 0, 0]

opt = minimize(timeseriesCVscore, x0=x,

              args=(data, mean_absolute_percentage_error, slen),

              method="TNC", bounds = ((0, 1), (0, 1), (0, 1))

             )

alpha_final, beta_final, gamma_final = opt.x

model = HoltWinters(data, slen = slen,

                   alpha = alpha_final,

                   beta = beta_final,

                   gamma = gamma_final,

                   n_preds = 100, scaling_factor = 3)

model.triple_exponential_smoothing()

看起來很不錯，模型捕捉了向上的趨勢和季節性尖峰，總體而言很好地擬合了數據。

也捕獲了一些異常

偏差隨著預測周期的推進而上升

平穩性

在開始建模之前，我們需要提一下時序的一個重要性質：平穩性（stationarity）。

如果過程是平穩的，那麼它的統計性質不隨時間而變，也就是均值和方差不隨時間改變（方差的恆定性也稱為同方差性），同時協方差函數也不取決於時間（應該只取決於觀測之間的距離）。Sean Abu的博客提供了一些可視化的圖片：

為什麼平穩性如此重要？我們假定未來的統計性質不會和現在觀測到的不同，在平穩序列上做出預測很容易。大多數時序模型多多少少建模和預測這些性質（例如均值和方差），這就是如果原序列不平穩，預測會出錯的原因。不幸的是，我們在教科書以外的地方見到的大多數時序都是不平穩的。不過，我們可以（並且應該）改變這一點。

知己知彼，百戰不殆。為了對抗不平穩性，我們首先需要檢測它。我們現在將查看下白噪聲和隨機遊走，並且了解下如何免費從白噪聲轉到隨機遊走，無需註冊和接受驗證簡訊。

白噪聲圖形：

white_noise = np.random.normal(size=1000)

with plt.style.context('bmh'):  

   plt.figure(figsize=(15, 5))

plt.plot(white_noise)

這一通過標準正態分布生成的過程是平穩的，以0為中心振蕩，偏差為1. 現在我們將基於這一過程生成一個新過程，其中相鄰值之間的關係為：xt = ρxt-1 + et

def plotProcess(n_samples=1000, rho=0):

   x = w = np.random.normal(size=n_samples)

   for t in range(n_samples):

       x[t] = rho * x[t-1] + w[t]

   with plt.style.context('bmh'):  

       plt.figure(figsize=(10, 3))

       plt.plot(x)

       plt.title("Rho {}\n Dickey-Fuller p-value: {}".format(rho, round(sm.tsa.stattools.adfuller(x)[1], 3)))

for rho in [0, 0.6, 0.9, 1]:

plotProcess(rho=rho)

第一張圖上你可以看到之前的平穩的白噪聲。第二張圖的ρ值增加到0.6，導致周期更寬了，但總體上還是平穩的。第三張圖更偏離均值0，但仍以其為中心振蕩。最後，ρ值為1時我們得到了隨機遊走過程——不平穩的時序。

這是因為達到閾值後，時序xt = ρxt-1 + et不再回歸其均值。如果我們從等式的兩邊減去xt-1，我們將得到xt - xt-1 = (ρ-1)xt-1 + et，其中等式左邊的表達式稱為一階差分（first difference）。如果ρ = 1，那麼一階差分將是平穩的白噪聲et。這一事實是時序平穩性的迪基-福勒檢驗（Dickey-Fuller test）的主要思想（檢驗是否存在單位根）。如果非平穩序列可以通過一階差分得到平穩序列，那麼這樣的序列稱為一階單整（integrated of order 1）序列。需要指出的是，一階差分並不總是足以得到平穩序列，因為過程可能是d階單整且d > 1（具有多個單位根），在這樣的情形下，需要使用增廣迪基-福勒檢驗（augmented Dickey-Fuller test）。

我們可以使用不同方法對抗不平穩性——多階差分，移除趨勢和季節性，平滑，以及Box-Cox變換或對數變換。

創建SARIMA擺脫不平穩性

現在，讓我們歷經使序列平穩的多層地獄，創建一個ARIMA模型。

def tsplot(y, lags=None, figsize=(12, 7), style='bmh'):

   """

       繪製時序及其ACF（自相關性函數）、PACF（偏自相關性函數），計算迪基-福勒檢驗

       y - 時序

       lags - ACF、PACF計算所用的時差

   """

   if not isinstance(y, pd.Series):

       y = pd.Series(y)

   with plt.style.context(style):    

       fig = plt.figure(figsize=figsize)

       layout = (2, 2)

       ts_ax = plt.subplot2grid(layout, (0, 0), colspan=2)

       acf_ax = plt.subplot2grid(layout, (1, 0))

       pacf_ax = plt.subplot2grid(layout, (1, 1))

       y.plot(ax=ts_ax)

       p_value = sm.tsa.stattools.adfuller(y)[1]

       ts_ax.set_title('Time Series Analysis Plots\n Dickey-Fuller: p={0:.5f}'.format(p_value))

       smt.graphics.plot_acf(y, lags=lags, ax=acf_ax)

       smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax)

       plt.tight_layout()

tsplot(ads.Ads, lags=60)

上：時序；左下：自相關性；右下：偏自相關性

出乎意料，初始序列是平穩的，迪基-福勒檢驗拒絕了單位根存在的零假設。實際上，從上面的圖形本身就可以看出這一點——沒有可見的趨勢，所以均值是恆定的，整個序列的方差也相對比較穩定。在建模之前我們只需處理季節性。為此讓我們採用「季節差分」，也就是對序列進行簡單的減法操作，時差等於季節周期。

ads_diff = ads.Ads - ads.Ads.shift(24)

tsplot(ads_diff[24:], lags=60)

好多了，可見的季節性消失了，然而自相關函數仍然有過多顯著的時差。為了移除它們，我們將取一階差分：從序列中減去自身（時差為1）

ads_diff = ads_diff - ads_diff.shift(1)

tsplot(ads_diff[24+1:], lags=60)

完美！我們的序列看上去是難以用筆墨形容的完美！在零周圍振蕩，迪基-福勒檢驗表明它是平穩的，ACF中顯著的尖峰不見了。我們終於可以開始建模了！

ARIMA系速成教程

我們將逐字母講解SARIMA(p,d,q)(P,D,Q,s)，季節自回歸移動平均模型（Seasonal Autoregression Moving Average model）：

讓我們把這4個字母組合起來：

AR(p) + MA(q) = ARMA(p,q)

這裡我們得到了自回歸移動平均模型！如果序列是平穩的，我們可以通過這4個字母逼近這一序列。

加上這一字母後我們得到了ARIMA模型，可以通過非季節性差分處理非平穩數據。

加上最後一個字母S後，我們發現這最後一個字母除了s之外，還附帶了三個額外參數——(P,D,Q)。

P —— 模型的季節分量的自回歸階數，同樣可以從PACF得到，但是這次需要查看季節周期長度的倍數的顯著時差的數量。例如，如果周期長度等於24，查看PACF發現第24個時差和第48個時差顯著，那麼初始P值應當是2.

Q —— 移動平均模型的季節分量的階數，初始值的確定和P同理，只不過使用ACF圖形。

D —— 季節性單整階數。等於1或0，分別表示是否應用季節差分。

了解了如何設置初始參數後，讓我們回過頭去重新看下最終的圖形：

p最有可能是4，因為這是PACF上最後一個顯著的時差，之後大多數時差變得不顯著。

d等於1，因為我們採用的是一階差分。

q大概也等於4，這可以從ACF上看出來。

P可能等於2，因為PACF上第24個時差和第48個時差某種程度上比較顯著。

D等於1，我們應用了季節差分。

Q大概是1，ACF上第24個時差是顯著的，而第48個時差不顯著。

現在我們打算測試不同的參數組合，看看哪個是最好的：

# 設定初始值和初始範圍

ps = range(2, 5)

d=1

qs = range(2, 5)

Ps = range(0, 3)

D=1

Qs = range(0, 2)

s = 24

# 創建參數所有可能組合的列表

parameters = product(ps, qs, Ps, Qs)

parameters_list = list(parameters)

len(parameters_list)的結果是54，也就是說，共有54種組合。

def optimizeSARIMA(parameters_list, d, D, s):

   """

   返回參數和相應的AIC的dataframe

       parameters_list - (p, q, P, Q)元組列表

       d - ARIMA模型的單整階

       D - 季節性單整階

       s - 季節長度

   """

      results = []

   best_aic = float("inf")

   for param in tqdm_notebook(parameters_list):

       # 由於有些組合不能收斂，所以需要使用try-except

       try:

           model=sm.tsa.statespace.SARIMAX(ads.Ads, order=(param[0], d, param[1]),

                                           seasonal_order=(param[3], D, param[3], s)).fit(disp=-1)

       except:

           continue

       aic = model.aic

       # 保存最佳模型、AIC、參數

       if aic < best_aic:

           best_model = model

           best_aic = aic

           best_param = param

       results.append([param, model.aic])

   result_table = pd.DataFrame(results)

   result_table.columns = ['parameters', 'aic']

   # 遞增排序，AIC越低越好

   result_table = result_table.sort_values(by='aic', ascending=True).reset_index(drop=True)

   return result_table

result_table = optimizeSARIMA(parameters_list, d, D, s)

# 設定參數為給出最低AIC的參數組合

p, q, P, Q = result_table.parameters[0]

best_model=sm.tsa.statespace.SARIMAX(ads.Ads, order=(p, d, q),

                                       seasonal_order=(P, D, Q, s)).fit(disp=-1)

print(best_model.summary())

讓我們查看下這一模型的殘餘分量（residual）：

tsplot(best_model.resid[24+1:], lags=60)

很明顯，殘餘是平穩的，沒有明顯的自相關性。

讓我們使用這一模型進行預測：

def plotSARIMA(series, model, n_steps):

   """

       繪製模型預測值與實際數據對比圖

       series - 時序數據集

       model - SARIMA模型

       n_steps - 預測未來的步數

   """

   data = series.copy()

   data.columns = ['actual']

   data['arima_model'] = model.fittedvalues

   # 平移s+d步，因為差分的緣故，前面的一些數據沒有被模型觀測到

   data['arima_model'][:s+d] = np.NaN

   forecast = model.predict(start = data.shape[0], end = data.shape[0]+n_steps)

   forecast = data.arima_model.append(forecast)

   # 計算誤差，同樣平移s+d步

   error = mean_absolute_percentage_error(data['actual'][s+d:], data['arima_model'][s+d:])

   plt.figure(figsize=(15, 7))

   plt.title("Mean Absolute Percentage Error: {0:.2f}%".format(error))

   plt.plot(forecast, color='r', label="model")

   plt.axvspan(data.index[-1], forecast.index[-1], alpha=0.5, color='lightgrey')

   plt.plot(data.actual, label="actual")

   plt.legend()

   plt.grid(True);

plotSARIMA(ads, best_model, 50)

最終我們得到了相當不錯的預測，模型的平均誤差率是4.01%，這非常非常好。但是為了達到這一精確度，在準備數據、使序列平穩化、暴力搜索參數上付出了太多。

在我的工作中，創建模型的指導原則常常是快、好、省。這意味著有些模型永遠不會用於生產環境，因為它們需要過長的時間準備數據（比如SARIMA），或者需要頻繁地重新訓練新數據（比如SARIMA），或者很難調整（比如SARIMA）。相反，我經常使用輕鬆得多的方法，從現有時序中選取一些特徵，然後創建一個簡單的線性回歸或隨機森林模型。又快又省。

也許這個方法沒有充分的理論支撐，打破了一些假定（比如，高斯-馬爾可夫定理，特別是誤差不相關的部分），但在實踐中，這很有用，在機器學習競賽中也相當常用。

特徵提取

很好，模型需要特徵，而我們所有的不過是1維時序。我們可以提取什麼特徵？

首先當然是時差。

窗口統計量：

窗口內序列的最大/小值

窗口的平均數/中位數

窗口的方差

等等

日期和時間特徵

目標編碼

其他模型的預測（不過如此預測的話會損失速度）

讓我們運行一些模型，看看我們可以從廣告序列中提取什麼

時序的時差

將序列往回移動n步，我們能得到一個特徵，其中時序的當前值和其t-n時刻的值對齊。如果我們移動1時差，並訓練模型預測未來，那麼模型將能夠提前預測1步。增加時差，比如，增加到6，可以讓模型提前預測6步，不過它需要在觀測到數據的6步之後才能利用。如果在這期間序列發生了根本性的變動，那麼模型無法捕捉這一變動，會返回誤差很大的預測。因此，時差的選取需要平衡預測的質量和時長。

data = pd.DataFrame(ads.Ads.copy())

data.columns = ["y"]

for i in range(6, 25):

data["lag_{}".format(i)] = data.y.shift(i)

from sklearn.linear_model import LinearRegression

from sklearn.model_selection import cross_val_score

# 5折交叉驗證

tscv = TimeSeriesSplit(n_splits=5)

def timeseries_train_test_split(X, y, test_size):

   test_index = int(len(X)*(1-test_size))

   X_train = X.iloc[:test_index]

   y_train = y.iloc[:test_index]

   X_test = X.iloc[test_index:]

   y_test = y.iloc[test_index:]

   return X_train, X_test, y_train, y_test

def plotModelResults(model, X_train=X_train, X_test=X_test, plot_intervals=False, plot_anomalies=False):  

   prediction = model.predict(X_test)

   plt.figure(figsize=(15, 7))

   plt.plot(prediction, "g", label="prediction", linewidth=2.0)

   plt.plot(y_test.values, label="actual", linewidth=2.0)

   if plot_intervals:

       cv = cross_val_score(model, X_train, y_train,

                                   cv=tscv,

                                   scoring="neg_mean_absolute_error")

       mae = cv.mean() * (-1)

       deviation = cv.std()

       scale = 1.96

       lower = prediction - (mae + scale * deviation)

       upper = prediction + (mae + scale * deviation)

       plt.plot(lower, "r--", label="upper bond / lower bond", alpha=0.5)

       plt.plot(upper, "r--", alpha=0.5)

       if plot_anomalies:

           anomalies = np.array([np.NaN]*len(y_test))

           anomalies[y_test<lower] = y_test[y_test<lower]

           anomalies[y_test>upper] = y_test[y_test>upper]

           plt.plot(anomalies, "o", markersize=10, label = "Anomalies")

   error = mean_absolute_percentage_error(prediction, y_test)

   plt.title("Mean absolute percentage error {0:.2f}%".format(error))

   plt.legend(loc="best")

   plt.tight_layout()

   plt.grid(True);

def plotCoefficients(model):

   """

       繪製模型排序後的係數

   """

   coefs = pd.DataFrame(model.coef_, X_train.columns)

   coefs.columns = ["coef"]

   coefs["abs"] = coefs.coef.apply(np.abs)

   coefs = coefs.sort_values(by="abs", ascending=False).drop(["abs"], axis=1)

   plt.figure(figsize=(15, 7))

   coefs.coef.plot(kind='bar')

   plt.grid(True, axis='y')

   plt.hlines(y=0, xmin=0, xmax=len(coefs), linestyles='dashed');

y = data.dropna().y

X = data.dropna().drop(['y'], axis=1)

# 保留30%數據用於測試

X_train, X_test, y_train, y_test = timeseries_train_test_split(X, y, test_size=0.3)

# 機器學習

lr = LinearRegression()

lr.fit(X_train, y_train)

plotModelResults(lr, plot_intervals=True)

plotCoefficients(lr)

模型預測值：綠線為預測值，藍線為實際值

模型係數

好吧，簡單的時差和線性回歸給出的預測質量和SARIMA差得不遠。有大量不必要的特徵，不過我們將在之後進行特徵選擇。現在讓我們繼續增加特徵！

我們將在數據集中加入小時、星期幾、是否周末三個特徵。為此我們需要轉換當前dataframe的索引為datetime格式，並從中提取hour和weekday。

data.index = data.index.to_datetime()

data["hour"] = data.index.hour

data["weekday"] = data.index.weekday

data['is_weekend'] = data.weekday.isin([5,6])*1

可視化所得特徵：

plt.figure(figsize=(16, 5))

plt.title("Encoded features")

data.hour.plot()

data.weekday.plot()

data.is_weekend.plot()

plt.grid(True);

藍線：小時；綠線：星期幾；紅色：是否周末

由於現有的變量尺度不同——時差的長度數千，類別變量的尺度數十——將它們轉換為同一尺度再合理不過，這樣也便於探索特徵重要性，以及之後的正則化。

from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()

y = data.dropna().y

X = data.dropna().drop(['y'], axis=1)

X_train, X_test, y_train, y_test = timeseries_train_test_split(X, y, test_size=0.3)

X_train_scaled = scaler.fit_transform(X_train)

X_test_scaled = scaler.transform(X_test)

lr = LinearRegression()

lr.fit(X_train_scaled, y_train)

plotModelResults(lr, X_train=X_train_scaled, X_test=X_test_scaled, plot_intervals=True)

plotCoefficients(lr)

測試誤差略有下降。從上面的係數圖像我們可以得出結論，weekday和is_weekend是非常有用的特徵。

目標編碼

我想介紹另一種類別變量編碼的變體——基於均值。如果不想讓成噸的獨熱編碼使模型暴漲，同時導致距離信息損失，同時又因為「0點 < 23點」之類的衝突無法使用實數值，那麼我們可以用相對易於解釋的值編碼變量。很自然的一個想法是使用均值編碼目標變量。在我們的例子中，星期幾和一天的第幾小時可以通過那一天或那一小時瀏覽的廣告平均數編碼。非常重要的是，確保均值是在訓練集上計算的（或者交叉驗證當前的折），避免模型偷窺未來。

def code_mean(data, cat_feature, real_feature):

   """

   返回一個字典，鍵為cat_feature的類別，

   值為real_feature的均值。

   """

   return dict(data.groupby(cat_feature)[real_feature].mean())

讓我們看下小時平均：

average_hour = code_mean(data, 'hour', "y")

plt.figure(figsize=(7, 5))

plt.title("Hour averages")

pd.DataFrame.from_dict(average_hour, orient='index')[0].plot()

plt.grid(True);

最後，讓我們定義一個函數完成所有的轉換：

def prepareData(series, lag_start, lag_end, test_size, target_encoding=False):

   data = pd.DataFrame(series.copy())

   data.columns = ["y"]

   for i in range(lag_start, lag_end):

       data["lag_{}".format(i)] = data.y.shift(i)

   data.index = data.index.to_datetime()

   data["hour"] = data.index.hour

   data["weekday"] = data.index.weekday

   data['is_weekend'] = data.weekday.isin([5,6])*1

   if target_encoding:

       test_index = int(len(data.dropna())*(1-test_size))

       data['weekday_average'] = list(map(

           code_mean(data[:test_index], 'weekday', "y").get, data.weekday))

       data["hour_average"] = list(map(

           code_mean(data[:test_index], 'hour', "y").get, data.hour))

       data.drop(["hour", "weekday"], axis=1, inplace=True)

   y = data.dropna().y

   X = data.dropna().drop(['y'], axis=1)

   X_train, X_test, y_train, y_test =\

   timeseries_train_test_split(X, y, test_size=test_size)

   return X_train, X_test, y_train, y_test

X_train, X_test, y_train, y_test =\

prepareData(ads.Ads, lag_start=6, lag_end=25, test_size=0.3, target_encoding=True)

X_train_scaled = scaler.fit_transform(X_train)

X_test_scaled = scaler.transform(X_test)

lr = LinearRegression()

lr.fit(X_train_scaled, y_train)

plotModelResults(lr, X_train=X_train_scaled, X_test=X_test_scaled,

                plot_intervals=True, plot_anomalies=True)

plotCoefficients(lr)

這裡出現了過擬合！Hour_average變量在訓練數據集上表現如此優異，模型決定集中全力在這個變量上——這導致預測質量下降。處理這一問題有多種方法，比如，我們可以不在整個訓練集上計算目標編碼，而是在某個窗口上計算，從最後觀測到的窗口得到的編碼大概能夠更好地描述序列的當前狀態。或者我們可以直接手工移除這一特徵，反正我們已經確定它只會帶來壞處。

X_train, X_test, y_train, y_test =\

prepareData(ads.Ads, lag_start=6, lag_end=25, test_size=0.3, target_encoding=False)

X_train_scaled = scaler.fit_transform(X_train)

X_test_scaled = scaler.transform(X_test)

正則化和特徵選取

正如我們已經知道的那樣，並不是所有的特徵都一樣健康，有些可能導致過擬合。除了手工檢查外我們還可以應用正則化。最流行的兩個帶正則化的回歸模型是嶺（Ridge）回歸和Lasso回歸。它們都在損失函數上施加了某種限制。

在嶺回歸的情形下，限制是係數的平方和，乘以正則化係數。也就是說，特徵係數越大，損失越大，因此優化模型的同時將儘可能地保持係數在較低水平。

因為限制是係數的平方和，所以這一正則化方法稱為L2。它將導致更高的偏差和更低的方差，所以模型的概括性會更好（至少這是我們希望發生的）。

第二種模型Lasso回歸，在損失函數中加上的不是平方和，而是係數絕對值之和。因此在優化過程中，不重要的特徵的係數將變為零，所以Lasso回歸可以實現自動特徵選擇。這類正則化稱為L1。

首先，確認下我們有特徵可以移除，也就是說，確實有高度相關的特徵：

plt.figure(figsize=(10, 8))

sns.heatmap(X_train.corr());

比某些現代藝術要漂亮

from sklearn.linear_model import LassoCV, RidgeCV

ridge = RidgeCV(cv=tscv)

ridge.fit(X_train_scaled, y_train)

plotModelResults(ridge,

                X_train=X_train_scaled,

                X_test=X_test_scaled,

                plot_intervals=True, plot_anomalies=True)

plotCoefficients(ridge)

我們可以很清楚地看到，隨著特徵在模型中的重要性的降低，係數越來越接近零（不過從未達到零）：

lasso = LassoCV(cv=tscv)

lasso.fit(X_train_scaled, y_train)

plotModelResults(lasso,

                X_train=X_train_scaled,

                X_test=X_test_scaled,

                plot_intervals=True, plot_anomalies=True)

plotCoefficients(lasso)

Lasso回歸看起來更保守一點，沒有將第23時差作為最重要特徵，同時完全移除了5項特徵，這提升了預測質量。

XGBoost

為什麼不試試XGBoost？

from xgboost import XGBRegressor

xgb = XGBRegressor()

xgb.fit(X_train_scaled, y_train)

plotModelResults(xgb,

                X_train=X_train_scaled,

                X_test=X_test_scaled,

plot_intervals=True, plot_anomalies=True)

我們的贏家出現了！在我們目前為止嘗試過的模型中，XGBoost在測試集上的誤差是最小的。

不過這一勝利帶有欺騙性，剛到手時序數據，馬上嘗試XGBoost也許不是什麼明智的選擇。一般而言，和線性模型相比，基於樹的模型難以應付數據中的趨勢，所以你首先需要從序列中去除趨勢，或者使用一些特殊技巧。理想情況下，平穩化序列，接著使用XGBoost，例如，你可以使用一個線性模型單獨預測趨勢，然後將其加入XGBoost的預測以得到最終預測。

我們熟悉了不同的時序分析和預測方法。很不幸，或者，很幸運，解決這類問題沒有銀彈。上世紀60年代研發的方法（有些甚至在19世紀就提出了）和LSTM、RNN（本文沒有介紹）一樣流行。這部分是因為時序預測任務和任何其他數據預測任務一樣，在許多方面都需要創造性和研究。儘管有眾多形式化的質量測度和參數估計方法，我們常常需要為每個序列搜尋並嘗試一些不同的東西。最後，平衡質量和成本很重要。之前提到的SARIMA模型是一個很好的例子，經過調節之後，它常常能生成出色的結果，但這需要許多小時、許多複雜技巧來處理數據，相反，10分鐘之內就可以創建好的簡單線性回歸模型卻能取得相當的結果。

杜克大學的高級統計預測課程的在線教材，其中介紹了多種平滑技術、線性模型、ARIMA模型的細節：https://people.duke.edu/~rnau/411home.htm

比較ARIMA和隨機森林預測H5N1高致病性禽流感爆發：https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-15-276

基於Python進行時序分析——從線性模型到圖模型，其中介紹了ARIMA模型家族，及其在建模財經指數上的應用：http://www.blackarbs.com/blog/time-series-analysis-in-python-linear-models-to-garch/11/1/2016

原文地址：https://medium.com/open-machine-learning-course/open-machine-learning-course-topic-9-time-series-analysis-in-python-a270cb05e0b3