時(shí)序分析實(shí)戰(zhàn)之SARIMA、Linear model介紹

1 數(shù)據(jù)準(zhǔn)備

首先引入相關(guān)的 statsmodels，包含統(tǒng)計(jì)模型函數(shù)（時(shí)間序列）。

# 引入相關(guān)的統(tǒng)計(jì)包
import warnings  # 忽略警告
warnings.filterwarnings('ignore')

import numpy as np  # 矢量和矩陣
import pandas as pd  # 表格和數(shù)據(jù)操作
import matplotlib.pyplot as plt
import seaborn as sns
from dateutil.relativedelta import relativedelta  # 有風(fēng)格地處理日期
from scipy.optimize import minimize  # 函數(shù)優(yōu)化
import statsmodels.formula.api as smf  # 統(tǒng)計(jì)與經(jīng)濟(jì)計(jì)量
import statsmodels.tsa.api as smt
import scipy.stats as scs
from itertools import product
from tqdm import tqdm_notebook
import statsmodels.api as sm

用真實(shí)的手機(jī)游戲數(shù)據(jù)作為樣例，研究每小時(shí)觀看的廣告數(shù)和每日所花的游戲幣。

# 1 如真實(shí)的手機(jī)游戲數(shù)據(jù)，將調(diào)查每小時(shí)觀看的廣告和每天花費(fèi)的游戲幣
ads = pd.read_csv(r'./test/ads.csv', index_col=['Time'], parse_dates=['Time'])
currency = pd.read_csv(r'./test/currency.csv', index_col=['Time'], parse_dates=['Time'])

2 穩(wěn)定性

建模前，先來了解一下穩(wěn)定性（stationarity)。

如果一個(gè)過程是平穩(wěn)的，這意味著它不會隨時(shí)間改變其統(tǒng)計(jì)特性，如均值和方差等等。

方差的恒常性稱為同方差，協(xié)方差函數(shù)不依賴于時(shí)間，它只取決于觀測值之間的距離。

非平穩(wěn)過程是指分布參數(shù)或者分布規(guī)律隨時(shí)間發(fā)生變化。也就是說，非平穩(wěn)過程的統(tǒng)計(jì)特征是時(shí)間的函數(shù)（隨時(shí)間變化）。

下面的紅色圖表不是平穩(wěn)的：

• 平均值隨時(shí)間增加

• 方差隨時(shí)間變化

• 隨著時(shí)間的增加，距離變得越來越近。因此，協(xié)方差不是隨時(shí)間而恒定的

為什么平穩(wěn)性如此重要呢? 通過假設(shè)未來的統(tǒng)計(jì)性質(zhì)與目前觀測到的統(tǒng)計(jì)性質(zhì)不會有什么不同，可以很容易對平穩(wěn)序列進(jìn)行預(yù)測。

大多數(shù)的時(shí)間序列模型，以這樣或那樣的方式，試圖預(yù)測那些屬性(例如均值或方差)。

如果原始序列不是平穩(wěn)的，那么未來的預(yù)測將是錯(cuò)誤的。

大多數(shù)時(shí)間序列都是非平穩(wěn)的，但可以(也應(yīng)該)改變這一點(diǎn)。

平穩(wěn)時(shí)間序列的類型：

• 平穩(wěn)過程（stationary process）：產(chǎn)生平穩(wěn)觀測序列的過程。
• 平穩(wěn)模型（stationary model）：描述平穩(wěn)觀測序列的模型。
• 趨勢平穩(wěn)（trend stationary）：不顯示趨勢的時(shí)間序列。
• 季節(jié)性平穩(wěn)（seasonal stationary）：不表現(xiàn)出季節(jié)性的時(shí)間序列。
• 嚴(yán)格平穩(wěn)（strictly stationary）：平穩(wěn)過程的數(shù)學(xué)定義，特別指觀測值的聯(lián)合分布不受時(shí)移的影響。

若時(shí)序中有明顯的趨勢和季節(jié)性，則對這些成分建模，將它們從觀察中剔除，然后用殘差訓(xùn)練建模。

平穩(wěn)性檢查方法（可以檢查觀測值和殘差）：

