1. 程式人生 > >時間序列的R語言實現_Part2

時間序列的R語言實現_Part2

這部分是用指數平滑法做的時間序列的R語言實現,建議先看看指數平滑演算法。

用指數平滑做預測

簡單指數平滑(Simple Exponential Smoothing

對可用加性模型描述的,非週期性的時間序列資料,可用簡單指數平滑來做短期的預測。指數平滑是根據平滑常熟α來做的,α取值在0-1的區間上,α越小越接近0,就表示做預測時對近期觀測所取的比重較大。

說明:指數平滑演算法的原理就是利用歷史觀測資料對未來做預測,α的取值決定著對近期和遠期觀測資料所取的權重。詳細的可以去了解該演算法。

下面是倫敦1813年到1912年的降雨量英尺數的時間序列資料:

計算機生成了可選文字: >raln<~Read100SCanltpm只(nh七七p://robjhynd現an.com/t3dldata/hurst/preclp幾.datn,sklp‘工)>rain3erles<一七s(raln,3tart=c(1813))

計算機生成了可選文字: 況異呂呂Sa一JQ)SU一伯」18201840186018801900Time

由圖可以看出,資料隨時間的隨機波動幅度是大致不變的,所以可以說該時間序列是穩定的。在

R中用簡單指數平滑做預測,我們可以用HoltWinters()方法,使用時需要設定兩個引數beta=FALSEgamma=FALSE

說明:HoltWinters有三個類似的引數,alpha,beta,gamma,這三個引數跟上面講到的簡單指數平滑的取值範圍一樣在0-1之間。可以簡單理解alpha是平滑指數,beta是趨勢指數,gamma是季節指數。rainseries時間序列沒有明顯上升或下降的趨勢,也沒有季節性的變化,所以這裡這兩個引數取false

計算機生成了可選文字: ralnserle3forecast3<一HoltWinters(rainserle3,beta'FALSE,qanu皿a=FALSE)rain3erle3foreCa弓t3>>Holt一Wlnter3expanentlal3moo七hinqwlthouttrendandwithaut3ea3onalcompan$C己11:Hal七Wln七ers(x=ralnserles,beta=FA毛SE,qa"lr.a=FA工SE)Smoothinqparalneter3:alpha:0.024工2工5工beta:FA工SEqa刀口毛a:FA工SECOeff1Clents:l,1]己24。67819

結果alpha很接近0,說明預測中對近期觀測資料取值權重較大。結果儲存在rainseriesforecasts這個list變數中,預測結果儲存在這個

list變數的fitted元素中,它的結果可以檢視到。

計算機生成了可選文字: >raln3erlesforecast3$flttedTllneserles:stdl七=181性End=1912Frequency=1xh己七level181性23。5600023。56000181523。6205性23。6205性181623。5780823。57808181723。7629023。76290181823。7601723。76017181923。7630623。76306182023。8269123。82691182123。7990023。79900182223。9893523。98935182323。9862323。98623

在圖中將原始時間序列和新的時間序列對照看:

計算機生成了可選文字: >plo七(ra工n3erlesforeca3七5)

計算機生成了可選文字: Holt·Winters朽Itering的的0門的入0入p山州一U-、paIU田SqO18201840186018801900T1me

黑色線是原始資料,紅色線是預測資料。檢驗結果的準確度,可以用SSE(誤差項平方和)的值來判斷。SSE也是list變數rainseriesforecasts中的一個元素。檢視SSE的值如下:

計算機生成了可選文字: >:aln3e:lesforecast3$SSE[1]1828.855

這個預測結果原始資料對比誤差項平方和是1828.855

上面例子中,HoltWinters()方法預設的預測僅覆蓋有原始資料的那個時間段,也就是1813年到1912年的降水量的時間序列。用R中的forecast包中的forecast.HoltWinters()方法可以來做這個預測。首先安裝forecast

包。安裝方法很簡單就不說了。安裝完成後載入forecast包。

計算機生成了可選文字: >lib:ary("fo:eca3t")載入需要的程輯包:載入程輯包:200Thefolla甘in口object3are址a3kedfram、packaqe:ba3er:a3.Date,a3.Date.nur吐erlc載入需要的程輯包:"imeDa七eThis15fo了eC己3t5.6警告信.息:1:程輯包、foreca:t’是用R版本3.1.2來建造的2:程輯包,zoo’是用R版本3.1.2來建造的3:程輯包、timeDate’是用R版本3.1.2來建造的>rain3erie3foreca3t32<一fo二ecast.HoltWin七e工3(rain,erle3foreca3t3,h,8)>工己InseriesfOleC己3t32POintFOreC己3tLO80Hl80LO95Hl95191324。6781919。1749330。181性516。2616933。09470191任2咬。6781919。1733330.1830516。2592433。0971519152任。6781919。1717330。18喂6516。2567933。09960191624.6781919.1701330.1862516。25性3433。1020性19172咬.6781919.1685330.1878516。2519033。10性49191824。6781919。1669喂30.1894516。2494533。1069性191924。6781919。1653曦30。1910516。2470133。10938192024。6781919。1637喂30。1926516。2445633。11182

上面forecast.HoltWinters()的兩個引數,第一個是上一步rainseriesforecasts用預設的HoltWinters()得到的結果,第二個h是你想要預測的期數,為8,所以得到的預測結果是從1912年往後推的8年即1913-1920。預測結果有5列資料,第一列Forecast是預測值,第二列第三列是80%的置信區間的下限和上限,第四列第五列是95%置信區間的下限和上限。這個預測結果用圖表展示出來如下:

計算機生成了可選文字: FOFeCastSfFOmHOltwinteFS況異帥C叨呂182018401860188019001920

藍線是從1913年到1920年的預測值,顏色較深的部分是預測值80%的置信區間,顏色較淡的藍色是預測值95%的置信區間。對未來資料進行的預測,因為沒有同期觀測資料,所以沒辦法用前面講到的SSE來檢測預測誤差。

這個例子中樣本的預測誤差存在forecast.HoltWinters()方法返回的list變數中的residuals(殘差)元素中。如果做預測模型不可改良,那預測誤差和連續預測結果不相關。也就是說如果預測誤差和預測結果間存在相關性,那所用的簡單指數平滑模型可以用其他預測方法優化。

R中提供了acf()方法可以檢視樣本預測誤差的相關性圖。若要定義我們想要檢視的最大滯後期數,可以定義acf()方法中的lag.max引數。

例如,計算滯後期在1-20時的樣本預測誤差的相關性,如下操作:

計算機生成了可選文字: >acf(rain3erle3foreca3t32Sre31dua13,laq.znax,20)

計算機生成了可選文字: SeriesrainseriesforecastsZ$residualsU-0哎入.0.(日}八之00】月1.a.L05

從上面的相關性圖可以看出來,在滯後期為3時的自相關結果接近意義界限。

說明:

acf()的說明,自相關公式:

計算機生成了可選文字: R(k)=EI(X‘一川(X'+。一屍)1a2

k是在acf方法中定義的lag.max的值,倫敦降雨量預測的例子中,k的取值就是1-20的範圍內。兩條藍色的虛線是意義界限,acf值在這個界限之間,可以認為相關性接近0,可以忽略,如果超出這個範圍就需要另外的判斷了。這個是我自己的理解,僅供參考。如果有問題,歡迎指正!

測試在1-20的延遲期中,是否有意義的非零相關值,我們可以用Ljung-Boxt測試。在R中,用Box.test()的方法。Box.test()方法中的lag引數用來定義我們想要檢視的最大延遲期。用上面倫敦降雨量的預測誤差,作如下操作:

計算機生成了可選文字: >Box.testlraln3erlesforeca3t32零re31dual3,lag=20,type,"Ljunq一Box.)Box一Ljunq七e3tda七a:X一sq林raln3erie3foreca3t32車re3idual3ared=17.任008,df=20,p一value=0.626己

Ljung-Box測試的統計值是17.4p值是0.6,所以樣本預測誤差的非零相關的可能很小。再次確認預測模型不可再改進,檢視預測誤差是不是以均值0和不變方差按正態分佈。

要檢視預測誤差是否有不變方差,可以被預測誤差的結果做一個時序圖:

計算機生成了可選文字: >plo七.ts(rain3eriesforecastsZ零re3idua13)

計算機生成了可選文字: O}的O的·slenp一seJ$閃s】se。。Jo』sa,工J。sule」18201840186018801900T1me

從上圖可以看出,儘管1820-1830年這前十年的波動幅度相比1840-1850期間的波動幅度要小一些,預測誤差的方差在時間上還是大致穩定的。

檢視預測誤差是否按0正態分佈,可以檢視預測誤差的直方圖和以0和相同標準誤差呈正態分佈的曲線圖,兩者對比檢視。還是同一個例子,需要自己寫一個R的方法plotForecastErrors()來實現可實現:

計算機生成了可選文字: >plo七FarecastError3<一func七ion(forecasterror3)+{+奔rnakeahlstoqrarnoftheforeca3七error3:+mybin31ze<一IQR(forecasterrors)/性+my3d<一sd(forecasterror3)+m飛皿In<一mln(forecasterror3)一mysd*5+my幾ax<一max(forecasterrors)+mysd*3+幸qeneratenormallydls七rlbu七edda七awi七hmean0ands七andarddevlationIny$+mynarln<一rnorm(10000,皿ean=O,sd碩y3d)+mymlnZ<一min(mynorm》+my皿axZ<一nax(mynorm)+if(mymlnZ<mymln){mymln<一myminZ}+If(m生皿axZ>m飛皿ax){皿玉皿ax<一mynlaxZ}+奔rnakearedhlstogramoftheforeca3七error3,wl七hthenormallydlstrlbut$+mybin3<一seq(mymin,mymax,myblnsize)+h13t(fa:ecasterrors,col=nredn,freq=FALSE,break3=Inyblns),奔freq=F及工SEen3urestheareaunder七hehis七oqr視=1+希qeneratenormallydis七rlbu七edda七awi七hmean0ands七andarddeviationIny$+myhlst<一hlst(mynorm,plo七=FALSE,break3=mybin3),奔plotthenormalcurvea3abluel工neon七opof七heh13七oqr柳offoreca3t車+palnt3(myhlst$現主ds,myhist$denslty,七yPe'"1",col="blue",lwd=2)+}

上面是plotForecastErrors()方法程式碼,行末$符號表示不換行,#開始的行表示是註釋。下面使用這個方法:

計算機生成了可選文字: >platFareca3tE:二or3(ra工n3er工e3foreca3t32導re31dual3)

計算機生成了可選文字: Histogramofforecasterrors剎}.0800寸0.0必isu。000.0一30一20一1001020foreCasterrorS

圖上能看出來,預測誤差基本上是以0為中心,大致正態分佈的,與正態分佈曲線相比可能有些向右偏,但右偏量相對小。所以,可以說預測誤差是以0為中心呈正態分佈的。

Ljung-Box測試的結果表明預測誤差的非零自相關的結果很少,它的分佈也滿足以0為中心呈正態分佈。至此,我們可以下結論,簡單指數平滑方法足夠用來做倫敦降雨量的預測模型,也不需要進一步的改進了。並且,80%95%的預測區間也是基於上述結論得到,可以說是有效的。

霍特雙引數指數平滑(Holt’s Exponential Smoothing)

霍特指數平滑法,包含兩個引數α和β。α平滑常數,β是趨勢常數。跟HoltWinters方法的引數一樣,這兩個引數的取值範圍也在0-1間。所以對滿足加性模型,存在增加趨勢或者遞減趨勢的非季節性的資料,可以選擇用霍特指數平滑法來做短期的預測。

霍特指數平滑法能預估在當前時間點的水平和趨勢。下面是一個非季節性的存在增減趨勢的時間序列例子。從1866年到1911年每年女性裙子邊緣直徑的資料:

計算機生成了可選文字: >3klrt3<一scan(.http://robjhynd九an.com/七3dldata/raberts/3klrt3.dat",3klp=5)Re己d性6lte幾3>3klrt33erie3<一t3(3kirt3,3tart,c(1866))

計算機生成了可選文字: 00600900卜S山一J山55七一藝S00918701880189019001910Time

可以看出女性裙子邊緣直徑從1866600左右增大到了1880年的1050,在這之後又降到了1911年的520

還是用R中的HoltWinters()方法,這裡我們需要用到alphabeta兩個引數,所以只需要設定gamma=FALSE就行。給女性裙子邊緣直徑的變化這個時間序列做預測模型過程如下:

計算機生成了可選文字: 3klrts3erle3foreca3ts<一Hol七Winte二3(3k工rt33e二ie3,ga刀”毛a=FALSE)SklltSSellesf0工eC己StS>>Holt一WintersexponentialSlnaa七hlnqwl七htrendandwl七hou七3ea3onalcomponen七.Cdll:Hol七Winters(x=sk王r七53erle5,q己nlt己=F幾工SE)Smoo七hingparalne七e了s:alpha:0.8383任8幾betd:1qa刀口厄a:FA乙SECoeff1Clen七呂:[,1]d529。308585bs。690任6任>3ki了七33e二工e3fo工ec己呂七5場55已11]1695弓.18

alphabeta的值分別為0.8381,都很大,說明時間序列水平和趨勢部分的預測值,對近期觀測資料所取的權重較大。這個結果從該時間序列隨時間的水平和趨勢變化都很大,就能很直觀看出來。改時間序列預測的誤差項平方和SSE結果是16954.18

檢視預測結果時間序列圖:

計算機生成了可選文字: 》plat(3klrt33erle3foreca3七3)

計算機生成了可選文字: Holt·Wintersnltering00600卜po州一U-、p山IUOSqO00舊18701880189019001910Time

上圖可以看出,除了預測結果有很小的滯後外,預測值時間序列和實際值序列很接近。

也可以自己設定初始的水平值和趨勢值,HoltWinters()方法中的l.startb.start這兩個引數用來設定初始值。初始水平值一般取第一個時間點的值,初始趨勢值則常取第二個值與第一個值的差值。這個例子中,分別是6089(617-608)

嘗試設定l.startb.start的值,再對女性裙子邊緣直徑時間序列做預測,結果如下,與之前的結果有了一些不同。

計算機生成了可選文字: >HoltWinte工SHOlt一Winte工3(skirt33erie3qanu毛a=FALSE,1.3ta二t=608b.3ta二t=9)exponen七王alsrnoo七hinqwl七h七rendandw工七hou七3ea昌analcalnponen七.Cdll:HaltWinters(x=昌klrt55erles,qan立皿a=FALSE,1.3七ar七=608,b.3tart=9)smao七hlnqparalne七er3:alpha:0.83性6775bet己:工qa刃山衛a:FA毛SECOeffiClentS:[,11d529。278637bs。670129

前面是同期的預測,同樣可以用forecast包中的forecast.HoltWinters()方法對未來做預測。女性裙子邊緣直徑資料是從1866年到1911年的,對此我們可以做從1911年之後的預測。如若想要1912年到193019期的預測資料,如下操作:

計算機生成了可選文字: sklrt33erle3fa了eca3t32<一foreca3t.HoltWinter3(3klrt33erie3foreca3t3,h=19)plat.farecast(3klrt33erlesforeca3t32)>>

計算機生成了可選文字: FOFeCastSfFOmHO!twinteFS00咱}00咱000r.1870188018901900191019201930

與之前用HoltWinters()方法做預測結果一樣,藍線是預測值,深色區域是預測結果的80%置信區間,淺色部分是95%的置信區間。

重複前面用相關性函式來看看是否需要優化模型的過程。

計算機生成了可選文字: >acf(3k工r七33erle3fareca3t32$re3idual3,laq.幾ax=20)>Box.te3t(3klrt33erle3foreca3tsZ$re31dua13,lag=20,type,"Ljunq一Bax")Box一Ljunqte3七data:X一3quSklrt3Serie3fO了eCast32gre3idUa1Sared,19.7312,df,20,p一value,0.4749

計算機生成了可選文字: SeriesskirtsseriesforecastsZ$residualsU-O<咱.0.八U,L15產.飛J‘馬」月..a.』05

在滯後期為5時,預測誤差的自相關結果超過了意義界限,所以,這裡我們需要進一步的檢測。檢視Box-Ljung測試的結果,可以看看到,p值為0.4749,這個值很小,樣本的預測誤差在1-20的滯後期內,存在非零自相關的可能很小。

對簡單指數平滑,還需要檢視預測誤差是不是存在不變方差,並以0為中心呈正態分佈。類似的過程,使用之前的plotForecastErrors()方法。

計算機生成了可選文字: plat.七3(3klrt33eriesforeca3tsZ寫reslduals)plotForeca3tE了二or3(3k工r七33e二工e3fo二ecast32寫re3idual3)>>

計算機生成了可選文字: O寸O入00入·slenpls。J$閃s}se。ejo』s。工J。55七一藝s18701880189019001910Time

計算機生成了可選文字: HistogramOfforecasterrorSO釗0.0OrO.O必一su.0000.0一100一5050foreCasterrorS

從上面兩個圖,能看出,預測誤差的方差隨時間的變化可以認為是大致穩定的,預測誤差也基本滿足正態分佈。所以,結論是,用霍特雙指數平滑模型對從1866年到1911年每年女性裙子邊緣直徑的時間序列做預測是合適的,預測結果經過檢驗也是有效的。

霍特季節性指數平滑(Holt-Winters Exponential Smoothing)

如果時間序列滿足增量模型,存在升降趨勢,並且是季節性的資料。這時,可用霍特季節性指數平滑法來做短期預測。

霍特季節性指數平滑法,包含三個引數α、β和γ。α平滑常數,β是趨勢常數,γ是季節常數。三個引數的取值範圍都是0-1。在R中的實現,還是使用HoltWinters()方法,這一次,它的三個類似引數,我們都需要用到。

使用的時間序列資料是前面取對數後的昆士蘭沙灘旅遊勝地的某一紀念品店的銷售資料。

計算機生成了可選文字: 屍rOr68s。工Ja,s。任工七一ua,,nos口。一19871989199,1993Time

計算機生成了可選文字: souvenir七lmeselie3foreca3t3<一Holtwin七er3(lag3ouvenirtl幾e3erle3)SOUVenlrtlnese了le3fOleCast3>>Holt一Win七ersexponen七ialsmoothinqwithtrendandaddltlveseasonalcolnponen$Cdll:HoltWinte:s(x=loqsauvenlr七1幾e3elle3)smoo七hingalpha:0be七a:Oq己nu肚a:0pdl己現e七eI3:。413418。9561275COeff1Clents:[,11a10。37661961bo。0299631931一0。8Q95206352一0。60576477530。011032385任一0。2任16055155一0。3593351756一0。18076683570。07788605580。101任7055590。096任93535100。051978265110。性17936373121。18088任23

計算機生成了可選文字: >3ouvenlrtl現e3e二le3foreca3t3容SSE11]2.011491

alphabetagamma的值,分別是0.413400.956alpha的值比較小,表明該時間序列的某一時間點的水平預測值,是基於近期觀測值和遠期觀測值。beta0,表明時間序列趨勢部分值不隨時間變化而改變的,也就是所有時間點上,趨勢的預測值都是初始值。gamma0.956,這說明季節性部分的預測值,對近期觀測值所取得權重很大。

檢視預測結果,如下:

計算機生成了可選文字: >plo七<souven工rtl瓦e3e二le3fo二eca3t3)

計算機生成了可選文字: Holt·Wintersnltering屍ro屍68p田劃一U-、paIUasqO1988198919901991199219931994T1me

從圖上,可以看出,霍特季節性指數平滑方法在做這個時間序列的季節性預測時很成功,預測結果與原始資料的11月銷售高峰很貼切。

接下來,對未來做預測,還是用forecast包。

計算機生成了可選文字: 3ouvenir七工瓦e3e二le3fo二eca3t32<一foreca3七.HoltWin七e二3(3auven工rt工幾e3e二王e3fareca3t3,h,曦8)plot.foreca3七(3auvenlrtlme3erle3foreca3t32)>>

計算機生成了可選文字: FOFeCastSfFOmHoltwinteFS只寸}內L囚r}}198819901992199419961998

