1 數據準備

首先引入相關的 statsmodels，包含統計模型函數（時間序列）。

# 引入相關的統計包
import warnings  # 忽略警告
warnings.filterwarnings('ignore')

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 scipy.stats as scs
from itertools import product
from tqdm import tqdm_notebook
import statsmodels.api as sm

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

# 1 如真實的手機游戲數據，將調查每小時觀看的廣告和每天花費的游戲幣
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 穩定性

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

如果一個過程是平穩的，這意味著它不會隨時間改變其統計特性，如均值和方差等等。

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

非平穩過程是指分布參數或者分布規律隨時間發生變化。也就是說，非平穩過程的統計特征是時間的函數（隨時間變化）。

下面的紅色圖表不是平穩的：

• 平均值隨時間增加

• 方差隨時間變化

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

為什么平穩性如此重要呢? 通過假設未來的統計性質與目前觀測到的統計性質不會有什么不同，可以很容易對平穩序列進行預測。

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

如果原始序列不是平穩的，那么未來的預測將是錯誤的。

大多數時間序列都是非平穩的，但可以(也應該)改變這一點。

平穩時間序列的類型：

• 平穩過程（stationary process）：產生平穩觀測序列的過程。
• 平穩模型（stationary model）：描述平穩觀測序列的模型。
• 趨勢平穩（trend stationary）：不顯示趨勢的時間序列。
• 季節性平穩（seasonal stationary）：不表現出季節性的時間序列。
• 嚴格平穩（strictly stationary）：平穩過程的數學定義，特別指觀測值的聯合分布不受時移的影響。

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

平穩性檢查方法（可以檢查觀測值和殘差）：

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

若時序的均值和方差相差過大，則有可能是非平穩序列，此時可以對觀測值取對數，將指數變化轉變為線性變化。然后再次查看取對數后的觀測值的均值和方差以及頻率圖。

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

接下來將學習如何檢測穩定性，從白噪聲開始。

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

標準正態分布產生的過程是平穩的，在0附近振蕩，偏差為1。現在，基于這個過程將生成一個新的過程，其中每個后續值都依賴于前一個值。

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，大的周期出現，但整體是靜止的。

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

第四張圖，ρ= 1，有一個隨機游走過程即非平穩時間序列。當達到臨界值時， 不返回其均值。從兩邊減去 ，得到 ，左邊的表達式被稱為 一階差分 。

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

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

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

可以使用不同的方法來去除非平穩性：各種順序差分、趨勢和季節性去除、平滑以及轉換如Box-Cox或對數轉換。

3 SARIMA

接下來開始建立ARIMA模型，在建模之前需要將非平穩時序轉換為平穩時序。

3.1 去除非穩定性

擺脫非平穩性，建立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檢驗拒絕了單位根存在的原假設（p=0）；序列是平穩的，沒有明顯的趨勢，所以均值是常數，方差很穩定。

唯一剩下的是季節性，必須在建模之前處理它。可以通過季節差分去除季節性，即序列本身減去一個滯后等于季節周期的序列。

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

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

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

第二張圖中，通過季節差分和一階差分得到的序列在0周圍震蕩。Dickey-Fuller試驗表明，ACF是平穩的，顯著峰值的數量已經下降，可以開始建模了。

3.2 建 SARIMA 模型

SARIMA：Seasonal Autoregression Integrated Moving Average model。

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

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

每個成分都對應著相應的參數。

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

趨勢參數 ，趨勢有三個參數，與ARIMA模型的參數一樣：

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

季節參數 ：

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

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

同樣，可以分析ACF和PACF圖，查看季節滯后時間步長的相關性來指定季節模型參數的值。

現在知道了如何設置初始參數，查看最終圖并設置參數。

上面倒數第一張圖就是最終圖：

• p:最可能是4，因為它是PACF上最后一個顯著的滯后，在這之后，大多數其他的都不顯著。
• d:為1，因為我們計算了一階差分
• q:應該在4左右，就像ACF上看到的那樣
• P:可能是2，因為24和48的滯后對PACF有一定的影響
• D:為1，因為計算了季節差分
• Q:可能1，ACF的第24個滯后顯著，第48個滯后不顯著。

下面看不同參數的模型表現如何：

# 建模 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()) # 打印擬合模型的摘要。這總結了所使用的系數值以及對樣本內觀測值的擬合技巧。

# 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.

很明顯，殘差是平穩的，不存在明顯的自相關。可用我們的模型來預測。

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)

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

這一過程的步驟總結如下：

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

4 線性模型

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

一些模型不能拿來就用，因為它們需要太多的數據準備時間(如SARIMA)，或者需要對新數據進行頻繁的訓練(SARIMA)，或者很難調參(SARIMA)。因此，從現有的時間序列中選擇幾個特征并構建一個簡單的線性回歸模型，或者一個隨機森林，通常要容易得多。

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

4.1 特征提取

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

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

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

4.1.1 滯后特征

將序列向后移動 n個時間點，得到一個特征列，其中時間序列的當前值與其在時間 t-n 時的值一致。

若進行一個 1 延遲移位，并對該特征訓練模型，則該模型能夠在觀察到該序列的當前狀態后提前預測1步。

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

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

# 其他模型：線性...
# 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)

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

接下來提取時間特征。

4.1.2 時間特征

# 提取時間特征 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]

因為現在的變量中有不同的尺度，滯后特征有數千，分類特征有數十，我們需要將它們轉換成相同的尺度，以探索特征的重要性，然后是正則化。

4.2 歸一化

# 特征的尺度不一樣，需要歸一化

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 目標編碼

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)

把所有的變換放在一個函數中：

# 將所有的數據準備結合到一起
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_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 正則化

并不是所有的特征都是有用的，有些可能會導致過擬合，有些應該被刪除。

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

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

在嶺回歸的情況下，這些約束是系數乘正則化系數的平方和。一個特征的系數越大，損失就越大。因此，盡量優化模型，同時保持系數相當低。L2正則化的結果，將有更高的偏差和更低的方差，因此模型泛化性能更好。

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

首先，確定有要刪除的特征，并且數據具有高度相關的特征。

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

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

# 套索回歸
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個特征，這提高了預測的質量。

下面嘗試用XGBoost建模。

5 Boosting

# 預測模型 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)

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

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

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