時(shí)間序列（time series）是一系列有序的數(shù)據(jù)。通常是等時(shí)間間隔的采樣數(shù)據(jù),。如果不是等間隔,，則一般會(huì)標(biāo)注每個(gè)數(shù)據(jù)點(diǎn)的時(shí)間刻度。

time series data mining 主要包括decompose（分析數(shù)據(jù)的各個(gè)成分,，例如趨勢(shì),，周期性），prediction（預(yù)測(cè)未來(lái)的值）,，classification（對(duì)有序數(shù)據(jù)序列的feature提取與分類(lèi)）,，clustering（相似數(shù)列聚類(lèi)）等。

這篇文章主要討論prediction（forecast,，預(yù)測(cè)）問(wèn)題,。 即已知?dú)v史的數(shù)據(jù)，如何準(zhǔn)確預(yù)測(cè)未來(lái)的數(shù)據(jù),。

 下面以time series 普遍使用的數(shù)據(jù) airline passenger為例,。 這是十一年的每月乘客數(shù)量，單位是千人次,。

原始數(shù)據(jù)（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

如果想嘗試其他的數(shù)據(jù)集,，可以訪(fǎng)問(wèn)這里：  https:///data/list/?q=provider:tsdl

可以很明顯的看出，airline passenger的數(shù)據(jù)是很有規(guī)律的,。

先從簡(jiǎn)單的方法說(shuō)起,。給定一個(gè)時(shí)間序列，要預(yù)測(cè)下一個(gè)的值是多少,，最簡(jiǎn)單的思路是什么呢,？

（1）mean（平均值）：未來(lái)值是歷史值的平均,。

（2）exponential smoothing （指數(shù)衰減）：當(dāng)去平均值得時(shí)候，每個(gè)歷史點(diǎn)的權(quán)值可以不一樣,。最自然的就是越近的點(diǎn)賦予越大的權(quán)重,。

或者,，更方便的寫(xiě)法,，用變量頭上加個(gè)尖角表示估計(jì)值

（3) snaive ： 假設(shè)已知數(shù)據(jù)的周期，那么就用前一個(gè)周期對(duì)應(yīng)的時(shí)刻作為下一個(gè)周期對(duì)應(yīng)時(shí)刻的預(yù)測(cè)值

（4） drift：飄移,，即用最后一個(gè)點(diǎn)的值加上數(shù)據(jù)的平均趨勢(shì)

介紹完最簡(jiǎn)單的算法,，下面開(kāi)始介紹兩個(gè)time series里面最火的兩個(gè)強(qiáng)大的算法： Holt-Winters 和 ARIMA。 上面簡(jiǎn)答的算法都是這兩個(gè)算法的某種特例,。

（5）Holt-Winters：  三階指數(shù)平滑

 Holt-Winters的思想是把數(shù)據(jù)分解成三個(gè)成分：平均水平（level）,，趨勢(shì)（trend），周期性（seasonality）,。R里面一個(gè)簡(jiǎn)單的函數(shù)stl就可以把原始數(shù)據(jù)進(jìn)行分解：

一階Holt—Winters假設(shè)數(shù)據(jù)是stationary的（靜態(tài)分布）,，即是普通的指數(shù)平滑。二階算法假設(shè)數(shù)據(jù)有一個(gè)趨勢(shì),，這個(gè)趨勢(shì)可以是加性的(additive,線(xiàn)性趨勢(shì)),，也可以是乘性的(multiplicative,非線(xiàn)性趨勢(shì))，只是公式里面一個(gè)小小的不同而已,。  三階算法在二階的假設(shè)基礎(chǔ)上,，多了一個(gè)周期性的成分。同樣這個(gè)周期性成分可以是additive和multiplicative的,。 舉個(gè)例子,，如果每個(gè)二月的人數(shù)都比往年增加1000人，這就是additive,；如果每個(gè)二月的人數(shù)都比往年增加120%,，那么就是multiplicative。

 R里面有Holt-Winters的實(shí)現(xiàn),，現(xiàn)在就可以用它來(lái)試試效果了,。我用前十年的數(shù)據(jù)去預(yù)測(cè)最后一年的數(shù)據(jù)。 性能衡量采用的是RMSE,。 當(dāng)然也可以采用別的metrics：

預(yù)測(cè)結(jié)果如下：

結(jié)果還是很不錯(cuò)的,。

（6） ARIMA： AutoRegressive Integrated Moving Average

ARIMA是兩個(gè)算法的結(jié)合：AR和MA。其公式如下：