• 看圖：繪制時(shí)序圖，看是否有明顯的趨勢或季節(jié)性，如繪制頻率圖，看是否呈現(xiàn)高斯分布（鐘形曲線）。
• 概括統(tǒng)計(jì)：看不同季節(jié)的數(shù)據(jù)或隨機(jī)分割或檢查重要的差分，如將數(shù)據(jù)分成兩部分，計(jì)算各部分的均值和方差，然后比較是否一樣或同一范圍內(nèi)
• 統(tǒng)計(jì)測試：選用統(tǒng)計(jì)測試檢查是否有趨勢和季節(jié)性

若時(shí)序的均值和方差相差過大，則有可能是非平穩(wěn)序列，此時(shí)可以對觀測值取對數(shù)，將指數(shù)變化轉(zhuǎn)變?yōu)榫€性變化。然后再次查看取對數(shù)后的觀測值的均值和方差以及頻率圖。

上面前兩種方法常常會欺騙使用者，因此更好的方法是用統(tǒng)計(jì)測試 sm.tsa.stattools.adfuller()。

接下來將學(xué)習(xí)如何檢測穩(wěn)定性，從白噪聲開始。

# 5經(jīng)濟(jì)計(jì)量方法（Econometric approach)
# ARIMA 屬于經(jīng)濟(jì)計(jì)量方法
# 創(chuàng)建平穩(wěn)序列
white_noise = np.random.normal(size=1000)
with plt.style.context('bmh'):
    plt.figure(figsize=(15,5))
    plt.plot(white_noise)

標(biāo)準(zhǔn)正態(tài)分布產(chǎn)生的過程是平穩(wěn)的，在0附近振蕩，偏差為1?，F(xiàn)在，基于這個(gè)過程將生成一個(gè)新的過程，其中每個(gè)后續(xù)值都依賴于前一個(gè)值。

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,5))
        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，大的周期出現(xiàn)，但整體是靜止的。

第三張圖，偏離0均值，但仍然在均值附近震蕩。

第四張圖，ρ= 1，有一個(gè)隨機(jī)游走過程即非平穩(wěn)時(shí)間序列。當(dāng)達(dá)到臨界值時(shí)， 不返回其均值。從兩邊減去 ，得到 ，左邊的表達(dá)式被稱為 一階差分 。

如果 ρ= 1，那么一階差分等于平穩(wěn)白噪聲 。這是 Dickey-Fuller時(shí)間序列平穩(wěn)性測試 (測試單位根的存在)背后的主要思想。

如果 可以用一階差分從非平穩(wěn)序列中得到平穩(wěn)序列，稱這些序列為1階積分 。 該檢驗(yàn)的零假設(shè)是時(shí)間序列是非平穩(wěn)的 ，在前三個(gè)圖中被拒絕，在最后一個(gè)圖中被接受。

1階差分并不總是得到一個(gè)平穩(wěn)的序列，因?yàn)檫@個(gè)過程可能是 d 階的積分，d > 1階的積分(有多個(gè)單位根)。在這種情況下，使用增廣的Dickey-Fuller檢驗(yàn)，它一次檢查多個(gè)滯后時(shí)間。

可以使用不同的方法來去除非平穩(wěn)性：各種順序差分、趨勢和季節(jié)性去除、平滑以及轉(zhuǎn)換如Box-Cox或?qū)?shù)轉(zhuǎn)換。

3 SARIMA

接下來開始建立ARIMA模型，在建模之前需要將非平穩(wěn)時(shí)序轉(zhuǎn)換為平穩(wěn)時(shí)序。

3.1 去除非穩(wěn)定性

擺脫非平穩(wěn)性，建立SARIMA（Getting rid of non-stationarity and building SARIMA）

def tsplot(y,lags=None,figsize=(12,7),style='bmh'):
"""
Plot time series, its ACF and PACF, calculate Dickey-Fuller test
y:timeseries
lags:how many lags to include in ACF,PACF calculation
"""
ifnot 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)

從圖中可以看出，Dickey-Fuller檢驗(yàn)拒絕了單位根存在的原假設(shè)（p=0）；序列是平穩(wěn)的，沒有明顯的趨勢，所以均值是常數(shù)，方差很穩(wěn)定。

