8时间序列回归模型

1、时间序列回归模型1 干预分析1.1 概念及模型Box 和 Tiao 引入的干预分析提供了对于干预影响时间序列的效果进行评估的一个框架，假 设干预是可以通过时间序列的均值函数或者趋势而对过程施加影响， 干预可以自然产生也可 以人为施加的，如国家的宏观调控等。其模型可以如下表示：其中mt代表均值的变化， Nt是ARIMA过程。1.2 干预的分类阶梯响应干预脉冲响应干预1.3干预的实例分析模型初探对数化航空客运里程的干预模型的估计> data(airmiles)> acf(as.vector(diff(diff(wi ndow(log(airmiles),e nd=c(2001,8),1

2、2),lag.max=48)#用window得到在911事件以前的未爱干预的时间序列子集Cm对暂用的模型进行诊断>fitmode<-arima(airmiles,order=c(0,1,1),seas on al=list(order=c(0,1,0)> tsdiag(fitmode)从诊断图可以看出存在三个异常点，acf在12阶存在高度相关因此在季节中加入MA( 1)系数。拟合带有干预信息的模型函数：arimax(x, order = c(0, 0, 0), seas onal = list(order = c(0, 0, 0), period=NA),xreg = NUL

3、L, i nclude.mea n = TRUE, tran sform.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 函数扩展了 arima 函数，可以处理时间序列中干扰分析及异常值。假设干扰影响过 程的均值，相对未受干扰的无价值函数的偏离用一些协变量的ARM

4、A滤波器的输出这种来表示，偏差被称作传递函数。 构造传递函数的协变量通过 xtransf 参数以矩阵或者 data.frame 的形式代入 arimax 函数。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(airmile

5、s)=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(airm

6、iles) = 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: ma1I911.1-MA0sma1Dec96Jan97Dec02I911-MA0I911.1-AR1-0.3825-0.64990.0989-0.06900.0810-0.09490.8139-0.2715s.e.0.09260.11890.02280.02180.

7、02020.04620.09780.0439sigmaA2 estimated as 0.0006721: log likelihood = 219.99, aic = -423.98 画图plot(log(airmiles),ylab="log(airmiles)") points(fitted(air.m1)2ood2W4Timepl-" 竹民 L 理匚Nin e11p=1*(seq(airmiles)=69)plot(ts(Ni ne11p*(-0.0949)+filter(Nin e11p,filter=.8139,method='recursiv

8、e',side=1)*(-0.2715),frequen cy=12,start=1996),type='h',ylab='9/11 Effects')abli ne(h=0)I9M!1!HBJDOO30022CMMThb从上图可以看出在2003年底后，911事件的影响效应才平息，航班客运量恢复了正常。2异常值AO 和 10。在时间序列中异常有两种，可加异常和新息异常，分别记2.1异常值示例模拟数据模拟一般的ARIMA( 1 , 0, 1),然后故意将第10个观测值变成异常值 10.> set.seed(12345)> y=arima.sim

9、(model=list(ar=0.8,ma=0.5), n. start=158 ,n=100)> yTime Series:Start = 1End = 100Freque ncy = 11 0.49180881 -0.22323665 -0.99151270 -0.73387818 -0.67750094-1.14472133 -2.14844671 -2.495307949 -1.50355358 -2.12615253 -0.55651713 0.41326344 0.518691291.86210605 2.19935472 2.6021016517 0.79130003 0.

10、26265426 2.93414857 3.99045889 3.608226781.17845765 -0.87682948 -1.2063779925 -1.39501221 -0.18832171 1.22999827 1.46814850 2.666474913.23417469 2.60349624 1.4951321533 1.48852142 0.95739219 1.30011654 1.73444053 2.848251033.73214655 4.23579456 3.3704979041 2.02783955 1.41218929 -0.29974176 -1.58712

11、591 -1.340808780.10747609 1.44651081 1.6780948749 -0.34663129 -0.50291459 0.01739605 -0.01426474 0.942172040.39046221 -0.39883530 1.6063891857 1.70668201 1.37518194 1.91824534 0.14254056 -2.88169481-3.30372327 -1.74068408 -3.2486805765 -3.89415683 -3.45920240 -1.11042078 0.67959744 0.670510840.44394

12、061 1.89536060 2.3606387373 2.00559443 0.86443324 0.46847572 0.72338498 1.602150981.25922277 1.53180859 0.9628977981 1.07712188 1.42386354 0.56318008 -0.46689543 -0.91861106-1.92947085 -2.18188785 -1.0275908789 2.31088272 3.13847319 3.01237881 3.43454807 2.315394942.44909873 2.91589141 1.1264890897

13、-0.08123871 0.44412579 0.26116418 -0.45815484> y10<-10模型初步判断> acf(y)stnet vw _5101520Lag> pacf(y)V40 s-*!1 «g¥>n> 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 o o o o o2 o o o o o o ooo o oooo3 o x o o o o ooo o oooo4 o x o

14、o o o ooo o oooo5 x x o o o o ooo o oooo6 x o o o o o ooo o oooo7 o x o o o o ooo o oooo从三个的结果来看，可以初步分析y是AR (1)模型对模型时行拟合> m1=arima(y,order=c(1,0,0)> m1Call:arima(x = y, order = c(1, 0, 0)Coefficie nts:ar1 in tercept0.54190.7096s.e. 0.08310.3603对模拟模型进行异常值探测> detectAO(m1),1,2,3ind 9.000000 10

15、.000000 11.000000Iambda2 -4.018412 9.068982 -4.247367> detectAO(m1,robust=F),1ind 10.000000lambda2 7.321709> detectlO(m1),1 ,2ind 10.000000 11.00000lambda1 7.782013 -4.67421AO探测结果认为第 9, 10, 11.可能出现异常值。10探测认为第10, 11可能出现了异常值。 由于检验统计量的最大取值出现在10且AO 10,所以更认为出现异常值在第10是AO异常考虑异常值的时间序列拟合> m2=arima(y

16、,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)Coefficie nts:ar1 in terceptAO0.80720.5698 10.9940s.e. 0.05700.5129 0.8012sigmaA2 estimated as 1.059: log likelihood = -145.29, aic = 296.58> detectAO(m2)1 "No AO dete

17、cted"> detectIO(m2)1 "No IO detected"比较有无异常值的两模型再次进行异常值探测时，没有发现异常值，验证最初序列异常出现在10的猜测对比模型1和2的拟合效果> tsdiag(m2)虽然模型二的残差通过引入异常值后正太性是显性的，但是其acf和P值结果显示引入 MA（1）是必要的。重新拟合适当模型> m3=arima(y,order=c(1,0,1),xreg=data.frame(AO=seq(y)=10)> detectAO(m3)1 "No AOdetected"> detec

18、tIO(m3)1 "No IO detected"> tsdiag(m3)Call:arima(x = y, order = c(1, 0, 1), xreg = data.frame(AO = seq(y) = 10)Coefficie nts:ar1 ma1 in tercept AO0.6596 0.61540.5850 11.1781s.e. 0.0799 0.07960.4132 0.4755sigmaA2 estimated as 0.793: log likelihood = -131.16, aic = 270.33II1*JD模型的拟合效果是显著提高

19、。Acf和P值检验也一步通过。> plot(y,type='b')> arrows(40,7,11,9.8,le ngth=0.8,a ngle=30)2.2另一个现实例子数据包中的co2>m1.co2=arima(co2,order=c(0,1,1),seas on al=list(order=c(0,1,1),period=12)> m1.co2Call:arima(x = co2, order = c(0, 1,1), seas onal = list(order = c(0, 1, 1), period = 12)Coefficie nts:ma

20、1 sma1-0.5792 -0.8206s.e. 0.0791 0.1137sigmaA2 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),seas on al=list(order=c(0,1,1),period= 12),io

21、=c(57)> m4.co2Call:arimax(x = co2, order = c(0, 1, 1), seas onal = list(order = c(0, 1,1), period = 12),io = c(57)Coefficie nts:ma1 sma1 10-57-0.5925 -0.8274 2.6770s.e. 0.0775 0.1016 0.7246 sigmaA2 estimated as 0.4869: log likelihood = -133.08, aic = 272.16模型显示AIC相比之前模型一更小了。而且10效应的P值=2.677/0.7246

22、 是显著的.3伪相关在时间序列中引入协变量，如非洲牧草产量通常与某些气候指标密切相关，在这种发问下在通过在时间序列模型中纳入相关的协变量，将有助于更好的了解基础过程以及得到更为准确的预测。3.1模拟数据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')Leg3.2奶产量与对数化发电量的伪相关data(milk)data(electricity)milk.electricity=

23、ts.i ntersect(milk,log(electricity)# 列合并在一个容器中。intersect 函数将多个时间序1094199&T99820002GD220042D06Tinwrnilk.eiettriciryrntDMr£oFLIm-±gccf(as.nu meric(milk.electricity,1),as .nu meric(milk.electricity,2), mai n='milk & electricity',ylab='CCF')milk & 曲euEcvtyLag两者相关性似乎

24、非常的强，但实际上这是因为他们的各自存在很强的自相关性。4预白化与随机回归对于具有强自相关的数据而言，很难评估两个过程之前是否存在依赖关系，因而，宜将x和y之间的线性关系关联从其各自相关关系中剥离出来。预白化正是为了达到此目的的一个有效工具。4.1牛奶与电量的CCF预白化校正再次分析两者的相关性，此时除了时滞-3具有边缘显著外，其他地方没有一个相关系数是显著的。幌动防震 这给出的35个样本互相关系娄中大约会出现1.75=35x0.05个虚假警报,即这个-3系数的显著可能就是一个虚假的信息。因此，牛奶与耗电量序列实际上是基本不 相关的。从而认为之前在原始数据序列中发现的强互相关是伪相关的。4.2

25、 Log （销售量）与价格数据的相关性分析预白化处理ULn-blrd口2D«UMlin预白化处理prewhite n(y=diff(bluebird),1,x=diff(bluebird),2,ylab='ccf)从CCF图可以看出两者之间只在时滞 0处是显著的。即价格与销售量之间存在着很强的同期 负相关关系。即当期提高价格将导致销售量的当期下降。4.2.2 一般线性回归分析> sales=bluebird,1> price=bluebird,2> chip.m1=lm(salesprice)> summary(chip.m1)Call:lm(form

26、ula = sales price)Residuals:Min 1Q Media n3Q Max-0.54950 -0.12373 0.00667 0.13136 0.45170Coefficie nts:Estimate Std. Error t value Pr(>|t|)(In tercept) price15.8900.217 73.22 <2e-16 *Sign if. codes: 00.001 ' * ' 0.01*' 0.050.1Residual sta ndard error: 0.188 on 102 degrees of freedo

27、mMultiple R-squared: 0.7926,Adjusted R-squared: 0.7906F-statistic: 389.9 on 1 and 102 DF, 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/MA0 1 2 3 4 5 6 7 8

28、 9 10 11 12 130 x x x x o o x x o o o o o o1 x o o x o o o o o o o o o o2 x x o x o o o o o o o o o o3 x x o x o o o o o o o o o o4 o x x o o o o o o o o o o o5 x x x o x o o o o o o o o o6 x x o x x x o o o o o o o o7 x o x o o o o o o o o o o oarma( 1，4)差的回归模型。Eacf推荐其残差包含一个以(1，4)为顶点为的零值三角形，从而表明其为

29、模型，因此可将对数化销售量拟合成对于价格序列的带有ARMA( 1, 4)误:423模拟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)Coefficie nts:ar1 ma1 ma2 ma3 ma4 in tercept price0.1989 -0.0554 0.2521 0.0735 0.526915.7792 -2.4234s.e.

30、0.1843 0.1660 0.0865 0.1084 0.13760.2166 0.1247 sigmaA2 estimated as 0.02556: log likelihood = 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