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

1、時(shí)間序列回歸模型1 干預(yù)分析1.1 概念及模型Box和Tiao引入的干預(yù)分析提供了對于干預(yù)影響時(shí)間序列的效果進(jìn)行評估的一個(gè)框架，假設(shè)干預(yù)是可以通過時(shí)間序列的均值函數(shù)或者趨勢而對過程施加影響，干預(yù)可以自然產(chǎn)生也可以人為施加的，如國家的宏觀調(diào)控等。其模型可以如下表示：其中mt代表均值的變化，Nt是ARIMA過程。1.2 干預(yù)的分類階梯響應(yīng)干預(yù)脈沖響應(yīng)干預(yù)1.3 干預(yù)的實(shí)例分析1.3.1 模型初探對數(shù)化航空客運(yùn)里程的干預(yù)模型的估計(jì)> data(airmiles)> acf(as.vector(diff(diff(window(log(airmiles),end=c(2001,8),12)

2、,lag.max=48)#用window得到在911事件以前的未愛干預(yù)的時(shí)間序列子集對暫用的模型進(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ù)。1.3.2 擬合帶有干預(yù)信息的模型函數(shù)：arimax(x, order = c(0, 0, 0), seasonal = list(order = c(0, 0, 0), period = NA), xreg = NU

3、LL, 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è)干擾影響過程的均值，相對未受干擾的無價(jià)值函數(shù)的偏離用一些協(xié)變量的ARMA濾波

4、器的輸出這種來表示，偏差被稱作傳遞函數(shù)。構(gòu)造傳遞函數(shù)的協(xié)變量通過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

5、=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),

6、 method = "ML", xtransf = data.frame(I911 = 1 * (seq(airmiles) = 69), I911 = 1 * (seq(airmiles) = 69), transfer = list(c(0, 0), c(1, 0)Coefficients: ma1 sma1 Dec96 Jan97 Dec02 I911-MA0 I911.1-AR1 I911.1-MA0 -0.3825 -0.6499 0.0989 -0.0690 0.0810 -0.0949 0.8139 -0.2715s.e. 0.0926 0.1189 0.02

7、28 0.0218 0.0202 0.0462 0.0978 0.0439sigma2 estimated as 0.0006721: log likelihood = 219.99, aic = -423.98畫圖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.

8、2715),frequency=12,start=1996),type='h',ylab='9/11 Effects')abline(h=0)從上圖可以看出在2003年底后，911事件的影響效應(yīng)才平息，航班客運(yùn)量恢復(fù)了正常。2 異常值在時(shí)間序列中異常有兩種，可加異常和新息異常，分別記AO和IO。2.1 異常值示例2.1.1 模擬數(shù)據(jù)模擬一般的ARIMA（1，0，1），然后故意將第10個(gè)觀測值變成異常值10.> set.seed(12345)> y=arima.sim(model=list(ar=0.8,ma=0.5),n.start=158,n=10

9、0)> yTime Series:Start = 1 End = 100 Frequency = 1 1 0.49180881 -0.22323665 -0.99151270 -0.73387818 -0.67750094 -1.14472133 -2.14844671 -2.49530794 9 -1.50355358 -2.12615253 -0.55651713 0.41326344 0.51869129 1.86210605 2.19935472 2.60210165 17 0.79130003 0.26265426 2.93414857 3.99045889 3.6082267

10、8 1.17845765 -0.87682948 -1.20637799 25 -1.39501221 -0.18832171 1.22999827 1.46814850 2.66647491 3.23417469 2.60349624 1.49513215 33 1.48852142 0.95739219 1.30011654 1.73444053 2.84825103 3.73214655 4.23579456 3.37049790 41 2.02783955 1.41218929 -0.29974176 -1.58712591 -1.34080878 0.10747609 1.44651

11、081 1.67809487 49 -0.34663129 -0.50291459 0.01739605 -0.01426474 0.94217204 0.39046221 -0.39883530 1.60638918 57 1.70668201 1.37518194 1.91824534 0.14254056 -2.88169481 -3.30372327 -1.74068408 -3.24868057 65 -3.89415683 -3.45920240 -1.11042078 0.67959744 0.67051084 0.44394061 1.89536060 2.36063873 7

