時間序列（一）：上手體驗

我們把按照時間次序排列的隨機變量序列

\[Y_0,\, Y_1,\, Y_2, \cdots \]

稱為時間序列（Time Series）。比如網站的PV、DAU，國家的GDP，股票的價格等。

這種特別的次序給模型提出了特別的挑戰，包含數據內的自相關性、不可交換性、以及數據和參數的不平穩性等。

時間序列里面的內容很多，小到接下來要講的平滑，大到金融里的混沌時間序列預測等。我准備花一些時間認真整理部分理論和應用，第一篇就分享下上手攻略，關鍵詞是 insight。

本文采用的數據集是 Peyton Manning 的 Wikipedia 頁面的 PV，下載地址是：https://github.com/facebook/prophet/tree/master/examples ，也可以在公眾號后台回復 數據集 下載。

df = pd.read_csv('.\\data\\example_wp_peyton_manning.csv')
df['ds']=pd.to_datetime(df['ds'])
df['y']=np.log(df['y'])# tricks
plt.plot(df['ds'],df['y'],'.')
ax.set_title('Raw Time Series')

1. 手動平滑時間序列和獲取趨勢項

假設時間序列有周期為 s 的周期性趨勢，一種最簡單的平滑方式就是移動平均法(moving average)，假設窗口設為7，則

\[s_t=\frac{1}{7}(y_t+y_{t-1}+y_{t-2}+\cdots +y_{t-6}) \]

其實就是股票中的5日均線、20日均線。

y1=pd.rolling_mean(df['y'],window=7)
# 如果要使得平滑后長度一致，可以設置參數,min_periods=1
y2=pd.rolling_mean(df['y'],window=365)
fig,ax=plt.subplots()
ax.plot(df['ds'],df['y'],'.',alpha=0.4,label='Raw')
ax.plot(df['ds'],y1,'-',label='Win:7')
ax.plot(df['ds'],y2,'-',label='Win:365')
ax.legend()
ax.set_title('Smoothing : Moving Average')

移動平均法只利用了前面n天的數據，而且權重一樣。一種更好的方法是指數平滑法（exponentially weighted moving average），越靠近當天的數據權重越高

\[s_0=y_0, \,\, s_t=\alpha y_t +(1-\alpha)y_{t-1} \]

從這個公式看不出指數在哪，但當我們把它展開后就很明顯了。

\[s_t=\alpha [y_t+(1-\alpha)y_{t-1} + \cdots + (1-\alpha)^{t-1}y_1] + (1-\alpha)^{t} y_0 \]

事實上 EMWA 是一種沒有常數項的ARIMA(0,1,1)模型，當把alpha的選取標准設為最小化 s_t 和 y_{t+1} 就是一個很簡單的時間序列預測模型。

股票中的MACD指標利用的就是指數平滑，其中的DIF線就是12日的指數平滑值 減去26日的指數平滑值。

y1=pd.ewma(df['y'],span=7)
y2=pd.ewma(df['y'],span=365)
fig,ax=plt.subplots()
ax.plot(df['ds'],df['y'],'.',alpha=0.4,label='Raw')
ax.plot(df['ds'],y1,'-',label='Span:7')
ax.plot(df['ds'],y2,'-',label='Span:365')
ax.legend()
ax.set_title('Smoothing : Exponentially Weighted Moving Average')

除此之外還有一些更好的方法，大家可以試試，比如 Holt-Winter三次指數平滑法等

2. 利用 Prophet 看趨勢和周期

時間序列經過合理的函數變換后都可以被認為是由三個部分疊加而成。分別是趨勢項部分、周期項部分和噪聲項部分

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

其中 s(t) 表示周期項，如 weekly seasonality（周一和周二是不一樣的）和 yearly seasonality（平時和寒暑假是不一樣的等）等。對於一些特別的場景，比如網站的DAU，還需要考慮節假日成分、特殊時間成分等。

我們可以利用 Facebook 開源的包 Prophet 來分解。

# 可以添加節假日參數，這樣分解會更准確
m=Prophet()
m.fit(df)
future=m.make_future_dataframe(periods=1)
forecast=m.predict(future)
forecast1=forecast.loc[:,['ds','yhat']]
#m.plot(forecast);
fig=m.plot_components(forecast,weekly_start=1)

Prophet 利用加法模型把序列分成了4個部分（如果給定節假日參數，則是5個部分）。我們可以大概看下它的預測效果，之后我准備花一篇文章專門梳理各個方法（包）的預測效果。

