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

2、dow得到在911事件以前的未愛干預的時間序列子集Cm對暫用的模型進行診斷>fitmode<-arima(airmiles,order=c(0,1,1),seasonal=list(order=c(0,1,0)>tsdiag(fitmode)從診斷圖可以看出存在三個異常點，acf在12階存在高度相關因此在季節中加入MA(1)系數。擬合帶有干預信息的模型函數：arimax(x,order=c(0,0,0),seasonal=list(order=c(0,0,0),period=NA),xreg=NULL,include.mean=TRUE,transform.pars=TRUE

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

4、ir.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')&g

5、t;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)

6、,c(1,0)Coefficients:ma1I911.1-MA0sma1Dec96Jan97Dec02I911-MA0I911.1-AR1-0.3825-0.64990.0989-0.06900.0810-0.09490.8139-0.27150.11890.02280.02180.02020.04620.09780.0439sigmaA2estimatedas0.0006721:loglikelihood=219.99,aic=-423.98畫圖plot(log(airmiles),ylab="log(airmiles)")points(fitted(air.m1)2o

7、od2W4Timepl-"竹民L理匚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)I9M!1!HBJDOO30022CMMThb從上圖可以看出在2003年底后，911事件的影響效應才平息，航班客運量恢復了正常。2異常值AO

8、和10。AO和10。在時間序列中異常有兩種，可加異常和新息異常，分別記2.1異常值示例模擬數據模擬一般的ARIMA(1,0,1),然后故意將第10個觀測值變成異常值10.> set.seed(12345)y=arima.sim(model=list(ar=0.8,ma=0.5),n.start=158,n=100)> yTimeSeries:Start=1End=100Frequency=197-0.081238710.444125790.26116418-0.45815484y10<-10模型初步判斷>acf(y)stnetvw_5101520Lag>pacf(

9、y)V40s-*!1«g¥>n>eacf(y)AR/MA1 oooooooooooooooooooooooooooo2 oxoooooooooooooxoooooooooooo3 xxooooooooooooxooooooooooooo4 oxoooooooooooo從三個的結果來看，可以初步分析y是AR(1)模型對模型時行擬合m1=arima(y,order=c(1,0,0)> m1Call:arima(x=y,order=c(1,0,0)Coefficients:ar1intercept對模擬模型進行異常值探測detectAO(m1),1,2,3Iam

10、bda2-4.0184129.068982-4.247367detectAO(m1,robust=F),1ind10.000000lambda27.321709detectlO(m1)AO探測結果認為第9,10,11.可能出現異常值。10探測認為第10,11可能出現了異常值。由于檢驗統計量的最大取值出現在10且AO10,所以更認為出現異常值在第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(

11、AO=seq(y)=10)Coefficients:ar1interceptAOsigmaA2estimatedas1.059:loglikelihood=-145.29,aic=296.58detectAO(m2)1"NoAOdetected"detectIO(m2)1"NoIOdetected"比較有無異常值的兩模型再次進行異常值探測時，沒有發現異常值，驗證最初序列異常出現在10的猜測對比模型1和2的擬合效果>tsdiag(m2)雖然模型二的殘差通過引入異常值后正太性是顯性的，但是其acf和P值結果顯示引入MA（1）是必要的。重新擬合適當模型&

12、gt; 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)Call:arima(x=y,order=c(1,0,1),xreg=data.frame(AO=seq(y)=10)Coefficients:ar1ma1interceptAOsigmaA2estimatedas0.793:loglikelihood=-131.16,aic=270.33IIII

13、1*JD模型的擬合效果是顯著提高。Acf和P值檢驗也一步通過。> plot(y,type='b')arrows(40,7,11,9.8,length=0.8,angle=30)2.2另一個現實例子數據包中的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:ma1sma1sigmaA

14、2estimatedas0.5446:loglikelihood=-139.54,aic=283.08detectAO(m1.co2)1"NoAOdetected"detectIO(m1.co2),1ind57.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,

15、1,1),period=12),io=c(57)Coefficients:ma1sma110-57A2estimatedas0.4869:loglikelihood=-133.08,aic=272.16模型顯示AIC相比之前模型一更小了。而且10效應的P值=2.677/0.7246是顯著的.3偽相關在時間序列中引入協變量，如非洲牧草產量通常與某些氣候指標密切相關，在這種發問下在通過在時間序列模型中納入相關的協變量，將有助于更好的了解基礎過程以及得到更為準確的預測。3.1模擬數據set.seed(12345)X=rnorm(105)Y=zlag(X,2)+.5*rnorm(105)X=ts(X-

16、(1:5),start=1,freq=1)Y=ts(Y-(1:5),start=1,freq=1)ccf(X,Y,ylab='CCF')LegLeg3.2奶產量與對數化發電量的偽相關data(milk)data(electricity)milk.electricity=ersect(milk,log(electricity)#列合并在一個容器中。intersect函數將多個時間序1094199&T99820002GD220042D06Tinwrnilk.eiettriciryrntDMr£oFLIm-±gccf(as.numeric(mi

17、lk.electricity,1),as.numeric(milk.electricity,2),main='milk&electricity',ylab='CCF')milk&曲euEcvtyLag兩者相關性似乎非常的強，但實際上這是因為他們的各自存在很強的自相關性。4預白化與隨機回歸對于具有強自相關的數據而言，很難評估兩個過程之前是否存在依賴關系，因而，宜將x和y之間的線性關系關聯從其各自相關關系中剝離出來。預白化正是為了達到此目的的一個有效工具。4.1牛奶與電量的CCF預白化校正再次分析兩者的相關性，此時除了時滯-3具有邊緣顯著外，其他地方

18、沒有一個相關系數是顯著的。幌動防震這給出的35個樣本互相關系婁中大約會出現1.75=35x0.05個虛假警報,即這個-3系數的顯著可能就是一個虛假的信息。因此，牛奶與耗電量序列實際上是基本不相關的。從而認為之前在原始數據序列中發現的強互相關是偽相關的。4.2 Log（銷售量）與價格數據的相關性分析預白化處理ULn-blrd口2D«UMlin預白化處理prewhiten(y=diff(bluebird),1,x=diff(bluebird),2,ylab='ccf)從CCF圖可以看出兩者之間只在時滯0處是顯著的。即價格與銷售量之間存在著很強的同期負相關關系。即當期提高價格將導致

19、銷售量的當期下降。一般線性回歸分析> sales=bluebird,1price=bluebird,2> chip.m1=lm(salesprice)summary(chip.m1)Call:lm(formula=salesprice)Residuals:Min1QMedian3QMaxCoefficients:EstimateStd.ErrortvaluePr(>|t|)(Intercept)price15.8900.21773.22<2e-16*Signif.codes:00.001'*'0.01*'0.050.1Residualstanda

20、rderror:0.188on102degreesoffreedomMultipleR-squared:0.7926,AdjustedR-squared:0.7906F-statistic:389.9on1and102DF,p-value:<2.2e-16>acf(residuals(chip.m1),ci.type='ma')由于回歸后的殘差自相關在四階是顯著的，因此我們要對其進行再一步的分析r«>idua)icriip.mi|-2.4890.126-19.75<2e-16*>eacf(residuals(chip.m1)AR/MAxx

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

22、chip.m2Call:arima(x=sales,order=c(1,0,4),xreg=data.frame(price)Coefficients:ar1ma1ma2ma3ma4interceptpriceA2estimatedas0.02556:loglikelihood=42.35,aic=-70.69結果表明ma1,ma3的系數并不顯著，即可認為其系數為0調整模型>chip.m3=arima(sales,order=c(1,0,4),xreg=data.frame(price),fixed=c(NA,0,NA,0,NA,NA,NA)#第一個NA指代AR1的系數，第一個0指ma1第二個NA指的是ma2第二個0指的是ma3的系數。第三個na指ma4,倒數第二個na是指截距項對應的系數，最后一個na指的是price對應的系數。>chip.m3Call:arima(x=sales,order=c(1,0,4),xreg=data.frame(price),fixed=c(NA,0,NA,0,NA,NA,NA)Coefficients:ar1ma1ma2ma3ma4inte