時間序列（二）：時序預測那些事兒

在上一篇文章中，我們簡略介紹了與時間序列相關的應用，這次我們聚焦於時間序列的預測，講講與之相關的那些事。

1. 效果評估

設 y 是時間序列的真實值， yhat 是模型的預測值。在分類模型中由於y是離散的，有很多維度可以去刻畫預測的效果。但現在的y是連續的，工具一下子就少了很多。時間序列里比較常用的是 MAPE（mean absolute percentage error） 和 RMSE （root mean square error）

\[MAPE = \frac{100\%}{n} \sum_{i=1}^{n} \big|\frac{y-yhat}{y}\big| \]

\[RMSE = \sqrt{\frac{\sum(y-yhat)^2}{n}} \]

個人比較喜歡MAPE，適合匯報，它揭示了預測值在真實值基礎上的波動范圍。比如，MAPE=10%，其代表着預測值的范圍大致是 0.9y 至 1.1y

這兩個指標或多或少都會有一些缺點，因此也有了各種變體,比如SMAPE, NRMSE等，感興趣的可以自己去Wiki看看。

\[SMAPE = \frac{200\%}{n} \sum_{i=1}^{n} \frac{|y-yhat|}{|y|+|yhat|} \]

2. 交叉驗證

因為時間序列的樣本之間是無法交換的，所以沒辦法隨機的像 KFold 一樣把數據集切分成若干份訓練集和測試集。但沒關系，方法還是有的。

一個比較好的思路是按照時間順序設置 cutoff time T_c :

這里有三個參數：

• horizon: 模型預測的范圍，如從cutoff點開始數未來30天

• period: 每兩個 cutoff 點之間的間隔，如60天

• initial: 用於訓練的日期范圍，如730天

假設我們有2015-01-01 至 2017-12-31 共三年的數據。且三個參數的值如下： horizon=30,period=120,initial=730,則我們有

\[K=\text{floor} \big(\frac{365*3 - initial -period}{period}\big) = 5 \]

5組訓練集和測試集。

有了交叉驗證，我們就能獲得評估模型效果和進行超參數優化了.

3. 時間序列預測模型

時間序列的預測有很多方法，本文不准備羅列的很細致，只展示三種模型：ARIMA,PROPHET,ETS,且都只用默認參數。在之后的一篇文章中，我會詳細講講prophet的應用。

時間序列的預測模型，不好說一個模型絕對的好，跟數據集的類型有很大關系。根據更新的頻率，可以分為Tick數據（小於1s,金融中的高頻交易）、分鍾數據、日數據、月數據等；根據波動情況，可以分為平穩序列和對日期/事件敏感的季節性序列等。一般數據的顆粒度越小，預測效果越差，因為隨機性太大了，想通過模型來預測

我找了一個對季節和節假日敏感的日數據來測試這幾個模型，時間范圍是5.3年。可以看到數據有一個整體的趨勢，也有季節性趨勢，還有那些尖峰對應的節假日/事件影響。

3.1 ARIMA 模型

在講ARIMA模型之前得先熟悉ARMA模型(autoregressive moving average model)

ARMA 模型是經典的預測模型，有着成熟的理論基礎，但條件也比較嚴格，就是要求時間序列是平穩的。這里講的平穩性條件(一般指弱平穩)是指時間序列的均值與時間無關和方差僅與時間差有關。

在ARMA模型中，時序可以由如下方程表示：

\[Y_t = c+\varepsilon_{t}+\sum_{i=1}^{p}\alpha_iY_{t-i} + \sum_{i=1}^{q}\theta_{i}\varepsilon_{t-i} \]

其中 \epsilon 是高斯白噪聲序列。

可以從線性系統的角度去看ARMA 模型，已證明其是一個離散線性時不變系統，系統的輸入是高斯白噪聲序列(標准的平穩序列)，輸出是我們需要的時序y。此時上述式子可以寫為

\[(1-\sum_{i=1}^{p} \alpha_i L^i)Y_t =(1+\sum_{i=1}^{q} \theta_i L^i )\varepsilon_t \]

其中 L 是單位平移算子(lag operator).