是白噪聲,，均值為0,， C是常數(shù)。  ARIMA的前半部分就是Autoregressive：,， 后半部分是moving average： ,。  AR實(shí)際上就是一個(gè)無(wú)限脈沖響應(yīng)濾波器（infinite impulse resopnse）, MA是一個(gè)有限脈沖響應(yīng)（finite impulse resopnse）,，輸入是白噪聲。

ARIMA里面的I指Integrated（差分）,。 ARIMA（p,d,q）就表示p階AR,，d次差分，q階MA,。  為什么要進(jìn)行差分呢,？ ARIMA的前提是數(shù)據(jù)是stationary的，也就是說(shuō)統(tǒng)計(jì)特性（mean,，variance,，correlation等）不會(huì)隨著時(shí)間窗口的不同而變化。用數(shù)學(xué)表示就是聯(lián)合分布相同：

當(dāng)然很多時(shí)候并不符合這個(gè)要求,，例如這里的airline passenger數(shù)據(jù),。有很多方式對(duì)原始數(shù)據(jù)進(jìn)行變換可以使之stationary：

（1） 差分，即Integrated,。 例如一階差分是把原數(shù)列每一項(xiàng)減去前一項(xiàng)的值,。二階差分是一階差分基礎(chǔ)上再來(lái)一次差分。這是最推薦的做法

（2）先用某種函數(shù)大致擬合原始數(shù)據(jù),，再用ARIMA處理剩余量,。例如，先用一條直線(xiàn)擬合airline passenger的趨勢(shì),，于是原始數(shù)據(jù)就變成了每個(gè)數(shù)據(jù)點(diǎn)離這條直線(xiàn)的偏移,。再用ARIMA去擬合這些偏移量。

（3）對(duì)原始數(shù)據(jù)取log或者開(kāi)根號(hào),。這對(duì)variance不是常數(shù)的很有效,。

如何看數(shù)據(jù)是不是stationary呢？這里就要用到兩個(gè)很常用的量了： ACF（auto correlation function）和PACF(patial auto correlation function),。對(duì)于non-stationary的數(shù)據(jù),，ACF圖不會(huì)趨向于0，或者趨向0的速度很慢,。 下面是三張ACF圖,，分別對(duì)應(yīng)原始數(shù)據(jù)，一階差分原始數(shù)據(jù),，去除周期性的一階差分?jǐn)?shù)據(jù)：

acf(train)
acf(diff(train,lag=1))
acf(diff(diff(train,lag=7)))

確保stationary之后,，下面就要確定p和q的值了。定這兩個(gè)值還是要看ACF和PACF：

確定好p和q之后,，就可以調(diào)用R里面的arime函數(shù)了,。 以上是ARIMA的基本概念，要深究它的話(huà)還是有很多內(nèi)容要補(bǔ)充的。 ARIMA更多表示為  ARIMA(p,d,q)(P,D,Q)[m] 的形式,，其中m指周期（例如7表示按周）,，p,d,q就是前面提的內(nèi)容，P,D,Q是在周期性方面對(duì)應(yīng)的p,d,q含義,。

值得一提的是,，R里面有兩個(gè)很強(qiáng)大的函數(shù)： ets 和 auto.arima。 用戶(hù)什么都不需要做,，這兩個(gè)函數(shù)會(huì)自動(dòng)挑選一個(gè)最恰當(dāng)?shù)乃惴ㄈシ治鰯?shù)據(jù),。

在R中各個(gè)算法的效果如下：

代碼如下：

passenger = read.csv('passenger.csv',header=F,sep=' ')
p<-unlist(passenger)

#把數(shù)據(jù)變成time series。  frequency=12表示以月份為單位的time series. start 表示時(shí)間開(kāi)始點(diǎn),，可以用c(a,b,...)表示,，  例如按月為單位,，標(biāo)準(zhǔn)的做法是 start=c(2011,1) 表示從2011年1月開(kāi)始
#如果要表示按天的,，建議用 ts(p,frequency=7,start=c(1,1))  很多人喜歡用 ts(p,frequency=365,start=(2011,1))但是這樣有個(gè)壞處就是沒(méi)有按星期對(duì)齊
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)

看到有人問(wèn)代碼中的rmse是怎么寫(xiě)的，其實(shí)‘a(chǎn)ccuracy()’ 函數(shù)已經(jīng)包含了各種評(píng)價(jià)指標(biāo)了,。這里貼上自己寫(xiě)的代碼：

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的實(shí)現(xiàn)相比是不是簡(jiǎn)化了很多
    errSum
}

rmse = function(pred,test)
{ 
    res<- sqrt(mean((pred-test)^2) )
    res
}