用R分析時間序列(time series)數據

時間序列（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
}