時間序列（time series）是一系列有序的資料。通常是等時間間隔的取樣資料。如果不是等間隔，則一般會標註每個資料點的時間刻度。

下面以time series 普遍使用的資料 airline passenger為例。 這是十一年的每月乘客數量，單位是千人次。

如果想嘗試其他的資料集，可以訪問這裡： https://datamarket.com/data/list/?q=provider:tsdl

可以很明顯的看出，airline passenger的資料是很有規律的。

time series data mining 主要包括decompose（分析資料的各個成分，例如趨勢，週期性），prediction（預測未來的值），classification（對有序資料序列的feature提取與分類），clustering（相似數列聚類）等。

這篇文章主要討論prediction（forecast，預測）問題。 即已知歷史的資料，如何準確預測未來的資料。

先從簡單的方法說起。給定一個時間序列，要預測下一個的值是多少，最簡單的思路是什麼呢？

（1）mean（平均值）：未來值是歷史值的平均。

（2）exponential smoothing （指數衰減）：當去平均值得時候，每個歷史點的權值可以不一樣。最自然的就是越近的點賦予越大的權重。

或者，更方便的寫法，用變數頭上加個尖角表示估計值

（3) snaive ： 假設已知資料的週期，那麼就用前一個週期對應的時刻作為下一個週期對應時刻的預測值

（4） drift：飄移，即用最後一個點的值加上資料的平均趨勢

介紹完最簡單的演算法，下面開始介紹兩個time series裡面最火的兩個強大的演算法： Holt-Winters 和 ARIMA。 上面簡答的演算法都是這兩個演算法的某種特例。

（5）Holt-Winters： 三階指數平滑

Holt-Winters的思想是把資料分解成三個成分：平均水平（level），趨勢（trend），週期性（seasonality）。R裡面一個簡單的函式stl就可以把原始資料進行分解：

一階Holt—Winters假設資料是stationary的（靜態分佈），即是普通的指數平滑。二階演算法假設資料有一個趨勢，這個趨勢可以是加性的(additive,線性趨勢)，也可以是乘性的(multiplicative,非線性趨勢)，只是公式裡面一個小小的不同而已。 三階演算法在二階的假設基礎上，多了一個週期性的成分。同樣這個週期性成分可以是additive和multiplicative的。 舉個例子，如果每個二月的人數都比往年增加1000人，這就是additive；如果每個二月的人數都比往年增加120%，那麼就是multiplicative。

R裡面有Holt-Winters的實現，現在就可以用它來試試效果了。我用前十年的資料去預測最後一年的資料。 效能衡量採用的是RMSE。 當然也可以採用別的metrics：

預測結果如下：

結果還是很不錯的。

（6） ARIMA： AutoRegressive Integrated Moving Average

ARIMA是兩個演算法的結合：AR和MA。其公式如下：

是白噪聲，均值為0， C是常數。 ARIMA的前半部分就是Autoregressive：

， 後半部分是moving average：

。 AR實際上就是一個無限脈衝響應濾波器（infinite impulse resopnse）, MA是一個有限脈衝響應（finite impulse resopnse），輸入是白噪聲。

ARIMA裡面的I指Integrated（差分）。 ARIMA（p,d,q）就表示p階AR，d次差分，q階MA。 為什麼要進行差分呢？ ARIMA的前提是資料是stationary的，也就是說統計特性（mean，variance，correlation等）不會隨著時間視窗的不同而變化。用數學表示就是聯合分佈相同：

當然很多時候並不符合這個要求，例如這裡的airline passenger資料。有很多方式對原始資料進行變換可以使之stationary：

（1） 差分，即Integrated。 例如一階差分是把原數列每一項減去前一項的值。二階差分是一階差分基礎上再來一次差分。這是最推薦的做法

（2）先用某種函式大致擬合原始資料，再用ARIMA處理剩餘量。例如，先用一條直線擬合airline passenger的趨勢，於是原始資料就變成了每個資料點離這條直線的偏移。再用ARIMA去擬合這些偏移量。

（3）對原始資料取log或者開根號。這對variance不是常數的很有效。

如何看資料是不是stationary呢？這裡就要用到兩個很常用的量了： ACF（auto correlation function）和PACF(patial auto correlation function)。對於non-stationary的資料，ACF圖不會趨向於0，或者趨向0的速度很慢。 下面是三張ACF圖，分別對應原始資料，一階差分原始資料，去除週期性的一階差分資料：

確保stationary之後，下面就要確定p和q的值了。定這兩個值還是要看ACF和PACF：

確定好p和q之後，就可以呼叫R裡面的arime函數了。 值得一提的是，R裡面有兩個很強大的函式： ets 和 auto.arima。 使用者什麼都不需要做，這兩個函式會自動挑選一個最恰當的演算法去分析資料。

在R中各個演算法的效果如下：

程式碼如下：

passenger = read.csv('passenger.csv',header=F,sep=' ')
p<-unlist(passenger)
pt<-ts(p,frequency=12,start=2001)
plot(pt)
train<-window(pt,start=2001,end=2011+11/12)test<-window(pt,start=2012)

library(forecast)
pred_meanf<-meanf(train,h=12)
rmse(test,pred_meanf$mean) #226.2657pred_naive<-naive(train,h=12)
rmse(pred_naive$mean,test)#102.9765pred_snaive<-snaive(train,h=12)
rmse(pred_snaive$mean,test)#50.70832pred_rwf<-rwf(train,h=12, drift=T)
rmse(pred_rwf$mean,test)#92.66636pred_ses <- ses(train,h=12,initial='simple',alpha=0.2)
rmse(pred_ses$mean,test) #89.77035pred_holt<-holt(train,h=12,damped=F,initial="simple",beta=0.65)
rmse(pred_holt$mean,test)#76.86677  without beta=0.65 it would be 84.41239pred_hw<-hw(train,h=12,seasonal='multiplicative')
rmse(pred_hw$mean,test)#16.36156fit<-ets(train)
accuracy(predict(fit,12),test) #24.390252pred_stlf<-stlf(train)
rmse(pred_stlf$mean,test)#22.07215plot(stl(train,s.window="periodic"))  #Seasonal Decomposition of Time Series by Loess

fit<-auto.arima(train)
accuracy(forecast(fit,h=12),test) #23.538735ma = arima(train, order = c(0, 1, 3),   seasonal=list(order=c(0,1,3), period=12))
p<-predict(ma,12)
accuracy(p$pred,test)  #18.55567BT = Box.test(ma$residuals, lag=30, type = "Ljung-Box", fitdf=2)