唯一剩下的是季節(jié)性，必須在建模之前處理它?？梢酝ㄟ^季節(jié)差分去除季節(jié)性，即序列本身減去一個(gè)滯后等于季節(jié)周期的序列。

ads_diff = ads.Ads-ads.Ads.shift(24) # 去除季節(jié)性
tsplot(ads_diff[24:], lags=60)

ads_diff = ads_diff - ads_diff.shift(1) # 去除趨勢
tsplot(ads_diff[24+1:], lags=60) # 最終圖

第一張圖中，隨著季節(jié)性的消失，自回歸好多了，但是仍存在太多顯著的滯后，需要?jiǎng)h除。首先使用一階差分，用滯后1從自身中減去時(shí)序。

第二張圖中，通過季節(jié)差分和一階差分得到的序列在0周圍震蕩。Dickey-Fuller試驗(yàn)表明，ACF是平穩(wěn)的，顯著峰值的數(shù)量已經(jīng)下降，可以開始建模了。

3.2 建 SARIMA 模型

SARIMA：Seasonal Autoregression Integrated Moving Average model。

是簡單自回歸移動(dòng)平均的推廣，并增加了積分的概念。

• AR (p): 利用一個(gè)觀測值和一些滯后觀測值之間的依賴關(guān)系的模型。模型中的最大滯后稱為p。要確定初始p，需要查看PACF圖并找到最大的顯著滯后，在此之后大多數(shù)其他滯后都變得不顯著。
• I (d): 利用原始觀測值的差值(如觀測值減去上一個(gè)時(shí)間步長的觀測值)使時(shí)間序列保持平穩(wěn)。這只是使該系列固定所需的非季節(jié)性差分的數(shù)量。在例子中它是1，因?yàn)槭褂昧艘浑A差分。
• MA (q):利用觀測值與滯后觀測值的移動(dòng)平均模型殘差之間的相關(guān)性的模型。目前的誤差取決于前一個(gè)或前幾個(gè)，這被稱為q。初始值可以在ACF圖上找到，其邏輯與前面相同

每個(gè)成分都對應(yīng)著相應(yīng)的參數(shù)。

SARIMA(p,d,q)(P,D,Q,s) 模型需要選擇趨勢和季節(jié)的超參數(shù)。

趨勢參數(shù) ，趨勢有三個(gè)參數(shù)，與ARIMA模型的參數(shù)一樣：

• p: 模型中包含的滯后觀測數(shù)，也稱滯后階數(shù)。
• d: 原始觀測值被差值的次數(shù)，也稱為差值度。
• q:移動(dòng)窗口的大小，也叫移動(dòng)平均的階數(shù)

季節(jié)參數(shù) ：

• S(s):負(fù)責(zé)季節(jié)性，等于單個(gè)季節(jié)周期的時(shí)間步長
• P:模型季節(jié)分量的自回歸階數(shù)，可由PACF推導(dǎo)得到。但是需要看一下顯著滯后的次數(shù)，是季節(jié)周期長度的倍數(shù)。如果周期等于24，看到24和48的滯后在PACF中是顯著的，這意味著初始P應(yīng)該是2。P=1將利用模型中第一次季節(jié)偏移觀測，如t-(s*1)或t-24；P=2，將使用最后兩個(gè)季節(jié)偏移觀測值t-(s * 1)， t-(s * 2)
• Q:使用ACF圖實(shí)現(xiàn)類似的邏輯
• D:季節(jié)性積分的階數(shù)（次數(shù)）。這可能等于1或0，取決于是否應(yīng)用了季節(jié)差分。

可以通過分析ACF和PACF圖來選擇趨勢參數(shù)，查看最近時(shí)間步長的相關(guān)性（如1，2，3）。

同樣，可以分析ACF和PACF圖，查看季節(jié)滯后時(shí)間步長的相關(guān)性來指定季節(jié)模型參數(shù)的值。

現(xiàn)在知道了如何設(shè)置初始參數(shù)，查看最終圖并設(shè)置參數(shù)。

