最近剛好要用到時間序列，所以就開始整理相關資訊。其實想要做預測，在網路上可以找到很多種機器學習的方法做預測，但是要考慮預測的準確性，就是要考慮根據資料的特性來選擇合適的演算法。時間序列簡單來說就是資料序列的變動是跟時間有關係的序列，最明顯的資料序列就是像股票(很多時間序列預測的大神都是在做程式交易的)。其實他主要的原理就是線性回歸(Regression)模型和相關性(Correlation)分析的組合應用(這邊不講理論，只介紹原理概念)。這邊只先討論僅使用該時間序列本身的先前值來預測未來值的這種單變量時間序列預測 - 就是基於自迴歸(Autoregression)分析與自相關(Autocorrelation)分析建立的預測模型。
傳統的時間序列預測方法通常是針對每月，每季，或者年度的資料作分析，但隨著AI的技術運用越來越廣泛，很多時候可能會需要每分鐘或者每小時的即時資料分析方式，傳統的演算法就不一定適用，就需要做調整。因此了解原理概念是很重要的，因為當你知道它演算法原理的時候，你就會知道它使用的侷限性，當遇到問題的時候，你才能知道怎麼調整。
1. 時間序列的基本概念一般來說，時間序列在概念上可拆解成四個項目: (1)趨勢(Trend)，(2)循環性(Cycle)，(3)季節性(Seasonal)，(4)不規則項(Irregular)。所謂的Trend 就是反映一個時間序列的長期演變方向；Cycle則是時間序列所表現的週期性波動(通常是指資料因為商業周期而產生的至少2年的長期震盪)，Seasonal是反應時間序列在不同年份的相同季節性表現的週期性變化；不規則項則是其他無法解釋的誤差或者是隨機變數。傳統時間序列常常是在針對以月或季為單位的經濟統計資料作預測，而這種經濟統計資料很容易有因為受到氣候，風俗習慣或假期等所謂的季節性變動。這些季節變動會增加對資料長期趨勢或循環變動的判讀困難度。因此需要將季節性變動的影響加以分離，再使用一般的預測模型(如ARIMA或SARIMA)來做預測，使資料能更準確定反應長期的趨勢。
常見的時間序列分解方式有以下兩種:
1-1. 加法模型(Additive Model):
- Y(t)=Trend(t)+Cycle(t)+Seasonal(t)+Irregular(t)
- 也有人將Trend(t)和Cycle(t)合在一起看，寫成Y(t)=Trend(t)+Seasonal(t)+Irregular(t)
- 適用於季節性的規模基本保持不變，不隨原始序列水平增減而變化(如圖1(a))。
1-2. 乘法模型(Multiplicative Model):
- Y(t)=Trend(t) x Cycle(t) x Seasonal(t) x Irregular(t)
- 也有人Trend(t)和Cycle(t)合在一起看，寫成Y(t)=Trend(t)+X Seasonal(t) X Irregular(t)
- 適用於季節性的規模會隨著原始序列水平成比例變化(如圖1(b))。
圖2(a)和(b)是典型的加法模型時間序列和乘法模型時間序列的透過python seasonal_decompose()作季節性變動分解的範例。常見的季節性變動分解方法如seasonal_decompose()，STL Decomposition或者由美國普查局和西班牙合作開發的X13-ARIMA-SEATS，其中STL Decomposition可以針對不同的時間粒度的資料作分析，而X13-ARIMA-SEATS只能針對每月或每季的資料作分析。圖3是 X13-ARIMA-SEATS季節性調整流程。
2. 資料的預處理
一般的機器學習方法都會先針對有缺失的資料去補齊資料，但時間序列的資料除了要補缺失的資料以外，針對離群值的資料或者工作日/假日對資料的干擾(及所謂的日曆效應)也要先做移除。從圖3，我們可以窺見X13-ARIMA-SEATS在時間序列的資料做預測之前，針對資料做了那些處理。X13-ARIMA-SEATS有提供回歸變數(Regressors)來將這些干擾做排除。
2-1. 離群值
一般的離群值種類分類如下圖5。因為傳統的預測模型(如ARIMA)沒有辦法分離這些離群值，因此這些離群值必須要在進行預測前先排除，才不會影響模型的建立。
2-2. 日曆效應(Calender Effect)
就是因為假日移動或者交易日數的多寡造成時間序列的波動。分為
- 移動假日(Moving-holiday)效應: 如中秋節落在9月的第3季或者10月的第3季，對資料造成的波動影響。
- 交易日數(Trading day)效應: 每月的資料中，有些周末日數僅8個，有些周末日數有10個，還有閏年平年對資料的波動影響等。
3. 經典預測模型 - ARIMA和SARIMA
3-1. ARIMA模型
ARIMA模型(Autoregressive Integrated Moving Average model)是經典的時間序列預測方法之一。ARIMA(p, d, q)中，AR是自回歸(Autoregression)-就是資料本身跟過去歷史紀錄做回歸分析，MA是移動平均(Moving Average)。
- p是自回歸項數: 表示y的當前值跟前p個歷史值有關係。p值通常是透過偏自相關PACF(Partial Autocorrelation Function)圖來分析決定。
- q是移動平均項數: 表示y的當前值和前q個歷史值AR預測誤差有關係。q值通常是透過ACF(Autocorrelation Function)圖來分析決定。
- d是讓序列平穩的最小差分階數: 非平穩序列可以透過差分得到平穩序列，但過度差分，會導致時間序列失去自相關性。差分的階數對應於數列趨勢的程度，第一階差分代表線性趨勢，第二階差分代表二次趨勢，依此類推。
圖5就是ARIMA的基本預測分析流程。一般找參數的做法大概就是肉眼觀察資料趨勢，用統計的檢驗技術去檢測參數是否合理，用Grid Search(就是暴力搜尋)去尋找合理的參數。一般來說，當資料沒有明顯的季節性，我們就會使用ARIMA。其中圖5的虛線方框的部分也可以用statsmodels所提供的arma_order_select_ic()或者pmdarima.arima.auto_arima()來尋找合適的p和q。不過這些函式庫找到的參數是否合理，就見仁見智。
3-1-1. 定態(或稱平穩性)(Stationary)
如圖5所示，在使用時間序列模型第一件事情就是要確認資料是否屬於訂態(Stationary)資料。所謂的定態資料就是指數列的平均數跟變異數不會隨時間有顯著差異。如果不是定態資料，這個預測模型所估計出來的結果就會不準確。因為時間序列模型的原理就是把資料設定是一個隨時間變動的線性回歸模型，因此需要使用很多統計的方法來檢驗模型是否符合線性回歸的基本假設。ADF檢驗(Augmented Dickey-Fuller Test)就是一個常見的定態性檢驗方式(又稱單根檢定(Unit Root Test))，一般來說，如果計算出來的p值小於0.05就是定態資料，如果大於0.05就是非定態資料。
以圖6(a)的時間序列為例，因為該時間序列沒有季節性(週期)的影響，因此使用ARIMA模型來做預測，我們將該資料皆割成兩份，前120筆資料拿來做訓練模型，後30筆資料拿來驗證預測的準確性(習慣上都是取用觀察值1/4長度的預測值去做預測分析，越長的預測值就越不準確)。我們用ADF測試去檢驗該時間序列，發現p值<0.05，是屬於定態的資料，因此不需要再做差分，故d=0。
3-1-2. 用ACF和PACF找p和q
接著將資料以ACF和PACF作分析如圖6(b)和圖6(c)。ACF是考慮資料所有週期的相關性，而PACF是只考慮特定周期的相關性。圖6(a)中的X軸就是該時間序列從現在跟過去歷史時間點的時間差距，通常以LAG稱之，LAG=1就代表跟現在的時間點差距一個時間單位(假設這是一個每月統計的時間序列，LAG=1就是跟現在的時間序列差距一個月)。藍色區間是誤差範圍，在藍色區間內的資料被認為不顯著，可忽略不計，習慣上都是選最後一個超過藍色區間的數值。
圖6(a)的ACF中，最後一個大於藍色區間藍柱的X值為4，故取q=4；而PACF中最後一個大於藍色區間藍柱的X為3，故取p=3，因此ARIMA(p,d,q)為ARIMA(3,0,4)。不過用ACF跟PACF去找p和q很容易誤判，因為還要考慮圖形式截尾(類似圖6(b))還是拖尾(類似圖6(c))，比較好的方法還是用AIC或者BIC等訊息準則來判斷。所以也有人是多取幾個超過藍色區間的數值當作備選，舉例如(2,0,3),(2,0,4),(3,0,3)和(3,0,4)，再利用尋找最小AIC的方式來選擇最好的(p,d,q)值來做預測。
3-1-2. 差分
圖7(a)的時間序列用ADF測試去檢驗該時間序列，發現p值>0.05，是屬於非定態的時間序列，因此需要再做差分如圖7(b)，直到ADF檢測其p值<0.05為止。經測試發現做了一階差分之後的資料就會變成訂態。所謂的一階差分就是如圖7(c)，把現在t時刻的時間序列跟前一個時間t-1的時間序列相減得到的差值。接著對一階差分的時間序列做ACF和PACF圖形分析，取幾個超過藍色區間的數值當作備選。舉例如圖7(d)可能的q為2,3,4，圖7(e)可能的p為1,2，再用AIC來做判斷。這邊你可能覺得圖7(d)中的10,11,12,13也可以考慮，為什麼不選，因為當ACF表現為10階以上，就可能是時間序列不平穩，就需要再做差分(雖然ADF檢測明明p值<0.5)。因為用ACF和PACF觀察p和q值很容易有爭議，我就比較傾向使用statsmodels所提供的arma_order_select_ic()函示庫來提供幾個可能的p和q，再用AIC準則做篩選。
3-2. SARIMA模型
SARIMA就是ARIMA的延伸，SARIMA的模型可以寫作SARIMA(p,d,q) (P,D,Q)m或者是SARIMA(p,d,q)(P,D,Q)s，也有人稱之為乘積ARIMA模型，其中(p,d,q)就是ARIMA的(p,d,q)，(P,D,Q)m考慮的是季節性(Seasonal)。可以使用美國普查局(Bureau of Census)所提供的X-13ARIMA-SEATS工具來尋找合適的參數，但X13是針對以月分或者是每季為單位的資料。
3-2-1.用ACF和PACF尋找(p,d,q)(P,D,Q)m
圖8是傳統SARIMA尋找(p,d,q)(P,D,Q)m的傳統步驟。其實從原始的時間序列(如圖8(a)跟他對應的ACF圖8(b)就可以發現到這個波的週期(m)是20，因為ACF在LAG=0，LAG=20，LAG=40等地方都有局部最大值，即原始的時間序列在t=0個資料點跟t=20的資料點是重疊的，且在t=20個資料點跟t=40的資料點是重疊的。而且用ADF檢查，發現他不是一個定態時間序列，這個時候就表示需要用季節性差分來消除季節的影響(Seasonal Effect)，才方便觀察參數。
季節性差分的做法就是如圖8(c)先做一階差分(就是把第t時間點的資料跟第t-1點時間的資料相減，即LAG=1個差分)，把線性增長的趨勢消除(差分後變成水平常數行的趨勢)，圖8(d)是一階差分的ACF，這個時候可以很明顯發現LAG=20仍有很明顯的PEAK。接著再做20階差分(把第t時間點的資料旱第t-20時間點的資料相減)，這個時候的時間序列看起來就會如圖8(e)，變成一個定態時間序列(時間序列資料的趨勢變成水平)，這個時候就可以使用圖8(f)的ACF和圖8(g)的PACF來尋找p,q,P,和Q。
從圖8(f)中，跟ARIMA一樣的方式找q，在LAG<20以內，最後一個大於藍色區域的點是2，因此q=2；而LAG=20的時候剛好有一個PEAK，因此Q=0(剛好在週期點上)。從圖9(g)中，在LAG<20以內，最後一個大於藍色區域的點是1，因此p=1；而LAG=21的時候剛好有一個PEAK，因此P=0(剛好在週期+1的點上)。而且在整個找參數的過程中，因為我們先做了一階差分(d=1)讓資料變平穩趨勢，再做一個12階差分(D=1)去消除季節性的影響。所以整個SARIMA模型可能的參數為(1,1,2)(1,1,0)20(用ACF和PACF找出的只是可能的參數組合，而最佳的組合還是要用AIC來做判斷比較好)。其實我個人比較推薦使用grid search的方式，以最小AIC值來尋找最佳的(p,d,q)(P,D,Q)(雖然暴力破解挺花時間的，但在那邊畫圖觀察參數也很麻煩，如果圖形一複雜，其實很難用畫圖分析；我都是先用畫圖找幾個可能的參數解，再混合grid search來找最佳的模型)。
3-2-2. 殘差值自相關性檢測。
圖9(a)就是SARIMA預測模型提供的分析結果，紅色方框就是他的AIC數值(不過因為我只有測試一組預測模型的參數，這個參數不一定是最好的參數)。在使用預測模型的時候，除了AIC值來判斷以外，也需要檢查模型的殘差值，殘查值就是模型擬合值(Fitted Value)跟實際觀察值的差，合理的殘差值應該看起來要像是一個高斯分布的白噪音(White Noise)。圖9(b)就是他的殘差值(Residual)分析結果，從圖10(b)的左上圖是說明這個殘差值均值為0，右上圖是說明這個殘差值是符合常態分佈，左下圖是說明殘差值接近對角分布的直線，右下圖是說明殘差值沒有明顯自相關(就是除了LAG=0以外，沒有數值超過藍色區域)，符合一個白噪音的標準。
圖10就是SARIMA預測模型所產生的預測結果，這邊可以發現，雖然經過了重重的檢查，但是產生的預測值跟實際的觀察值仍然還是有誤差。不過，這個結果可能是由很多因素造成，有可能是觀察值的資料點數不夠(因為我只拿了120個資料點作訓練模型，理論上是要越多點越好，一般都是300個點以上)；也有可能是預測模型不好，因為我並不是拿AIC數字最小的模型。但還是可以發現他的預測其實還是有把過去資料的波形展現出來，只是剛好當下的波形並未跟隨著過去的波動，這也剛好說明不管你怎麼預測，現實狀況還是有不預期的結果發生。雖然預測未必會準確，但是預測還是提供了一定的參考資訊，讓使用者可以做提前準備。
參考資料
- https://www.ibm.com/docs/zh-tw/spss-modeler/SaaS?topic=SS3RA7_sub/modeler_mainhelp_client_ddita/clementine/timeseries_arima_criteria.html
- https://facebook.github.io/prophet/
- https://github.com/linkedin/greykitehttps://otexts.com/fpp2/intro.html
- http://biostat.tmu.edu.tw/oldFile/enews/ep_download/22rb.pdf
- https://ws.ndc.gov.tw/Download.ashx?u=LzAwMS9hZG1pbmlzdHJhdG9yLzEwL3JlbGZpbGUvMC8xMDczOC9lZDVjMzNjZS0yNmIxLTRjY2ItYTI0ZC1jZjcyYzA3M2I4ZGQucGRm&n=MeiHuueBo%2BaZr%2Bawo%2BaMh%2BaomeWto%2BevgOiqv%2BaVtOaWueazleS5i%2BeglOaekC5wZGY%3D&icon=..pdf
- https://www.stat.gov.tw/public/Attachment/8725164742ECW2CXEJ.pdf
- https://www.ibm.com/docs/zh-tw/spss-modeler/SaaS?topic=SS3RA7_sub/modeler_mainhelp_client_ddita/clementine/timeseries_arima_criteria.html
- https://ithelp.ithome.com.tw/articles/10252815
- https://www.allenchi.cn/notes/sarima/
- https://www.jianshu.com/p/f9e4cfc69e12
- https://www.math.pku.edu.cn/teachers/lidf/course/fts/ftsnotes/html/_ftsnotes/fts-seas.html