下面對預測結果做檢驗,看預測模型是否需要改進。同樣的方法計算相關性和做Ljung-Box檢驗。過程及結果如下:

計算機生成了可選文字: >acf(3ouvenirti皿e3erie3foreca3t32車re3idual3,laq.rnax=20)>Box.te3七(呂ouven工r七工現e3er王e3foreca3七32車re3工dua15,la於20,七卯e="Ljunq一Box時)Bax一Ljungte3tdd七d:X一呂qu昌auvenl工七l兀eserlesfo工ec已5t52s二es工du已15ared=17.530任,df=20,p一value=0.6183

計算機生成了可選文字: SeriessouvenirtimeseriesforecastsZ$residuals甲〔〕?O寸.0U-0哎卿〔〕9《習釗.0.0_00_5Lag

樣本的預測誤差的自相關結果,在1-20的滯後期中,沒有超出意義界限。並且Ljung-Boxp值小,說明在1-20的滯後期中,存在非零自相關的可能性很小。

然後檢驗預測誤差是有不變方差,以0為中心呈正態分佈。過程結果如下:

計算機生成了可選文字: plo七.t3(souvenirtine3erie3foreca3t32$re51duals)plo七Foreca3七Error3(3ouvenlr七工meserlesforecas七sZSresiduals)>>

計算機生成了可選文字: 寸.0入.00.0入.0·寸O·竺mnpls。J$閃s}se。。Jo』s。一J。s。任一七一u。,nos1988198919901991199219931994TIme

計算機生成了可選文字: HistogramofforeCasterrorSp介〕吵〔〕9〔〕們.入O入們.}OL必一suao一1_0一05O_005,0foreCasterrors

可看出,預測誤差的方差可認為是隨時間穩定的,其分佈也符合正態分佈。所以這個例子中,使用霍特季節性指數平滑所做的預測結果是有效的。