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

1、時間序列回歸模型1干預(yù)分析1.1概念及模型Box和Tiao引入的干預(yù)分析提供了對于干預(yù)影響時間序列的效果進行評估的一個框架，假設(shè)干預(yù)是可以通過時間序列的均值函數(shù)或者趨勢而對過程施加影響干預(yù)可以自然產(chǎn)生也可 以人為施加的，如國家的宏觀調(diào)控等。其模型可以如下表示:其中mt代表均值的變化，Nt是ARIMA過程。1.2干預(yù)的分類階梯響應(yīng)干預(yù)脈沖響應(yīng)干預(yù)圈鞘1】-4有關(guān)蛛期哺庵干頊的一些常皿程郭(堀帶有一個時間申位的BE5R)1.3干預(yù)的實例分析 模型初探對數(shù)化航空客運里程的干預(yù)模型的估計現(xiàn)在回到每月航空客運里程的數(shù)據(jù).如前所述.，20019 JJ的恐怖襲擊事件使航空客運排徊于蕭條之中，該干預(yù)效應(yīng)可用在

2、明。1年9月有脈沖輸入的AR (D過程來表系*這一意 外事件神航空客運量即時造成了一種強烈的激冷效應(yīng).因此對此干預(yù)致成(第11效應(yīng)建 模如下昭=沒+土五P戶只：中 T代表引年9月.在這一去示中，他代表叩酎的9】效應(yīng).1時, 庭代去9/11娥應(yīng)對其后*個月份所造成的影響這里逐需要魂定基礎(chǔ)天擾過程的季節(jié) ARIMA結(jié)枸.基于珂于頊數(shù)據(jù)，暫用個ARIMA (0, b 1)X(0, 1.田兒模型表示未受擾 的過程！參見圖表115data(airmiles)acf(as.vector(diff(diff(window(log(airmiles),end=c(2001,8),12),lag.max=48)

3、#用window得到在911事件以前的未愛干預(yù)的時間序列子集對暫用的模型進行診斷fitmode tsdiag(fitmode)5IQH從診斷圖可以看出存在三個異常點，acf在12階存在高度相關(guān)因此在季節(jié)中加AMA（ 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(

4、CSS-ML, ML, CSS), n.cond, optim.control = list(),kappa = 1e+06, io = NULL, xtransf, transfer = NULL)arimax函數(shù)擴展了 arima函數(shù)，可以處理時間序列中干擾分析及異常值。假設(shè)干擾影響過 程的均值，相對未受干擾的無價值函數(shù)的偏離用一些協(xié)變量的ARMA濾波器的輸出這種來 表示，偏差被稱作傳遞函數(shù)。構(gòu)造傳遞函數(shù)的協(xié)變量通過xtransf參數(shù)以矩陣或者data.frame的形式代入arimax函數(shù)。air.m1=arimax(log(airmiles),order=c(0,1,1),seasona

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

6、, 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 = lis

7、t(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人2 estimated as 0.0006721: log likelihood = 219.99, aic = -423.98畫圖plot(log(airmiles),ylab=log(airmiles)poin

8、ts(fitted(air.m1)+ (6 寸60.0rdltl uz)smo_d(69R(S 蘭 Esb srIHdtltl-uz s品 QSIa?852日 i ik S 521. F2I-seiuiiiB|OD|8&LttrEO-(OHSu=qe (-sbyzLJJIr(966I0sCNIH&uenbE 、(LnIR.0r(srDE-HpoqE、6rnl8.Hwd y10 acf(y)Lag pacf(y) eacf(y)AR/MA0 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 o1 o o o o o o o o o

9、o o o o oo o o o o o o o o o o o o oo x o o o o o o o o o o o oo x o o o o o o o o o o o ox x o o o o o o o o o o o ox o o o o o o o o o o o o oo x o o o o o o o o o o o o從三個的結(jié)果來看，可以初步分析y是AR (1)模型對模型時行擬合m1=arima(y,order=c(1,0,0)m1Call:arima(x = y, order = c(1, 0, 0)Coefficients:ar1 intercept0.54190.

10、7096s.e. 0.08310.3603對模擬模型進行異常值探測detectAO(m1),1,2,3ind 9.000000 10.000000 11.000000lambda2 -4.018412 9.068982 -4.247367detectAO(m1,robust=F),1ind 10.000000lambda2 7.321709detectIO(m1),1,2ind 10.000000 11.00000lambda1 7.782013 -4.67421AO探測結(jié)果認為第9，10，11.可能出現(xiàn)異常值。IO探測認為第10，11可能出現(xiàn)了異常 值。由于檢驗統(tǒng)計量的最大取值出現(xiàn)在10且A

11、OIO，所以更認為出現(xiàn)異常值在第10是AO異?？紤]異常值的時間序列擬合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:ar1 intercept AO0.80720.5698 10.9940s.e. 0.05700.51290.8012 sigmaA2 estimated as 1.059: log likelihood = -145.29, aic

12、= 296.58detectAO(m2)1 No AO detecteddetectIO(m2)1 No IO detected比較有無異常值的兩模型再次進行異常值探測時，沒有發(fā)現(xiàn)異常值，驗證最初序列異常出現(xiàn)在10的猜測對比模型1和2的擬合效果tsdiag(m2) tsdiag(m1)M澆他ns禰ffiia史AMM面訓(xùn)用汁吊湖削吊昏、B湖Macf *口 P面薜湘削斗史A MA (1)湖&圈g。 SOI晉雌胃m3 Har-ma(ybrderHC(lo、l)xregHdaor.frame(AOHSeq(y)HH 10)de-ec_AO(m3)E-ZO 0 deieciecrde-ec_IO(m3)

13、E=NOIO deieciecrV-l-sd-ag(m3)m3ca=az.maxH y、order H cp。、1)、xreg H daor.frame(AO H seq(y) 10)Coefficients:ar1ma1interceptAO0.65960.61540.585011.1781s.e. 0.07990.07960.41320.4755sigmaA2 estimated as 0.793: log likelihood = -131.16, aic = 270.33模型的擬合效果是顯著提高。Acf和P值檢驗也一步通過。plot(y,type=b)arrows(40,7,11,9.

14、8,length=0.8,angle=30)2.2另一個現(xiàn)實例子數(shù)據(jù)包中的co2m1.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:ma1 sma1-0.5792 -0.8206s.e. 0.07910.1137sigmaA2 estimated as 0.5446: log likelih

15、ood = -139.54, aic = 283.08detectAO(m1.co2)1 No AO detecteddetectIO(m1.co2),1ind 57.000000lambda1 3.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(5

16、7)Coefficients:ma1 sma1 IO-57-0.5925 -0.8274 2.6770s.e. 0.07750.1016 0.7246sigma人2 estimated as 0.4869: log likelihood = -133.08, aic = 272.16模型顯示AIC相比之前模型一更小了。而且9效應(yīng)的P值=2.677/0.7246是顯著的.3偽相關(guān)在時間序列中引入?yún)f(xié)變量，如非洲牧草產(chǎn)量通常與某些氣候指標(biāo)密切相關(guān)，在這種發(fā)問下在 通過在時間序列模型中納入相關(guān)的協(xié)變量，將有助于更好的了解基礎(chǔ)過程以及得到更為準(zhǔn)確 的預(yù)測。3.1模擬數(shù)據(jù)set.seed(12345)X=

17、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)X&Y-15-10-5051015Lag從ccf中可以看出兩樣本在滯后2期存在明顯的相關(guān)性。3.2奶產(chǎn)量與對數(shù)化發(fā)電量的偽相關(guān)data(milk)data(electricity)milk.electricity二ersect(milk,log(electricity)# intersect 函數(shù)將多個時間序列合并 在一個容器中。milk,eiecVicl(y

18、ccf(as.numeric(milk.electricity,1),as.numeric(milk.electricity,2), main = milk & electricity,ylab=CCF)mikk & eh&ctrtclly兩者相關(guān)性似乎非常的強，但實際上這是因為他們的各自存在很強的自相關(guān)性。4預(yù)白化與隨機回歸對于具有強自相關(guān)的數(shù)據(jù)而言，很難評估兩個過程之前是否存在依賴關(guān)系，因而，宜將x 和y之間的線性關(guān)系關(guān)聯(lián)從其各自相關(guān)關(guān)系中剝離出來。預(yù)白化正是為了達到此目的的一個 有效工具。4.1牛奶與電量的CCF預(yù)白化校正data(milk)me.dif=ersect(dif

19、f(diff(milk,12),diff(diff(log(electricity),12)prewhiten(as.vector(me.dif,1),as.vector(me.dif,2),ylab=CCf)再次分析兩者的相關(guān)性，此時除了時滯-3具有邊緣顯著外，其他地方?jīng)]有一個相關(guān)系數(shù)是 顯著的?；蟿臃勒疬@給出的35個樣本互相關(guān)系婁中大約會出現(xiàn)1.75 = 35x0.05個虛假警 報，即這個-3系數(shù)的顯著可能就是一個虛假的信息。因此，牛奶與耗電量序列實際上是基 本不相關(guān)的。從而認為之前在原始數(shù)據(jù)序列中發(fā)現(xiàn)的強互相關(guān)是偽相關(guān)的。4.2 Log (銷售量)與價格數(shù)據(jù)的相關(guān)性分析預(yù)白化處理plot

20、(bluebird,yax.flip=T)#畫兩者的時間序列對比圖UuE-bhdline預(yù)白化處理prewhiten(y=diff(bluebird),1,x=diff(bluebird),2,ylab=ccf)從CCF圖可以看出兩者之間只在時滯0處是顯著的。即價格與銷售量之間存在著很強的同 期負相關(guān)關(guān)系。即當(dāng)期提高價格將導(dǎo)致銷售量的當(dāng)期下降。4.2.2 一般線性回歸分析sales二bluebird,1price二bluebird,2chip.m1 = lm(salesprice)summary(chip.m1)Call:lm(formula = sales price)Residuals:M

21、in1Q Median3QMax-0.54950 -0.12373 0.00667 0.13136 0.45170Coefficients:Estimate Std. Error t value Pr(|t|)(Intercept) 15.8900.21773.222e-16 *price-2.4890.126 -19.752e-16 *Signif. codes: 0 *0.001 * 0.01 * 0.05 . 0.1Residual standard error: 0.188 on 102 degrees of freedomMultiple R-squared: 0.7926,Adju

22、sted R-squared: 0.7906F-statistic: 389.9 on 1 and 102 DF, p-value: acf(residuals(chip.m1),ci.type=ma)由于回歸后的殘差自相關(guān)在四階是顯著的，因此我們要對其進行再一步的分析eacf(residuals(chip.m1)AR/MA0 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 ox o o x o o o o o o o o o ox x o x o o o o o o o o o ox x o x o o o o o o o o

23、 o oo x x o o o o o o o o o o ox x x o x o o o o o o o o ox x o x x x o o o o o o o ox o x o o o o o o o o o o oEacf推薦其殘差包含一個以(1，4)為頂點為的零值三角形，從而表明其為arma(1，4) 模型，因此可將對數(shù)化銷售量擬合成對于價格序列的帶有ARMA(1，4)誤差的回歸模型。模ARMA(1,4)初探chip.m2=arima(sales,order=c(1,0,4),xreg=data.frame(price)chip.m2Call: arima(x = sales, o

24、rder = c(1, 0, 4), xreg = data.frame(price)Coefficients:ar1 ma1 ma2 ma3 ma4 intercept price 0.1989 -0.0554 0.2521 0.0735 0.526915.7792 -2.4234s.e. 0.18430.1660 0.0865 0.1084 0.13760.21660.1247sigma人2 estimated as 0.02556: log likelihood = 42.35, aic = -70.69 結(jié)果表明ma1,ma3的系數(shù)并不顯著，即可認為其系數(shù)為0調(diào)整模型chip.m3=a

25、rima(sales,order=c(1,0,4),xreg=data.frame(price),fixed=c(NA,0,NA,0,N A,NA,NA)#第一個NA指代AR1的系數(shù)，第一個0指ma1第二個NA指的是ma2第二 個0指的是ma3的系數(shù)。第三個na指ma4，倒數(shù)第二個na是指截距項對應(yīng)的系數(shù)，最 后一個na指的是price對應(yīng)的系數(shù)。 chip.m3Call:arima(x = sales, order = c(1, 0, 4), xreg = data.frame(price), fixed = c(NA,0, NA, 0, NA, NA, NA)Coefficients:ar1 ma1 ma2 ma3 ma4 intercept price0.14440 0.26760 0.521015.8396