上面倒數(shù)第一張圖就是最終圖：

• p:最可能是4，因?yàn)樗荘ACF上最后一個(gè)顯著的滯后，在這之后，大多數(shù)其他的都不顯著。
• d:為1，因?yàn)槲覀冇?jì)算了一階差分
• q:應(yīng)該在4左右，就像ACF上看到的那樣
• P:可能是2，因?yàn)?4和48的滯后對PACF有一定的影響
• D:為1，因?yàn)橛?jì)算了季節(jié)差分
• Q:可能1，ACF的第24個(gè)滯后顯著，第48個(gè)滯后不顯著。

下面看不同參數(shù)的模型表現(xiàn)如何：

# 建模 SARIMA
# setting initial values and some bounds for them
ps = range(2,5)
d=1
qs=range(2,5)
Ps=range(0,2)
D=1
Qs=range(0,2)
s=24#season length

# creating list with all the possible combinations of parameters
parameters=product(ps,qs,Ps,Qs)
parameters_list = list(parameters)
print(parameters)
print(parameters_list)
print(len(parameters_list))
# 36

def optimizeSARIMA(parameters_list, d,D,s):
    """
    Return dataframe with parameters and corresponding AIC
    parameters_list:list with (p,q,P,Q) tuples
    d:integration order in ARIMA model
    D:seasonal integration order
    s:length of season
    """
    results = []
    best_aic = float('inf')

    for param in tqdm_notebook(parameters_list):
        # we need try-exccept because on some combinations model fails to converge
        try:
            model = sm.tsa.statespace.SARIMAX(ads.Ads, order=(param[0], d,param[1]),
                                              seasonal_order=(param[2], D,param[3], s)).fit(disp=-1)
        except:
            continue
        aic = model.aic
        # saving best model, AIC and parameters
        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']
    # sorting in ascending order, the lower AIC is - the better
    result_table = result_table.sort_values(by='aic',ascending=True).reset_index(drop=True)

    return result_table

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

result_table.head()

# set the parameters that give the lowerst 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()) # 打印擬合模型的摘要。這總結(jié)了所使用的系數(shù)值以及對樣本內(nèi)觀測值的擬合技巧。

# inspect the residuals of the model
tsplot(best_model.resid[24+1:], lags=60)

parameters          aic
0  (2, 3, 1, 1)  3888.642174
1  (3, 2, 1, 1)  3888.763568
2  (4, 2, 1, 1)  3890.279740
3  (3, 3, 1, 1)  3890.513196
4  (2, 4, 1, 1)  3892.302849

                                Statespace Model Results
==========================================================================================
Dep. Variable:                                Ads   No. Observations:                  216
Model:             SARIMAX(2, 1, 3)x(1, 1, 1, 24)   Log Likelihood               -1936.321
Date:                            Sun, 08 Mar 2020   AIC                           3888.642
Time:                                    23:06:23   BIC                           3914.660
Sample:                                09-13-2017   HQIC                          3899.181
                                     - 09-21-2017
Covariance Type:                              opg
==============================================================================
                 coef    std err          z      P >|z|      [0.0250.975]
------------------------------------------------------------------------------
ar.L1          0.79130.2702.9280.0030.2621.321
ar.L2         -0.55030.306-1.7990.072-1.1500.049
ma.L1         -0.73160.262-2.7930.005-1.245-0.218
ma.L2          0.56510.2822.0050.0450.0131.118
ma.L3         -0.18110.092-1.9640.049-0.362-0.000
ar.S.L24       0.33120.0764.3510.0000.1820.480
ma.S.L24      -0.76350.104-7.3610.000-0.967-0.560
sigma2      4.574e+075.61e-098.15e+150.0004.57e+074.57e+07
===================================================================================
Ljung-Box (Q):                       43.70   Jarque-Bera (JB):                10.56
Prob(Q):                              0.32   Prob(JB):                         0.01
Heteroskedasticity (H):               0.65   Skew:                            -0.28
Prob(H) (two-sided):                  0.09   Kurtosis:                         4.00
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
[2] Covariance matrix is singular or near-singular, with condition number 8.82e+31. Standard errors may be unstable.