有定理指出，對於一個離散線性平移不變系統而言，當輸入是平穩序列時，可以證明輸出也是平穩的，且輸入序列和輸出序列是平穩相關的。（這里面的東西還比較復雜，本科學的，現在已經忘了差不多了，感興趣的同學可以自己去看看線性系統的理論）

商業中的數據大都不滿足平穩性條件的，為了解決這個問題，所以就衍生了ARIMA 模型，中間的字母I代表的是差分。通過一些差分/季節性差分操作，我們能消除掉大部分趨勢，使得原有不平穩的序列變得平穩（當然也不絕對，一般需要在差分后作平穩性檢驗）。

這里我設置交叉驗證的參數如下：horizon,initial,period=30,1500,30 , 並生成12組訓練集和測試集用來評估模型效果。

import pandas as pd
import numpy as np
from pyecharts import Bar
from pyecharts import Line
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
%matplotlib inline
from pyramid.arima import auto_arima
from fbprophet import Prophet
import ts_model_selection ##自己寫的輪子，還沒整理好，之后會放出來

df=pd.read_excel('.\\data\\tehr.xlsx')
df['ds']=pd.to_datetime(df['ds'])
df['logy']=np.log(df['y'])

horizon,initial,period=30,1500,30
df_train_list,df_test_list=ts_model_selection.train_test_split(\
    df,horizon='{} days'.format(horizon),initial='{} days'.format(initial),period='{} days'.format(period))

# 記錄各個模型的分數
result={} #method,mape,cv

arima 有很多需要調參的地方，如p(AR),d(差分),q(MA)等。在 statsmodels 包中就有很多 ARIMA 相關的函數，但遺憾都需要手動調參。無奈找了一個仿照R中auto.arima寫的包，還不錯，名字是 pyramid-arima。接下來我們就用它來預測。

forecasts=pd.DataFrame(columns=['y','yhat','k'])
for k in range(len(df_train_list)):
    #print(k)
    df_train,df_test=df_train_list[k],df_test_list[k]
    stepwise_fit = auto_arima(df_train['y'],seasonal=True,trace=False,error_action='ignore',suppress_warnings=True,stepwise=True)
    yhat = stepwise_fit.predict(n_periods=horizon)
    y=np.array(df_test['y'])
    forecasts=pd.concat([forecasts,pd.DataFrame({'y':y,'yhat':yhat,'k':k})],axis=0)
result['auto.arima']={}
result['auto.arima']['mape']=np.mean(np.abs((forecasts['y']-forecasts['yhat'])/forecasts['y']))## 0.285
result['auto.arima']['rmse']=np.sqrt(np.sum((forecasts['y']-forecasts['yhat'])**2)/len(forecasts['y'])) ## 166943
print(result['auto.arima'])
forecasts_arima=forecasts.copy()#備份用

看了下，模型的MAPE是 0.285 ，RMSE 高達166943. 可能是由於數據過長，再加上沒有調參，所以導致嚴重過擬合了。我也不打算深究了，感興趣的可以細致的調一下。

3.2 Prophet 模型

facebook 在設計 prophet 之初，綜合考慮了與時序預測相關的主要因素，具體包含：

• 季節性趨勢，具體包含 yearly、weekly等
• 節假日和特殊事件的影響
• changepoint
• outlier

所以在商業上的應用場景會比較廣，據說對於日數據的效果還不錯。從大框架上講，其是一個加法模型（經過log處理也可以變成一個乘法模型），一個時序會被分解為三部分之和

\[y(t) = g(t) +s(t) +h(t)+\varepsilon_t \]

其中g(t)是趨勢項， s(t) 是周期項，如 weekly 和 yearly 等，h(t)是節假日趨勢

對於趨勢項 g(t) 部分，可以由帶changepoint的線性模型或者logistic growth (高中學的物種增長S型曲線？)模型來擬合：

\[\text{logistic:}\quad g(t)=\frac{C(t)}{1+exp(-k(t-m))} \]

\[\text{linear:}\quad g(t)=kt+m \]

其中 C 是carrying capacity，指明環境能容許的最大值，比如手機銷售量最大也不會超過人口太多。k 是增長率。m是offset parameter。當加了changepoint 影響后，方程會有一些變化，具體可參見論文。