12、3 2.00559443 0.86443324 0.46847572 0.72338498 1.60215098 1.25922277 1.53180859 0.96289779 81 1.07712188 1.42386354 0.56318008 -0.46689543 -0.91861106 -1.92947085 -2.18188785 -1.02759087 89 2.31088272 3.13847319 3.01237881 3.43454807 2.31539494 2.44909873 2.91589141 1.12648908 97 -0.08123871 0.444125

13、79 0.26116418 -0.45815484> y10<-102.1.2 模型初步判斷> acf(y)> pacf(y)> eacf(y)AR/MA 0 1 2 3 4 5 6 7 8 9 10 11 12 130 x x o o o o o o o o o o o o 1 o o o o o o o o o o o o o o 2 o o o o o o o o o o o o o o 3 o x o o o o o o o o o o o o 4 o x o o o o o o o o o o o o 5 x x o o o o o o o o o o

14、o o 6 x o o o o o o o o o o o o o 7 o x o o o o o o o o o o o o 從三個(gè)的結(jié)果來看，可以初步分析y是AR（1）模型2.1.3 對模型時(shí)行擬合> m1=arima(y,order=c(1,0,0)> m1Call:arima(x = y, order = c(1, 0, 0)Coefficients: ar1 intercept 0.5419 0.7096s.e. 0.0831 0.36032.1.4 對模擬模型進(jìn)行異常值探測> detectAO(m1) ,1 ,2 ,3ind 9.000000 10.000000

15、11.000000lambda2 -4.018412 9.068982 -4.247367> detectAO(m1,robust=F) ,1ind 10.000000lambda2 7.321709> detectIO(m1) ,1 ,2ind 10.000000 11.00000lambda1 7.782013 -4.67421AO探測結(jié)果認(rèn)為第9，10，11.可能出現(xiàn)異常值。IO探測認(rèn)為第10，11可能出現(xiàn)了異常值。由于檢驗(yàn)統(tǒng)計(jì)量的最大取值出現(xiàn)在10且AOIO，所以更認(rèn)為出現(xiàn)異常值在第10是AO異常2.1.5 考慮異常值的時(shí)間序列擬合> m2=arima(y,order

16、=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: ar1 intercept AO 0.8072 0.5698 10.9940s.e. 0.0570 0.5129 0.8012sigma2 estimated as 1.059: log likelihood = -145.29, aic = 296.58> detectAO(m2)1 "No AO detected

17、"> detectIO(m2)1 "No IO detected"2.1.6 比較有無異常值的兩模型再次進(jìn)行異常值探測時(shí)，沒有發(fā)現(xiàn)異常值，驗(yàn)證最初序列異常出現(xiàn)在10的猜測對比模型1和2的擬合效果> tsdiag(m2)> tsdiag(m1)雖然模型二的殘差通過引入異常值后正太性是顯性的，但是其acf和P值結(jié)果顯示引入MA（1）是必要的。2.1.7 重新擬合適當(dāng)模型> m3=arima(y,order=c(1,0,1),xreg=data.frame(AO=seq(y)=10)> detectAO(m3)1 "No AO d

18、etected"> detectIO(m3)1 "No IO detected"> tsdiag(m3)> m3Call:arima(x = y, order = c(1, 0, 1), xreg = data.frame(AO = seq(y) = 10)Coefficients: ar1 ma1 intercept AO 0.6596 0.6154 0.5850 11.1781s.e. 0.0799 0.0796 0.4132 0.4755sigma2 estimated as 0.793: log likelihood = -131.16,

19、 aic = 270.33模型的擬合效果是顯著提高。Acf和P 值檢驗(yàn)也一步通過。> plot(y,type='b')> arrows(40,7,11,9.8,length=0.8,angle=30)2.2 另一個(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), peri

20、od = 12)Coefficients: ma1 sma1 -0.5792 -0.8206s.e. 0.0791 0.1137sigma2 estimated as 0.5446: log likelihood = -139.54, aic = 283.08> detectAO(m1.co2)1 "No AO detected"> detectIO(m1.co2) ,1ind 57.000000lambda1 3.752715擬合含有新息異常的模型> m4.co2=arimax(co2,order=c(0,1,1),seasonal=list(order

21、=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: ma1 sma1 IO-57 -0.5925 -0.8274 2.6770s.e. 0.0775 0.1016 0.7246sigma2 estimated as 0.4869: log likelihood = -133.08, aic = 272.16模型顯示AIC相比之前模型一更小了

22、。而且IO效應(yīng)的P 值=2.677/0.7246是顯著的.3 偽相關(guān)在時(shí)間序列中引入?yún)f(xié)變量，如非洲牧草產(chǎn)量通常與某些氣候指標(biāo)密切相關(guān)，在這種發(fā)問下在通過在時(shí)間序列模型中納入相關(guān)的協(xié)變量，將有助于更好的了解基礎(chǔ)過程以及得到更為準(zhǔn)確的預(yù)測。3.1 模擬數(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)性。3.2 奶產(chǎn)量與對數(shù)化發(fā)電

23、量的偽相關(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)性。4 預(yù)白化與隨機(jī)回歸對于具有強(qiáng)自相關(guān)的數(shù)據(jù)而言，很難評估兩個(gè)

24、過程之前是否存在依賴關(guān)系，因而，宜將x和y之間的線性關(guān)系關(guān)聯(lián)從其各自相關(guān)關(guān)系中剝離出來。預(yù)白化正是為了達(dá)到此目的的一個(gè)有效工具。4.1 牛奶與電量的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

25、個(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)的。4.2 Log（銷售量）與價(jià)格數(shù)據(jù)的相關(guān)性分析4.2.1 預(yù)白化處理plot(bluebird,yax.flip=T)#畫兩者的時(shí)間序列對比圖預(yù)白化處理prewhiten(y=diff(bluebird),1,x=diff(bluebird),2,ylab='ccf')從CCF圖可以看出兩者之間只在時(shí)滯0處是顯著的。即價(jià)格與銷售量之間存在著很強(qiáng)的同期負(fù)相關(guān)關(guān)系。即當(dāng)期提高價(jià)

26、格將導(dǎo)致銷售量的當(dāng)期下降。4.2.2 一般線性回歸分析> sales=bluebird,1> price=bluebird,2> chip.m1=lm(salesprice)> summary(chip.m1)Call:lm(formula = sales price)Residuals: Min 1Q Median 3Q Max -0.54950 -0.12373 0.00667 0.13136 0.45170 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 15.890 0.217

27、73.22 <2e-16 *price -2.489 0.126 -19.75 <2e-16 *-Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Residual standard error: 0.188 on 102 degrees of freedomMultiple R-squared: 0.7926,Adjusted R-squared: 0.7906 F-statistic: 389.9 on 1 and 102 DF, p-value: < 2.2e-16> acf(residuals(chip.m1),ci.t

28、ype='ma')由于回歸后的殘差自相關(guān)在四階是顯著的，因此我們要對其進(jìn)行再一步的分析> eacf(residuals(chip.m1)AR/MA 0 1 2 3 4 5 6 7 8 9 10 11 12 130 x x x x o o x x o o o o o o 1 x o o x o o o o o o o o o o 2 x x o x o o o o o o o o o o 3 x x o x o o o o o o o o o o 4 o x x o o o o o o o o o o o 5 x x x o x o o o o o o o o o 6 x

29、x o x x x o o o o o o o o 7 x o x o o o o o o o o o o oEacf推薦其殘差包含一個(gè)以（1，4）為頂點(diǎn)為的零值三角形，從而表明其為arma（1，4）模型，因此可將對數(shù)化銷售量擬合成對于價(jià)格序列的帶有ARMA（1，4）誤差的回歸模型。4.2.3 模擬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

30、)Coefficients: ar1 ma1 ma2 ma3 ma4 intercept price 0.1989 -0.0554 0.2521 0.0735 0.5269 15.7792 -2.4234s.e. 0.1843 0.1660 0.0865 0.1084 0.1376 0.2166 0.1247sigma2 estimated as 0.02556: log likelihood = 42.35, aic = -70.69結(jié)果表明ma1,ma3的系數(shù)并不顯著，即可認(rèn)為其系數(shù)為04.2.4 調(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)對應(yīng)的系數(shù)，最后一個(gè)na指的是price對應(yīng)的系數(shù)。> chip.m3Call:arima(x = sale