8时间序列回归模型——R实现.doc

word时间序列回归模型1 干预分析1.1 概念与模型Box和Tiao引入的干预分析提供了对于干预影响时间序列的效果进展评估的一个框架，假设干预是可以通过时间序列的均值函数或者趋势而对过程施加影响，干预可以自然产生也可以人为施加的，如国家的宏观调控等其模型可以如下表示：其中mt代表均值的变化，Nt是ARIMA过程1.2 干预的分类阶梯响应干预脉冲响应干预1.3 干预的实例分析1.3.1 模型初探对数化航空客运里程的干预模型的估计> data(airmiles)> acf(as.vector(diff(diff(window(log(airmiles),end=c(2001,8)),12))),lag.max=48)#用window得到在911事件以前的未爱干预的时间序列子集对暂用的模型进展诊断>fitmode<-arima(airmiles,order=c(0,1,1),seasonal=list(order=c(0,1,0)))> tsdiag(fitmode)从诊断图可以看出存在三个异常点，acf在12阶存在高度相关因此在季节中参加MA〔1〕系数1.3.2 拟合带有干预信息的模型函数：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函数扩展了arima函数，可以处理时间序列中干扰分析与异常值。

假设干扰影响过程的均值，相对未受干扰的无价值函数的偏离用一些协变量的ARMA滤波器的输出这种来表示，偏差被称作传递函数构造传递函数的协变量通过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(airmiles)==12),Jan97=1*(seq(airmiles)==13),Dec02=1*(seq(airmiles)==84)),method='ML')> Call: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:画图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/11 Effects')abline(h=0)从上图可以看出在2003年底后，911事件的影响效应才平息，航班客运量恢复了正常。

2 异常值在时间序列中异常有两种，可加异常和新息异常，分别记AO和IO2.1 异常值示例2.1.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)> yTime Series:Start = 1 End = 100 Frequency = 1 > y[10]<-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 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 从三个的结果来看，可以初步分析y是AR〔1〕模型2.1.3 对模型时行拟合> m1=arima(y,order=c(1,0,0))> m1Call:arima(x = y, order = c(1, 0, 0))Coefficients: ar1 intercept2.1.4 对模拟模型进展异常值探测> detectAO(m1) [,1] [,2] [,3]> detectAO(m1,robust=F) [,1]> detectIO(m1) [,1] [,2]AO探测结果认为第9，10，11.可能出现异常值。

IO探测认为第10，11可能出现了异常值由于检验统计量的最大取值出现在10且AO〉IO，所以更认为出现异常值在第10是AO异常2.1.5 考虑异常值的时间序列拟合> 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 AO> detectAO(m2)[1] "No AO detected"> detectIO(m2)[1] "No IO detected"2.1.6 比拟有无异常值的两模型再次进展异常值探测时，没有发现异常值，验证最初序列异常出现在10的猜想比照模型1和2的拟合效果> tsdiag(m2)> tsdiag(m1)虽然模型二的残差通过引入异常值后正太性是显性的，但是其acf和P值结果显示引入MA〔1〕是必要的2.1.7 重新拟适宜当模型> m3=arima(y,order=c(1,0,1),xreg=data.frame(AO=seq(y)==10))> detectAO(m3)[1] "No AO detected"> 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模型的拟合效果是显著提高。

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))> Call:arima(x = co2, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))Coefficients: ma1 sma1sigma^2 estimated as 0.5446: log likelihood = -139.54, aic = > detectAO(m1.co2)[1] "No AO detected"> detectIO(m1.co2) [,1]拟合含有新息异常的模型> m4.co2=arimax(co2,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12),io=c(57))> Call:arimax(x = co2, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12), io = c(57))Coefficients: ma1 sma1 IO-57sigma^2 estimated as 0.4869: log likelihood = -133.08, aic = 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')从ccf中可。