對於周期項 s(t) 部分，可以由Fourier series 來逼近

\[s(t)=\sum_{n=1}^{N}\big( a_{n}\cos(\frac{2\pi nt}{P}+b_{n}\sin(\frac{2\pi nt}{P})) \big) \]

這里參數P的設置很有意思，facebook沒有粗暴的直接用n，而是分別對yearly data，weekly data等作擬合。例如當P=7時， 序列 s(t) 的周期有 {7/n}

對於 節假日/事件 項 h(t) 部分，可以由示性函數來逼近

\[h(t)=Z(t)\cdot\boldsymbol{\kappa}=[\boldsymbol{1}(t\in D_1),\cdots,\boldsymbol{1}(t \in D_{L})]\cdot \boldsymbol{\kappa} \]

其中 D_i 是指節假日i的日期集合，Z(t) 是一個矩陣，\kappa 是一個正態分布，擬合每一個節假日的影響程度

prophet 有很多參數可調，這里我們先都用默認的。

forecasts=pd.DataFrame(columns=['y','yhat','k'])

horizon,initial,period=30,1500,30
df_train_list,df_test_list=ts_model_selection.train_test_split(\
    df,horizon='{} days'.format(horizon),initial='{} days'.format(initial),period='{} days'.format(period))
for k in range(len(df_train_list)):
    df_train,df_test=df_train_list[k],df_test_list[k]
    m=Prophet()
    m.fit(df_train[['ds','logy']].rename(columns={'logy':'y'}))
    yhat=np.array(m.predict(df_test[['ds']])['yhat'])
    y=np.array(df_test['logy'])
    y,yhat=np.exp(y),np.exp(yhat)
    forecasts=pd.concat([forecasts,pd.DataFrame({'y':y,'yhat':yhat,'k':k})],axis=0)

result['prophet']={}
result['prophet']['mape']=np.mean(np.abs((forecasts['y']-forecasts['yhat'])/forecasts['y'])) #平均MAPE ,0.126
result['prophet']['rmse']=np.sqrt(np.sum((forecasts['y']-forecasts['yhat'])**2)/len(forecasts['y'])) ##137291
print(result['prophet'])
forecasts_prophet=forecasts.copy()#備份用

這里可以看到我把Y做了log化處理，這是因為log后，相當於把加法模型變成了乘法模型，對於我用的這份數據，乘法模型會更好一些。

可以看到模型的效果是 MAPE =0.13, RMSE=137293, 相比ARIMA要好很多了。

3.3 其他模型

除下上述兩個，還有一些傳統的可以用來預測時序的模型。本文就不再一一細講，感興趣的可以去看看R中的forecast包及其文檔，介紹的很詳細。

找了很久都沒找到好的python包，所以這里提供一個在python中調用R中的ETS模型的code。

# 不保證這段代碼能在任何環境中運行，需要大家自己去調
import os
os.environ['R_HOME'] = 'C:\Program Files\R\R-3.5.0'
os.environ['R_USER'] = 'C:\Anaconda3\Lib\site-packages\rpy2'

import rpy2.robjects as robjects
from rpy2.robjects.packages import importr
ts=robjects.r('ts')
# import R's "base" package
base = importr('base')

# import R's "utils" package
utils = importr('utils')
forecast=importr('forecast',lib_loc = "C:/Users/gason/Documents/R/win-library/3.5")

from rpy2.robjects import pandas2ri
pandas2ri.activate()

df_train=df_train_list[1].loc[:,['y']]
rdata=ts(df_train)

horizon,initial,period=30,1500,30
df_train_list,df_test_list=ts_model_selection.train_test_split(\
    df,horizon='{} days'.format(horizon),initial='{} days'.format(initial),period='{} days'.format(period))

forecasts=pd.DataFrame(columns=['y','yhat','k'])
for k in range(len(df_train_list)):
    #print(k)
    df_train,df_test=df_train_list[k],df_test_list[k]
    rdata=ts(df_train.loc[:,'y'])
    yhat=forecast.forecast(forecast.ets(rdata),h=horizon)
    yhat = np.asarray(yhat[1])
    y=np.array(df_test['y'])
    forecasts=pd.concat([forecasts,pd.DataFrame({'y':y,'yhat':yhat,'k':k})],axis=0)

