用R語(yǔ)言做時(shí)間序列分析

1、用R語(yǔ)言做時(shí)間序列分析時(shí)間序列（time series ）是一系列有序的數(shù)據(jù)。通常是等時(shí)間間隔的采樣數(shù)據(jù)。如果不是 等間隔， 則一般會(huì)標(biāo)注每個(gè)數(shù)據(jù)點(diǎn)的時(shí)間刻度。下面以time series普遍使用的數(shù)據(jù)airline passenger為例。這是一年的每月乘客數(shù)量，單位是千人次。Time200220042006200820102012Time200220042006200820102012如果想嘗試其他的數(shù)據(jù)集，可以訪問(wèn)這里：/data/list/?q=provider:tsdl可以很明顯的看出，airli ne passe nger的數(shù)據(jù)是很有規(guī)律的。time series data min

2、ing主要包括decompose（分析數(shù)據(jù)的各個(gè)成分，例如I趨勢(shì)，周期性），prediction（預(yù)測(cè)未來(lái)的值），classificatio n（對(duì)有序數(shù)據(jù)序列的feature提取與分類），clusteri ng （相似數(shù)列聚類）等。這篇文章主要討論prediction （ forecast，預(yù)測(cè)）問(wèn)題。即已知?dú)v史的數(shù)據(jù)，如何準(zhǔn)確預(yù)測(cè)未來(lái)的數(shù)據(jù)。先從簡(jiǎn)單的方法說(shuō)起。給定一個(gè)時(shí)間序列，要預(yù)測(cè)下一個(gè)的值是多少，最簡(jiǎn)單的思路是什么呢？（1）mean （平均值）：未來(lái)值是歷史值的平均。Xr+l =（2） exponential smoothing（指數(shù)衰減）：當(dāng)去平均值得時(shí)候，每個(gè)歷史點(diǎn)的權(quán)值可 以不

3、一樣。最自然的就是越近的點(diǎn)賦予越大的權(quán)重。=aX+ + ct% + a3X3 + 或者，更方便的寫法，用變量頭上加個(gè)尖角表示估計(jì)值Xt+1 = aXt + （1 - a）X（3）sn aive :假設(shè)已知數(shù)據(jù)的周期，那么就用前一個(gè)周期對(duì)應(yīng)的時(shí)刻作為下一個(gè)周期對(duì) 應(yīng)時(shí)刻的預(yù)測(cè)值（4）drift :飄移，即用最后一個(gè)點(diǎn)的值加上數(shù)據(jù)的平均趨勢(shì)tXt+h|t =禺+占 2 心-斗-X= xt +占（身如Tt 二 介紹完最簡(jiǎn)單的算法，下面開(kāi)始介紹兩個(gè)time series里面最火的兩個(gè)強(qiáng)大的算法：Holt-Winters 和ARIMA。上面簡(jiǎn)答的算法都是這兩個(gè)算法的某種特例。（5） Holt-Winte

4、rs :三階指數(shù)平滑Holt-Winters的思想是把數(shù)據(jù)分解成三個(gè)成分：平均水平（level ），趨勢(shì)（trend ），周 期性（seasonality ）。R里面一個(gè)簡(jiǎn)單的函數(shù)stl就可以把原始數(shù)據(jù)進(jìn)行分解:200220042006200820102012200220042006200820102012time一階Holt - Winters假設(shè)數(shù)據(jù)是stationary 的（靜態(tài)分布），即是普通的指數(shù)平滑。二階算法假設(shè)數(shù)據(jù)有一個(gè)趨勢(shì)，這個(gè)趨勢(shì)可以是加性的（additive，線性趨勢(shì)），也可以是乘性的（multiplicative,非線性趨勢(shì)），只是公式里面一個(gè)小小的不同而已。三階算法在二