df1=pd.merge(forecast.loc[:,['ds','yhat']],df,on='ds',how='inner')
# 計算RMSE
from sklearn import metrics
rmse=np.sqrt(metrics.mean_squared_error(df1['yhat'],df1['y']))
fig,ax=plt.subplots()
ax.plot(df1['ds'],df1['y'],'.',label='Raw')
ax.plot(df1['ds'],df1['yhat'],'-',alpha=0.8,label='Prophet')
ax.legend()
ax.set_title('Prophet Predict (RMSE= {:.2f})'.format(rmse));
fig.savefig('.\\_images\\prophet_predict.png',dpi=500)

PS： 因為我們用的是訓練集，所以這個RMSE並不能用來評估預測的效果。

3. 從頻域看可能存在的周期

我們還可以利用Fourier Transform 在頻域里看看時間序列。給定一個函數 f(x), 則其 傅里葉變換可以表示為：

\[\hat{f}(x)=\int f(x)e^{-2\pi i x \xi} dx \]

w = np.fft.fft(df['y']-y4)
n=len(w)
w=w[1:]
power = np.abs(w[:int(n/2)])
nyquist = 1/2
freq = np.arange(int(n/2))/(n/2)*nyquist
fig,ax=plt.subplots()
ax.plot(freq,power)
ax.set_xlim(0,0.2)
ax.set_title('Power Spectral');
ax.plot(1/365,power[7],'o',color='red')
ax.text(1/365,power[7],'   Hz: 1/365')
ax.plot(1/7,power[416],'o',color='red')
ax.text(1/7,power[416],'   Hz:1/7')
fig.savefig('.\\_images\\power_spectral.png',dpi=500)

功率譜上的兩個峰值就是時間序列存在的周期。

4. 從時頻域看異常值

小波變換可以把時間序列直接在時頻域進行分解成高頻部分和低頻部分。低頻部分包含了序列的大部分信息，高頻部分則包含了一些細節信息。現在的JPEG2000壓縮標准就是基於小波變換設計的，如果將一張圖片進行小波變換，則低頻部分跟原圖差別很小，高頻部分則大概能看出圖片的輪廊，通過一些方法就能進行圖片邊界的檢測了，大家有興趣可以試試。

首先回顧一下小波變換的相關理論. 詳細的我就不講了，比較復雜。設 \(f(x)\in L^2(\mathbb{R})\)，對於任意精度 \(\varepsilon\)，我們都能找到一個 \(j\) 使得

\[\| f-f_j\|_{L^2} \leq \varepsilon,\quad f_j \in V_j, \]

\[f_j(x)=\sum_{k\in \mathbb{Z}}c_{j,k}\varphi_{j,k}(x) \]

其中 \(\{\varphi_{j,k}(x)\}\) 是具有緊支撐的函數族，比如 db1,db2,...

import pywt
(cA,cD)=pywt.dwt(df['y'],'db4')
n=len(df)
cAn=pywt.upcoef('a',cA,'db4',take=n)
cDn=pywt.upcoef('d',cD,'db4',take=n)
fig,[ax0,ax1,ax2]=plt.subplots(3,1)
ax0.plot(df['ds'],df['y'],'.')
ax0.set_xlabel('Raw')

ax1.plot(df['ds'],cAn,'.')
ax1.set_xlabel('Low Frequency')

ax2.plot(df['ds'],cDn,'.')
ax2.set_xlabel('High Frequency')

fig.tight_layout()
#fig.savefig('.\\_images\\小波分解.png',dpi=500)

fig,ax=plt.subplots()
cDn[np.abs(cDn)<3*np.std(cDn)]=0
ax.plot(df['ds'],df['y'],'.')
tmp=df.copy()
tmp['cDn']=cDn
tmp=tmp.loc[tmp['cDn']>0,['ds','y']]
ax.plot(tmp['ds'],tmp['y'],'o',alpha=0.4,color='red')
ax.set_xlabel('outlier mark')
fig.tight_layout()
#fig.savefig('.\\_images\\異常值標記.png',dpi=500)

接下來的計划還沒想好，大概率有一篇講ARIMA系列模型，一篇評估現有各種包的預測效果，一篇講 Change Point ，一篇將 RNN 和 LSTM。

參考文獻

[1]. 指數平滑法, http://connor-johnson.com/2014/02/01/smoothing-with-exponentially-weighted-moving-averages/

[2]. 指數平滑法, https://www.jianshu.com/p/6fb0408b3f54

---