8時(shí)間序列回歸模型-R實(shí)現(xiàn)

..時(shí)間序列回歸模型干預(yù)分析概念及模型Box和Tiao引入的干預(yù)分析提供了對(duì)于干預(yù)影響時(shí)間序列的效果進(jìn)行評(píng)估的一個(gè)框架,假設(shè)干預(yù)是可以通過(guò)時(shí)間序列的均值函數(shù)或者趨勢(shì)而對(duì)過(guò)程施加影響,干預(yù)可以自然產(chǎn)生也可以人為施加的,如國(guó)家的宏觀調(diào)控等。其模型可以如下表示：其中mt代表均值的變化,Nt是ARIMA過(guò)程。干預(yù)的分類(lèi)階梯響應(yīng)干預(yù)脈沖響應(yīng)干預(yù)干預(yù)的實(shí)例分析模型初探對(duì)數(shù)化航空客運(yùn)里程的干預(yù)模型的估計(jì)>data<airmiles>>acf<as.vector<diff<diff<window<log<airmiles>,end=c<2001,8>>,12>>>,lag.max=48>#用window得到在911事件以前的未愛(ài)干預(yù)的時(shí)間序列子集對(duì)暫用的模型進(jìn)行診斷>fitmode<-arima<airmiles,order=c<0,1,1>,seasonal=list<order=c<0,1,0>>>>tsdiag<fitmode>從診斷圖可以看出存在三個(gè)異常點(diǎn),acf在12階存在高度相關(guān)因此在季節(jié)中加入MA〔1系數(shù)。擬合帶有干預(yù)信息的模型函數(shù)：arimax<x,order=c<0,0,0>,seasonal=list<order=c<0,0,0>,period=NA>,xreg=NULL,include.mean=TRUE,transform.pars=TRUE,fixed=NULL,init=NULL,method=c<"CSS-ML","ML","CSS">,n.cond,optim.control=list<>,kappa=1e+06,io=NULL,xtransf,transfer=NULL>arimax函數(shù)擴(kuò)展了arima函數(shù),可以處理時(shí)間序列中干擾分析及異常值。假設(shè)干擾影響過(guò)程的均值,相對(duì)未受干擾的無(wú)價(jià)值函數(shù)的偏離用一些協(xié)變量的ARMA濾波器的輸出這種來(lái)表示,偏差被稱(chēng)作傳遞函數(shù)。構(gòu)造傳遞函數(shù)的協(xié)變量通過(guò)xtransf參數(shù)以矩陣或者data.frame的形式代入arimax函數(shù)。air.m1=arimax<log<airmiles>,order=c<0,1,1>,seasonal=list<order=c<0,1,1>,period=12>,xtransf=data.frame<I911=1*<seq<airmiles>==69>,I911=1*<seq<airmiles>==69>>,transfer=list<c<0,0>,c<1,0>>,xreg=data.frame<Dec96=1*<seq<airmiles>==12>,Jan97=1*<seq<airmiles>==13>,Dec02=1*<seq<airmiles>==84>>,method='ML'>>air.m1Call:arimax<x=log<airmiles>,order=c<0,1,1>,seasonal=list<order=c<0,1,1>,period=12>,xreg=data.frame<Dec96=1*<seq<airmiles>==12>,Jan97=1*<seq<airmiles>==13>,Dec02=1*<seq<airmiles>==84>>,method="ML",xtransf=data.frame<I911=1*<seq<airmiles>==69>,I911=1*<seq<airmiles>==69>>,transfer=list<c<0,0>,c<1,0>>>Coefficients:ma1sma1Dec96Jan97Dec02I911-MA0I911.1-AR1I911.1-MA0-0.3825-0.64990.0989-0.06900.0810-0.09490.8139-0.2715s.e.0.09260.11890.02280.02180.02020.04620.09780.0439sigma^2estimatedas0.0006721:loglikelihood=219.99,aic=-423.98畫(huà)圖plot<log<airmiles>,ylab="log<airmiles>">points<fitted<air.m1>>Nine11p=1*<seq<airmiles>==69>plot<ts<Nine11p*<-0.0949>+filter<Nine11p,filter=.8139,method='recursive',side=1>*<-0.2715>,frequency=12,start=1996>,type='h',ylab='9/11Effects'>abline<h=0>從上圖可以看出在20XX底后,911事件的影響效應(yīng)才平息,航班客運(yùn)量恢復(fù)了正常。異常值在時(shí)間序列中異常有兩種,可加異常和新息異常,分別記AO和IO。異常值示例模擬數(shù)據(jù)模擬一般的ARIMA〔1,0,1,然后故意將第10個(gè)觀測(cè)值變成異常值10.>set.seed<12345>>y=arima.sim<model=list<ar=0.8,ma=0.5>,n.start=158,n=100>>yTimeSeries:Start=1End=100Frequency=1[1]0.49180881-0.22323665-0.99151270-0.73387818-0.67750094-1.14472133-2.14844671-2.49530794[9]-1.50355358-2.12615253-0.556517130.413263440.518691291.862106052.199354722.60210165[17]0.791300030.262654262.934148573.990458893.608226781.17845765-0.87682948-1.20637799[25]-1.39501221-0.188321711.229998271.468148502.666474913.234174692.603496241.49513215[33]1.488521420.957392191.300116541.734440532.848251033.732146554.235794563.37049790[41]2.027839551.41218929-0.29974176-1.58712591-1.340808780.107476091.446510811.67809487[49]-0.34663129-0.502914590.01739605-0.014264740.942172040.39046221-0.398835301.60638918[57]1.706682011.375181941.918245340.14254056-2.88169481-3.30372327-1.74068408-3.24868057[65]-3.89415683-3.45920240-1.110420780.679597440.670510840.443940611.895360602.36063873[73]2.005594430.864433240.468475720.723384981.602150981.259222771.531808590.96289779[81]1.077121881.423863540.56318008-0.46689543-0.91861106-1.92947085-2.18188785-1.02759087[89]2.310882723.138473193.012378813.434548072.315394942.449098732.915891411.12648908[97]-0.081238710.444125790.26116418-0.45815484>y[10]<-10模型初步判斷>acf<y>>pacf<y>>eacf<y>AR/MA0123456789101112130xxoooooooooooo1oooooooooooooo2oooooooooooooo3oxoooooooooooo4oxoooooooooooo5xxoooooooooooo6xooooooooooooo7oxoooooooooooo從三個(gè)的結(jié)果來(lái)看,可以初步分析y是AR〔1模型對(duì)模型時(shí)行擬合>m1=arima<y,order=c<1,0,0>>>m1Call:arima<x=y,order=c<1,0,0>>Coefficients:ar1intercept0.54190.7096s.e.0.08310.3603對(duì)模擬模型進(jìn)行異常值探測(cè)>detectAO<m1>[,1][,2][,3]ind9.00000010.00000011.000000lambda2-4.0184129.068982-4.247367>detectAO<m1,robust=F>[,1]ind10.000000lambda27.321709>detectIO<m1>[,1][,2]ind10.00000011.00000lambda17.782013-4.67421AO探測(cè)結(jié)果認(rèn)為第9,10,11.可能出現(xiàn)異常值。IO探測(cè)認(rèn)為第10,11可能出現(xiàn)了異常值。由于檢驗(yàn)統(tǒng)計(jì)量的最大取值出現(xiàn)在10且AO〉IO,所以更認(rèn)為出現(xiàn)異常值在第10是AO異?？紤]異常值的時(shí)間序列擬合>m2=arima<y,order=c<1,0,0>,xreg=data.frame<AO=seq<y>==10>>>m2Call:arima<x=y,order=c<1,0,0>,xreg=data.frame<AO=seq<y>==10>>Coefficients:ar1interceptAO0.80720.569810.9940s.e.0.05700.51290.8012sigma^2estimatedas1.059:loglikelihood=-145.29,aic=296.58>detectAO<m2>[1]"NoAOdetected">detectIO<m2>[1]"NoIOdetected"比較有無(wú)異常值的兩模型再次進(jìn)行異常值探測(cè)時(shí),沒(méi)有發(fā)現(xiàn)異常值,驗(yàn)證最初序列異常出現(xiàn)在10的猜測(cè)對(duì)比模型1和2的擬合效果>tsdiag<m2>>tsdiag<m1>雖然模型二的殘差通過(guò)引入異常值后正太性是顯性的,但是其acf和P值結(jié)果顯示引入MA〔1是必要的。重新擬合適當(dāng)模型>m3=arima<y,order=c<1,0,1>,xreg=data.frame<AO=seq<y>==10>>>detectAO<m3>[1]"NoAOdetected">detectIO<m3>[1]"NoIOdetected">tsdiag<m3>>m3Call:arima<x=y,order=c<1,0,1>,xreg=data.frame<AO=seq<y>==10>>Coefficients:ar1ma1interceptAO0.65960.61540.585011.1781s.e.0.07990.07960.41320.4755sigma^2estimatedas0.793:loglikelihood=-131.16,aic=270.33模型的擬合效果是顯著提高。Acf和P值檢驗(yàn)也一步通過(guò)。>plot<y,type='b'>>arrows<40,7,11,9.8,length=0.8,angle=30>另一個(gè)現(xiàn)實(shí)例子數(shù)據(jù)包中的co2>m1.co2=arima<co2,order=c<0,1,1>,seasonal=list<order=c<0,1,1>,period=12>>>m1.co2Call:arima<x=co2,order=c<0,1,1>,seasonal=list<order=c<0,1,1>,period=12>>Coefficients:ma1sma1-0.5792-0.8206s.e.0.07910.1137sigma^2estimatedas0.5446:loglikelihood=-139.54,aic=283.08>detectAO<m1.co2>[1]"NoAOdetected">detectIO<m1.co2>[,1]ind57.000000lambda13.752715擬合含有新息異常的模型>m4.co2=arimax<co2,order=c<0,1,1>,seasonal=list<order=c<0,1,1>,period=12>,io=c<57>>>m4.co2Call:arimax<x=co2,order=c<0,1,1>,seasonal=list<order=c<0,1,1>,period=12>,io=c<57>>Coefficients:ma1sma1IO-57-0.5925-0.82742.6770s.e.0.07750.10160.7246sigma^2estimatedas0.4869:loglikelihood=-133.08,aic=272.16模型顯示AIC相比之前模型一更小了。而且IO效應(yīng)的P值=2.677/0.7246是顯著的.偽相關(guān)在時(shí)間序列中引入?yún)f(xié)變量,如非洲牧草產(chǎn)量通常與某些氣候指標(biāo)密切相關(guān),在這種發(fā)問(wèn)下在通過(guò)在時(shí)間序列模型中納入相關(guān)的協(xié)變量,將有助于更好的了解基礎(chǔ)過(guò)程以及得到更為準(zhǔn)確的預(yù)測(cè)。模擬數(shù)據(jù)set.seed<12345>X=rnorm<105>Y=zlag<X,2>+.5*rnorm<105>X=ts<X[-<1:5>],start=1,freq=1>Y=ts<Y[-<1:5>],start=1,freq=1>ccf<X,Y,ylab='CCF'>從ccf中可以看出兩樣本在滯后2期存在明顯的相關(guān)性。奶產(chǎn)量與對(duì)數(shù)化發(fā)電量的偽相關(guān)data<milk>data<electricity>milk.electricity=ersect<milk,log<electricity>>#intersect函數(shù)將多個(gè)時(shí)間序列合并在一個(gè)容器中。ccf<as.numeric<milk.electricity[,1]>,as.numeric<milk.electricity[,2]>,main='milk&electricity',ylab='CCF'>兩者相關(guān)性似乎非常的強(qiáng),但實(shí)際上這是因?yàn)樗麄兊母髯源嬖诤軓?qiáng)的自相關(guān)性。預(yù)白化與隨機(jī)回歸對(duì)于具有強(qiáng)自相關(guān)的數(shù)據(jù)而言,很難評(píng)估兩個(gè)過(guò)程之前是否存在依賴(lài)關(guān)系,因而,宜將x和y之間的線(xiàn)性關(guān)系關(guān)聯(lián)從其各自相關(guān)關(guān)系中剝離出來(lái)。預(yù)白化正是為了達(dá)到此目的的一個(gè)有效工具。牛奶與電量的CCF預(yù)白化校正>data<milk>>me.dif=ersect<diff<diff<milk,12>>,diff<diff<log<electricity>,12>>>>prewhiten<as.vector<me.dif[,1]>,as.vector<me.dif[,2]>,ylab='CCf'>再次分析兩者的相關(guān)性,此時(shí)除了時(shí)滯-3具有邊緣顯著外,其他地方?jīng)]有一個(gè)相關(guān)系數(shù)是顯著的?；蟿?dòng)防震這給出的35個(gè)樣本互相關(guān)系婁中大約會(huì)出現(xiàn)1.75=35x0.05個(gè)虛假警報(bào),即這個(gè)-3系數(shù)的顯著可能就是一個(gè)虛假的信息。因此,牛奶與耗電量序列實(shí)際上是基本不相關(guān)的。從而認(rèn)為之前在原始數(shù)據(jù)序列中發(fā)現(xiàn)的強(qiáng)互相關(guān)是偽相關(guān)的。Log〔銷(xiāo)售量與價(jià)格數(shù)據(jù)的相關(guān)性分析預(yù)白化處理plot<bluebird,yax.flip=T>#畫(huà)兩者的時(shí)間序列對(duì)比圖預(yù)白化處理prewhiten<y=diff<bluebird>[,1],x=diff<bluebird>[,2],ylab='ccf'>從CCF圖可以看出兩者之間只在時(shí)滯0處是顯著的。即價(jià)格與銷(xiāo)售量之間存在著很強(qiáng)的同期負(fù)相關(guān)關(guān)系。即當(dāng)期提高價(jià)格將導(dǎo)致銷(xiāo)售量的當(dāng)期下降。一般線(xiàn)性回歸分析>sales=bluebird[,1]>price=bluebird[,2]>chip.m1=lm<sales~price>>summary<chip.m1>Call:lm<formula=sales~price>Residuals:Min1QMedian3QMax-0.54950-0.123730.006670.131360.45170Coefficients:EstimateStd.ErrortvaluePr<>|t|><Intercept>15.8900.21773.22<2e-16***price-2.4890.126-19.75<2e-16***---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1Residualstandarderror:0.188on102degreesoffreedomMultipleR-squared:0.7926, AdjustedR-squared:0.7906F-statistic:389.9on1and102DF,p-value:<2.2e-16>acf<residuals<chip.m1>,ci.type='ma'>由于回歸后的殘差自相關(guān)在四階是顯著的,因此我們要對(duì)其進(jìn)行再一步的分析>eacf<residuals<chip.m1>>AR/MA0123456789101112130xxxxooxxoooooo1xooxoooooooooo2xxoxoooooooooo3xxoxoooooooooo4oxxooooooooooo5xxxoxooooooooo6xxoxxxoooooooo7xoxoooooooooooEacf推薦其殘差包含一個(gè)以〔1,4為頂點(diǎn)為的零值三角形,從而表明其為arma〔1,4模型,因此可將對(duì)數(shù)化銷(xiāo)售量擬合成對(duì)于價(jià)格序列的帶有ARMA〔1,4誤差的回歸模型。模擬ARMA〔1,4初探>chip.m2=arima<sales,order=c<1,0,4>,xreg=data.frame<price>>>chip.m2Call:arima<x=sales,order=c<1,0,4>,xreg=data.frame<price>>Coefficients:ar1ma1ma2ma3ma4interceptprice0.1989-0.05540.25210.07350.526915.7792-2.4234s.e.0.18430.16600.08650.10840.13760.21660.1247sigma^2estimatedas0.02556:loglikelihood=42.35,aic=-70.69結(jié)果表明ma1,ma3的系數(shù)并不顯著,即可認(rèn)為其系數(shù)為0調(diào)整模型>chip.m3=arima<sales,order=c<1,0,4>,xreg=data.frame<price>,fixed=c<NA,0,NA,0,NA,NA,NA>>#第一個(gè)NA指代AR1的系數(shù),第一個(gè)0指ma1第二個(gè)NA指的是ma2第二個(gè)0指的是ma3的系數(shù)。第三個(gè)na指ma4,倒數(shù)第二個(gè)na是指截距項(xiàng)對(duì)應(yīng)的系數(shù),最后一個(gè)na指的是price對(duì)應(yīng)的系數(shù)。>chip.m3Call:arima<x=sales,order=c<1,0,4>,xreg=data.frame<price>,fixed=c<NA,0,N