時間序列（time series）是一系列有序的數據。通常是等時間間隔的采樣數據。如果不是等間隔，則一般會標注每個數據點的時間刻度。
time series data mining 主要包括decompose（分析數據的各個成分，例如趨勢，周期性），prediction（預測未來的值），classification（對有序數據序列的feature提取與分類），clustering（相似數列聚類）等。
這篇文章主要討論prediction（forecast，預測）問題。 即已知歷史的數據，如何准確預測未來的數據。
下面以time series 普遍使用的數據 airline passenger為例。 這是十一年的每月乘客數量，單位是千人次。
原始數據（passenger.csv）：
112 115 145 171 196 204 242 284 315 340 360 417 118 126 150 180 196 188 233 277 301 318 342 391 132 141 178 193 236 235 267 317 356 362 406 419 129 135 163 181 235 227 269 313 348 348 396 461 121 125 172 183 229 234 270 318 355 363 420 472 135 149 178 218 243 264 315 374 422 435 472 535 148 170 199 230 264 302 364 413 465 491 548 622 148 170 199 242 272 293 347 405 467 505 559 606 136 158 184 209 237 259 312 355 404 404 463 508 119 133 162 191 211 229 274 306 347 359 407 461 104 114 146 172 180 203 237 271 305 310 362 390 118 140 166 194 201 229 278 306 336 337 405 432
如果想嘗試其他的數據集，可以訪問這里： https://datamarket.com/data/list/?q=provider:tsdl
可以很明顯的看出，airline passenger的數據是很有規律的。
先從簡單的方法說起。給定一個時間序列，要預測下一個的值是多少，最簡單的思路是什么呢？
（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圖，分別對應原始數據，一階差分原始數據，去除周期性的一階差分數據：
acf(train) acf(diff(train,lag=1)) acf(diff(diff(train,lag=7)))
確保stationary之后，下面就要確定p和q的值了。定這兩個值還是要看ACF和PACF：
確定好p和q之后，就可以調用R里面的arime函數了。 以上是ARIMA的基本概念，要深究它的話還是有很多內容要補充的。 ARIMA更多表示為 ARIMA(p,d,q)(P,D,Q)[m] 的形式，其中m指周期（例如7表示按周），p,d,q就是前面提的內容，P,D,Q是在周期性方面對應的p,d,q含義。
值得一提的是，R里面有兩個很強大的函數： ets 和 auto.arima。 用戶什么都不需要做，這兩個函數會自動挑選一個最恰當的算法去分析數據。
在R中各個算法的效果如下：
代碼如下：
passenger = read.csv('passenger.csv',header=F,sep=' ')
p<-unlist(passenger)
#把數據變成time series。 frequency=12表示以月份為單位的time series. start 表示時間開始點，可以用c(a,b,...)表示， 例如按月為單位，標准的做法是 start=c(2011,1) 表示從2011年1月開始
#如果要表示按天的，建議用 ts(p,frequency=7,start=c(1,1)) 很多人喜歡用 ts(p,frequency=365,start=(2011,1))但是這樣有個壞處就是沒有按星期對齊 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.2657 pred_naive<-naive(train,h=12) rmse(pred_naive$mean,test)#102.9765 pred_snaive<-snaive(train,h=12) rmse(pred_snaive$mean,test)#50.70832 pred_rwf<-rwf(train,h=12, drift=T) rmse(pred_rwf$mean,test)#92.66636 pred_ses <- ses(train,h=12,initial='simple',alpha=0.2) rmse(pred_ses$mean,test) #89.77035 pred_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.41239 pred_hw<-hw(train,h=12,seasonal='multiplicative') rmse(pred_hw$mean,test)#16.36156 fit<-ets(train) accuracy(predict(fit,12),test) #24.390252 pred_stlf<-stlf(train) rmse(pred_stlf$mean,test)#22.07215 plot(stl(train,s.window="periodic")) #Seasonal Decomposition of Time Series by Loess fit<-auto.arima(train) accuracy(forecast(fit,h=12),test) #23.538735 ma = 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.55567 BT = Box.test(ma$residuals, lag=30, type = "Ljung-Box", fitdf=2)
看到有人問代碼中的rmse是怎么寫的，其實‘accuracy()’ 函數已經包含了各種評價指標了。這里貼上自己寫的代碼：
wape = function(pred,test) { len<-length(pred) errSum<-sum(abs(pred[1:len]-test[1:len])) corSum<-sum(test[1:len]) result<-errSum/corSum result } mae = function(pred,test) { errSum<-mean(abs(pred-test)) #注意 和wape的實現相比是不是簡化了很多 errSum } rmse = function(pred,test) { res<- sqrt(mean((pred-test)^2) ) res }