很明顯，殘差是平穩(wěn)的，不存在明顯的自相關(guān)?？捎梦覀兊哪Ｐ蛠眍A(yù)測。

def plotSARIMA(series, model, n_steps):
    """
    plot model vs predicted values
    series:dataset with timeseries
    model:fitted SARIMA model
    n_steps:number of steps to predict in the future
    """
    # adding model values
    data = series.copy()
    data.columns = ['actual']
    data['arima_model']=model.fittedvalues
    # making a shift on s+d steps, because these values were unobserved by the model due the differentiating
    data['arima_model'][:s+d]=np.nan

    # forecasting on n_steps forward
    forecast = model.predict(start=data.shape[0],end=data.shape[0]+n_steps)
    forecast = data.arima_model.append(forecast)
    # calculate error, again having shifted on s+d steps from the beginning
    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)

最后，得到了非常充分的預(yù)測。模型平均誤差是3.94%，非常好。但準(zhǔn)備數(shù)據(jù)、使序列平穩(wěn)和選擇參數(shù)的總成本可能不值得這么精確。

這一過程的步驟總結(jié)如下：

• 模式識別：使用圖表和匯總統(tǒng)計(jì)數(shù)據(jù)來確定趨勢、季節(jié)性和自回歸元素，從而了解差分的大小和所需的滯后的大小。
• 參數(shù)估計(jì)：利用擬合程序求出回歸模型的系數(shù)。
• 模型檢查：利用剩余誤差的圖和統(tǒng)計(jì)檢驗(yàn)來確定模型沒有捕捉到的時(shí)間結(jié)構(gòu)的數(shù)量和類型。

4 線性模型

通常，在工作中必須以快速、好作為指導(dǎo)原則來建立模型。

一些模型不能拿來就用，因?yàn)樗鼈冃枰嗟臄?shù)據(jù)準(zhǔn)備時(shí)間(如SARIMA)，或者需要對新數(shù)據(jù)進(jìn)行頻繁的訓(xùn)練(SARIMA)，或者很難調(diào)參(SARIMA)。因此，從現(xiàn)有的時(shí)間序列中選擇幾個(gè)特征并構(gòu)建一個(gè)簡單的線性回歸模型，或者一個(gè)隨機(jī)森林，通常要容易得多。

這種方法沒有理論支持，并且打破了幾個(gè)假設(shè)(例如高斯-馬爾科夫定理，尤其是誤差不相關(guān)的情況下)，但是在實(shí)踐中非常有用，經(jīng)常在機(jī)器學(xué)習(xí)競賽中使用。

4.1 特征提取

模型需要特征，但只有一維的時(shí)間序列。可以提取那些特征？

• 時(shí)間序列窗口統(tǒng)計(jì)的滯后 ：一個(gè)窗口中序列的最大值/最小值、一個(gè)窗口中值、窗口方差等。
• 日期和時(shí)間特征 ：一小時(shí)的第幾分鐘，一天的第幾小時(shí)，一周的第幾天等等
• 特別的活動(dòng) ：將其表示為布爾特征(注意，這種方式可能失去預(yù)測的速度)。讓我們運(yùn)行一些方法，看看我們可以從ads時(shí)間序列數(shù)據(jù)中提取什么。

接下來用這些方法提取特征。

4.1.1 滯后特征

將序列向后移動(dòng) n個(gè)時(shí)間點(diǎn)，得到一個(gè)特征列，其中時(shí)間序列的當(dāng)前值與其在時(shí)間 t-n 時(shí)的值一致。

若進(jìn)行一個(gè) 1 延遲移位，并對該特征訓(xùn)練模型，則該模型能夠在觀察到該序列的當(dāng)前狀態(tài)后提前預(yù)測1步。

增加滯后，比如增加到6，將允許模型提前6步做出預(yù)測，將使用6步前觀察到的數(shù)據(jù)。

如果在這段未觀察到的時(shí)間內(nèi)，有什么東西從根本上改變了這個(gè)序列，那么模型就不會捕捉到這些變化，并且會返回帶有很大誤差的預(yù)測。因此，在初始滯后選擇過程中，必須在最優(yōu)預(yù)測質(zhì)量與預(yù)測長度之間找到一個(gè)平衡點(diǎn)。