5、階的假設(shè)基礎(chǔ)上，多了一個(gè)周期性的成分。同樣這個(gè)周期性成分可以是additive和multiplicative的。舉個(gè)例子，如果每個(gè)二月的人數(shù)都比往年增加1000人，這就是additive ；如果每個(gè)二月的人數(shù)都比往年增加120%，那么就是multiplicative。Holt-WintersExponential smoothing數(shù))tq=CtZf + ( Q)S,才 A t)Double exponential smoothingsi = tI齡二 art + (1 -+ bi)5 = J 1 工 0 /3(5f -5(_1) + (1 FtA = St + mbtTriple expon

6、ential smoothing死=R。卸二ft + (1 n)(的-i + 如-)L二#(弱一5(_i) + (1 -如-iQ = 7+(i-7)q-l斤 + m =臨 + 77也)G-1+l+g-i) TTind JR里面有Holt-Winters的實(shí)現(xiàn)，現(xiàn)在就可以用它來(lái)試試效果了。我用前十年的數(shù)據(jù)去預(yù)測(cè)最后一年的數(shù)據(jù)。性能衡量采用的是RMSE。當(dāng)然也可以采用別的metrics :ftlMAE =廣 1 工 SYt -ftlt=iMSE =n(yft=iMMSE =n(yft=iMRMSE =、rt=LMAPE = 100n_1 工預(yù)測(cè)結(jié)果如下:rpIM 一 t=iForecasts fr

7、om Holt-Wintersmultiplieath/e method結(jié)果還是很不錯(cuò)的。ARIMA : AutoRegressive Integrated Moving AverageARIMA是兩個(gè)算法的結(jié)合：AR和MA。其公式如下：血=c +iQ-i是白噪聲，均值為0, C是常數(shù)。ARIMA的前半部分就是Autoregressive :_ ，后半部分是moving _ j, t ,_Vs? n-average :_。AR實(shí)際上就是一個(gè)無(wú)限脈沖響應(yīng)濾波器(infinite impulse resopnse ) , MA 是個(gè)有限脈沖響應(yīng)(finite impulse resopnse ),

8、輸入是白噪聲。ARIMA里面的I指Integrated（差分）。ARIMA （ p,d,q ）就表示p階AR, d次差分，q階MA。為什么要進(jìn)行差分呢？ARIMA的前提是數(shù)據(jù)是stationary的，也就是說(shuō)統(tǒng)計(jì)特性（mean，varianee，correlation等）不會(huì)隨著時(shí)間窗口的不同而變化。用數(shù)學(xué)表 示就是聯(lián)合 分布相同：+T? * -1 ”） F（赴 i S * i 龍起開(kāi)當(dāng)然很多時(shí)候并不符合這個(gè)要求，例如這里的airli ne passe nger數(shù)據(jù)。有很多方式對(duì)原始數(shù)據(jù)進(jìn)行變換可以使之 statio nary :（1）差分，即Integrated 。例如一階差分是把原數(shù)列每一

9、項(xiàng)減去前一項(xiàng)的值。二階差 分是一階差分基礎(chǔ)上再來(lái)一次差分。這是最推薦的做法（2） 先用某種函數(shù)大致擬合原始數(shù)據(jù)，再用ARIMA處理剩余量。例如，先用一條直線擬合airline passenger的趨勢(shì)，于是原始數(shù)據(jù)就變成了每個(gè)數(shù)據(jù)點(diǎn)離這條直線的偏移。再用ARIMA去擬合這些偏移量。（3）對(duì)原始數(shù)據(jù)取log或者開(kāi)根號(hào)。這對(duì)varianee不是常數(shù)的很有效。如何看數(shù)據(jù)是不是stationary呢？這里就要用到兩個(gè)很常用的量了：ACF（autoeorrelation function ）禾口 PACF （patial auto eorrelation function）。對(duì)于non-stationa

10、ry的數(shù)據(jù)，ACF圖不會(huì)趨向于0 ,或者趨向0的速度很慢。下面是三張ACF圖，分別對(duì)應(yīng)原始數(shù)據(jù)，一階差分原始數(shù)據(jù)，去除周期性的一階差分?jǐn)?shù)據(jù)：acfLagdiff act00051 015diff act andemave seasonalityLag確保stationary 之后，下面就要確定p和q的值了。定這兩個(gè)值還是要看ACF和PACF :MF形狀才旨數(shù)寰減至0正員交替MF形狀才旨數(shù)寰減至0正員交替，衰減至01個(gè)或者多個(gè)倩大于6其余大致為0 在某個(gè)點(diǎn)之后開(kāi)始指數(shù)衰減至0 所有點(diǎn)都罡0一固定點(diǎn)上有大于。的值不衰減至0ModelAR modek 用 PACF去確認(rèn)p的倩MA model, q等

11、于知:F中大于0的個(gè)數(shù)AR和MA模型的結(jié)合原始數(shù)據(jù)完全隨機(jī)二原始數(shù)據(jù)有周期性二原始數(shù)據(jù)不罡占tmtiongryl確定好p和q之后，就可以調(diào)用R里面的arime函數(shù)了。值得一提的是，R里面有兩個(gè) 很強(qiáng)大的函數(shù)：ets和 auto.arima。用戶什么都不需要做，這兩個(gè)函數(shù)會(huì)自動(dòng)挑選一個(gè) 最恰當(dāng)?shù)乃惴ㄈシ治鰯?shù)據(jù)。在R中各個(gè)算法的效果如下:meanf(train, h=12226.2657naive(tfainf h 12)102.9765snaiveftrain, h-12)50.70832rwfftrairt, h-12, drift=T)92.66636sesftrainA ti=12,ini

12、tial=lsimple,f alpha=0.2)89.77035holtftrain, h=12/damped=F/ initiakAsiimpfe”, beta=0.65)76.86677hw(train, h=12,seasonal-multiplicative1)16.36156ets(train)24390252stlf(train|22.07215auto.arimaftrain)23.538735arimpftrain, order = c(0? 1,3), seasonal=list(order=c(0,lT3), period=12)1835567代碼如下：ISpasseng

13、er.csv ,header=F,sep=)passenger = read.csv( pv-unlist(passenger) pt-ts(p,frequency= plot(pt)train-window(pt,start=12 ,start= 2001 )2012 )passenger = read.csv( pv-unlist(passenger) pt-ts(p,frequency= plot(pt)train-window(pt,start=12 ,start= 2001 )2012 )library(forecast)pred_meanf-meanf(train,h= rmse(

14、 test ,pred_meanf$mean) #12)226.265712)226.2657rmse(pred_naive$mean,pred_snaive-snaive(train,h=rmse(pred_snaive$mean,12)test )# 102.976512) test )# 50. 7083212)test )# 102.976512) test )# 50. 70832rmse(pred_rwf$mean, test )# 92.66636pred_ses - ses(train,h=12 ,initial= simple ,alpha= 0. 2)rmse(pred_s

15、es$mean, test ) # 89.77035rmse(pred_holt$mean.test )# 76.86677 without beta=0. 65 it would be84.41239pred_hw-hw(train,h=rmse(pred_hw$mean12 ,seasonal= multiplicativetest )# 16. 36156fitv-ets(train) accuracy(predict(fit.12), test ) # 24. 390252pred_stlf-stlf(train)rmse(pred_stlf$mean,test )# 22.07215plot(stl(train,s.window=periodic ) #SeasonalDecomposition of Time Seriesby Loess fitv-auto.arima(train)accuracy(forecast(fit,h=12), test ) # 23.538735accuracy(forecast(fit,h=12), test ) # 23.538735ma = arima(train, order =