result['ets']={}
result['ets']['mape']=np.mean(np.abs((forecasts['y']-forecasts['yhat'])/forecasts['y']))
result['ets']['rmse']=np.sqrt(np.sum((forecasts['y']-forecasts['yhat'])**2)/len(forecasts['y']))
print(result['ets'])

4. prphet 的調參

prophet 可調的參數有很多，個人比較比較調的有：

• holidays: 需要注意的是，不是所有的節假日或時間都需要標記，以及windows的設置很考究；
• growth: 可選 linear 模型或者 logistic模型。個人意見，數據長度不長的話，直接用linear就好
• seasonality_prior_scale: 季節性趨勢的強弱，默認為10
• holidays_prior_scale: 節假日趨勢的強弱，默認為10
• changepoint_prior_scale:節changepoit的強弱，默認為0.05

對於分類模型，我們可以直接用sklearn的grid_search來調參，但對於時序，由於交叉驗證的問題，我們沒辦法直接用。這里提供了一個函數供參考。

from fbprophet import Prophet
from fbprophet.diagnostics import cross_validation
from sklearn.model_selection import ParameterGrid
from sklearn.model_selection import ParameterSampler

def ts_evaluation(df,param,horizon=30,period=120,initial=1095,exp=True):
    '''
    利用交叉驗證評估效果   
    '''

    #param={'holidays':holidays,'growth':'linear','seasonality_prior_scale':50,'holidays_prior_scale':20}
    m=Prophet(**param)
    m.fit(df)
    forecasts=m.predict(df[['ds']])
    forecasts['y']=df['y']
    df_cv = cross_validation(m, horizon='{} days'.format(horizon), period='{} days'.format(period), initial='{} days'.format(initial))
    if exp:
        df_cv['yhat']=np.exp(df_cv['yhat'])
        df_cv['y']=np.exp(df_cv['y'])
    mape=np.mean(np.abs((df_cv['y'] - df_cv['yhat']) / df_cv['y']))
    rmse=np.sqrt(np.mean((df_cv['y'] - df_cv['yhat'])**2))
    scores={'mape':mape,'rmse':rmse}
    return scores

def ts_grid_search(df,holidays,param_grid=None,cv_param=None,RandomizedSearch=True,random_state=None):
    '''網格搜索
    時間序列需要特殊的交叉驗證

    df:   
    holidays: 需要實現調好  

    '''

    df=df.copy()
    if param_grid is None:
        param_grid={'growth':['linear']
        ,'seasonality_prior_scale':np.round(np.logspace(0,2.2,10))
        ,'holidays_prior_scale':np.round(np.logspace(0,2.2,10))
        ,'changepoint_prior_scale':[0.005,0.01,0.02,0.03,0.05,0.008,0.10,0.13,0.16,0.2]
        }



    if RandomizedSearch:
        param_list=list(ParameterSampler(param_grid,n_iter=10,random_state=random_state))
    else:
        param_list=list(ParameterGrid(param_grid))

    if cv_param is None:
        cv_param={'horizon':30,'period':120,'initial':1095}


    scores=[]
    for i,param in enumerate(param_list):
        print('{}/{}:'.format(i,len(param_list)),param)
        param.update({'holidays':holidays})
        scores_tmp=ts_evaluation(df,param,exp=True,**cv_param)        
        param.pop('holidays')
        tmp=param.copy()
        tmp.update({'mape':scores_tmp['mape'],'rmse':scores_tmp['rmse']})
        scores.append(tmp)                       
        print('mape : {:.5f}%'.format(100*scores_tmp['mape']))


    scores=pd.DataFrame(scores)

    best_param_=scores.loc[scores['mape'].argmin(),:].to_dict()
    best_scores_=best_param_['mape']
    best_param_.pop('mape')
    best_param_.pop('rmse')


    return best_param_,best_scores_,scores

當然不是說上來直接暴力調參就好，我看了下，對於這份數據而言，changepoint_prior_scale 的影響最大，在一定范圍內，越大效果越好。

下面附上我最后的效果。

---