# 其他模型：線性...
# Creating a copy of the initial dataframe to make various transformations
data = pd.DataFrame(ads.Ads.copy())
data.columns=['y']

# Adding the lag of the target variable from 6 setps back up to 24
for i in range(6,25):
    data['lag_{}'.format(i)]=data.y.shift(i)

# take a look at the new dataframe
data.tail(7)

y     lag_6    ...       lag_23    lag_24
Time                                     ...
2017-09-2117:00:00151790132335.0    ...     146215.0139515.0
2017-09-2118:00:00155665146630.0    ...     142425.0146215.0
2017-09-2119:00:00155890141995.0    ...     123945.0142425.0
2017-09-2120:00:00123395142815.0    ...     101360.0123945.0
2017-09-2121:00:00103080146020.0    ...      88170.0101360.0
2017-09-2122:00:0095155152120.0    ...      76050.088170.0
2017-09-2123:00:0080285151790.0    ...      70335.076050.0
[7 rows x 20 columns]

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score

# for time-series cross-validation set 5 folds
tscv = TimeSeriesSplit(n_splits=5)

def timeseries_train_test_split(X,y,test_size):
    """
    Perform train-test split with respect to time series structure
    """
    # get the index after which test set starts
    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

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

# reserve 30% of data for testing
X_train, X_test, y_train, y_test =timeseries_train_test_split(X,y,test_size=0.3)

#-------------------------------------------------------------------------------------
# machine learning in two lines
lr = LinearRegression()
lr.fit(X_train, y_train)

#-------------------------------------------------------------------------------------
def plotModelResults(model, X_train=X_train, X_test=X_test, plot_intervals=False, plot_anomalies=False):
    """
    Plot modelled vs fact values, prediction intervals and anomalies
    """
    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):
    """
    PLots sorted coefficient values of the model
    """

    coefs = pd.DataFrame(model.coef_, X_train.columns)
    print()
    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')

#-------------------------------------------------------------------------------------
plotModelResults(lr, plot_intervals=True)
plotCoefficients(lr)

簡單的滯后和線性回歸的預(yù)測在質(zhì)量方面與SARIMA差不多。有許多不必要的特征，稍后將進(jìn)行特征選擇。

接下來提取時(shí)間特征。

4.1.2 時(shí)間特征

# 提取時(shí)間特征 hour、day of week、is_weekend
data.index = pd.to_datetime(data.index)
data['hour'] = data.index.hour
data['weekday'] = data.index.weekday
data['is_weekend'] = data.weekday.isin([5,6])*1
data.tail()

#-------------------------------------------------------------------------------------
# 可視化特征
plt.figure(figsize=(16,5))
plt.title('Encoded features')
data.hour.plot()
data.weekday.plot()
data.is_weekend.plot()
plt.grid(True)

y     lag_6     ...      weekday  is_weekend
Time                                      ...
2017-09-2119:00:00155890141995.0     ...            30
2017-09-2120:00:00123395142815.0     ...            30
2017-09-2121:00:00103080146020.0     ...            30
2017-09-2122:00:0095155152120.0     ...            30
2017-09-2123:00:0080285151790.0     ...            30
[5 rows x 23 columns]

因?yàn)楝F(xiàn)在的變量中有不同的尺度，滯后特征有數(shù)千，分類特征有數(shù)十，我們需要將它們轉(zhuǎn)換成相同的尺度，以探索特征的重要性，然后是正則化。

4.2 歸一化

# 特征的尺度不一樣，需要?dú)w一化

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)

4.3 目標(biāo)編碼

def code_mean(data, cat_feature, real_feature):
    """
    Returns a dictionary where keys are unique categories of the cat_feature,
    and values are means over 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)

把所有的變換放在一個(gè)函數(shù)中：

# 將所有的數(shù)據(jù)準(zhǔn)備結(jié)合到一起
def prepareData(series, lag_start, lag_end, test_size, target_encoding=False):
    """
    series: pd.DataFrame or dataframe with timeseries
    lag_start: int, initial step back in time to slice target variable
               example-lag_start=1 means that the model will see yesterday's values to predict today
    lag_end: int, finial step back in time to slice target variable
             example-lag_end=4 means that the model will see up to 4 days back in time to predict today
    test_size:float, size of the test dataset after train/test split as percentage of dataset
    target_encoding:boolean, if True -  add target averages to dataset
    """

    # copy of the initial dataset
    data = pd.DataFrame(series.copy())
    data.columns = ['y']

    # lags of series
    for i in range(lag_start, lag_end):
        data['lag_{}'.format(i)]=data.y.shift(i)

    # datatime features
    data.index = pd.to_datetime(data.index)
    data['hour'] = data.index.hour
    data['weekday'] =data.index.weekday
    data['is_weekend']=data.weekday.isin([5,6])*1

    if target_encoding:
        # calculate averages on train set only
        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))

        # drop encoded variables
        data.drop(['hour', 'weekday'], axis=1, inplace=True)

    # train-test split
    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)
# plt.xticks(rotation=45, fontsize=7)

從圖中可以看出，存在過擬合。Hour_average在訓(xùn)練數(shù)據(jù)集中系數(shù)過大，以至于模型決定都集中在它上面。結(jié)果，預(yù)測的質(zhì)量下降了。這個(gè)問題可以用多種方法解決：例如，我們可以計(jì)算目標(biāo)編碼而不是整個(gè)訓(xùn)練集，而是針對某個(gè)窗口。這樣，來自最后一個(gè)觀察到的窗口的編碼將很可能更好地描述當(dāng)前的序列狀態(tài)。或者，可以手動(dòng)刪除它，因?yàn)槲覀兇_信在這種情況下它只會使事情變得更糟。

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)

4.4 正則化

并不是所有的特征都是有用的，有些可能會導(dǎo)致過擬合，有些應(yīng)該被刪除。

除了人工檢查，還可以采用正則化。

兩個(gè)有名的正則化回歸模型是嶺回歸和套索回歸 ^[2]^ 。它們都給損失函數(shù)增加了一些約束。

在嶺回歸的情況下，這些約束是系數(shù)乘正則化系數(shù)的平方和。一個(gè)特征的系數(shù)越大，損失就越大。因此，盡量優(yōu)化模型，同時(shí)保持系數(shù)相當(dāng)?shù)?。L2正則化的結(jié)果，將有更高的偏差和更低的方差，因此模型泛化性能更好。

第二個(gè)回歸模型是Lasso回歸，它將損失函數(shù)添加到系數(shù)的絕對值中，而不是平方。因此，在優(yōu)化過程中，不重要的特征的系數(shù)可能變成零，從而允許自動(dòng)選擇特征。這種正則化類型稱為L1。

首先，確定有要?jiǎng)h除的特征，并且數(shù)據(jù)具有高度相關(guān)的特征。

# 若過擬合，對特征進(jìn)行正則化
# 先看一下特征之間的相關(guān)性
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)

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)

可以清楚地看到，隨著它們在模型中的重要性下降，一些系數(shù)正越來越接近于零(盡管它們從未真正達(dá)到零)。

# 套索回歸
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)

套索回歸被證明是更保守的。它從最重要的特征中去除了23的延遲，并完全去掉了5個(gè)特征，這提高了預(yù)測的質(zhì)量。

下面嘗試用XGBoost建模。

5 Boosting

# 預(yù)測模型 Boosting
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)

這是到目前為止測試過的所有模型中測試集上誤差最小的。但是，這具有欺騙性。得到時(shí)間序列數(shù)據(jù)，不適合先用xgboost。

通常，與線性模型相比，基于樹的模型處理數(shù)據(jù)趨勢的能力較差。在這種情況下，您必須先停止使用序列，或者使用一些技巧來實(shí)現(xiàn)這個(gè)魔術(shù)。

理想情況下，您可以使該系列平穩(wěn)，然后使用XGBoost。例如，可以使用線性模型預(yù)測趨勢，然后加上xgboost的預(yù)測以獲得最終預(yù)測。